How to test the effect of covariate (Ex. Age) on a significant PERMANOVA test of beta diversity?

Question

[ST]

Hi,

Can anyone please help me how to analyze the effect of covariate (Ex. Age) on a significant PERMANOVA (Adonis) test of gut microbiome beta diversity? For my dataset, I ran the following scripts:

Distance:

dist <- phyloseq::distance(ps, method="unifrac")

Adonis test:

adonis(dist ~ DiseaseSeverity , data=meta, permutations=1000)

Results:

Terms added sequentially (first to last)

                                  Df SumsOfSqs  MeanSqs      F.Model      R2         Pr(>F)  
DiseaseSeverity     2    0.3698         0.18489       1.4086      0.07449   0.019 *

Residuals                35    4.5941        0.13126                         0.92551         
Total                        37    4.9639                                               1.00000         
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

So, my question is how should I formulate the script to test the effect of age on the above results? Should I do dist ~ DiseaseSeverity + AgeGroup (or dist ~ AgeGroup + DiseaseSeverity) OR dist ~ DiseaseSeverity*AgeGroup?

Thanks a bunch!