The effects of host species and sexual dimorphism differ among root, leaf and flower microbiomes of wild strawberries in situ

Plant-associated microbiomes profoundly influence host interactions with below- and aboveground environments. Characterizing plant-associated microbiomes in experimental settings have revealed important drivers of microbiota assemblies within host species. However, it remains unclear how important these individual drivers (e.g., organ type, host species, host sexual phenotype) are in structuring the patterns of plant–microbiota association in the wild. Using 16s rRNA sequencing, we characterized root, leaf and flower microbiomes in three closely related, sexually polymorphic Fragaria species, in the broadly sympatric portion of their native ranges in Oregon, USA. Taking into account the potential influence of broad-scale abiotic environments, we found that organ type explained the largest variation of compositional and phylogenetic α- and β-diversity of bacterial communities in these wild populations, and its overall effect exceeded that of host species and host sex. Yet, the influence of host species increased from root to leaf to flower microbiomes. We detected strong sexual dimorphism in flower and leaf microbiomes, especially in host species with the most complete separation of sexes. Our results provide the first demonstration of enhanced influence of host species and sexual dimorphism from root to flower microbiomes, which may be applicable to many other plants in the wild.

Despite substantial overlap in OTUs, AG and BG microbiomes exhibited distinct community structures (Figs 2a,b and S2). Hierarchical clustering analysis (Fig. 2a) clearly separated root and leaf/flower microbiomes regardless of host species, and also showed that OTUs of high frequency in roots were not the same ones in leaves and flowers. The nonmetric multidimensional scaling (NMDS) of Bray-Curtis dissimilarity also identified AG and BG organs as the primary source of variation in compositional β-diversity of microbial communities (Fig. 2b). Similar to Bray-Curtis NMDS, the separation of root and leaf/flower microbiomes along the dominant ordination axes was revealed by principal coordinates analyses (PCoAs) of phylogenetic β-diversity metrics (abundance-weighted UniFrac distance, and inter-community MPD, betaMPD; Fig. S2).
Consistent with the qualitative inference from NMDS and PCoAs, permutational multivariate analyses of variance (PERMANOVA) showed that organ type accounted for the largest source of variation in microbial β-diversity among all tested predictors (Table S2), when using Bray-Curtis (19.4% of variation, F = 7.924, df = 2, P = 0.001) and UniFrac distance (21.1%, F = 8.838, P = 0.001). PERMANOVA of betaMPD agreed with the other two β-diversity metrics that microbial community structures were significantly affected by organ type (F = 2.130, df = 2, P = 0.001; Table S2), albeit organ type (the main effect) and its interaction with species accounted for a similar amount of variation (6.2% and 7.0%, respectively).
Supporting the PERMANOVA results, generalized linear models (GLMs) with family-wise error rate (FWER) control in mvabund 39 (denoted as FWER-GLMs) also showed that organ type strongly predicted multivariate The seven climatic variables include temperature (mean, tmean; minimum, tmin; maximum, tmax), mean dewpoint temperature (tdmean), precipitation (ppt) and vapor pressure deficit (minimum, vpdmin; maximum, vpdmin). The first two principal components visualized here (PC1 and PC2, denoted as PC1.clim and PC2.clim) were used in statistical models to control for the effects of abiotic environments. (c-e) The relative abundances of major bacterial phyla in root, leaf and flower microbiomes, respectively. Dots represent individual microbial communities; means and error bars (2 s.e.m.) are indicated. (f) OTU overlap among root, leaf and flower microbiomes of all three host species.
Scientific RepoRts | (2018) 8:5195 | DOI:10.1038/s41598-018-23518-9 abundances of microbial communities (deviance = 19,793, df = 2, P = 0.001; Table S2). In addition to organ type, host species (deviance = 19,793, P = 0.001) and sex (deviance = 12,619, P = 0.006) were also identified as having significant impact on microbial communities, unlike in PERMANOVA (Table S2). The difference between Bray-Curtis PERMANOVA and FWER-GLMs suggested that OTUs responding to host species or sex likely had low or modest variability in abundances between groups of interest, which failed to be detected by less sensitive PERMANOVA 40 , but collectively these OTUs contributed to differential microbial communities. By contrast, OTUs of large among-organ variability were likely involved in distinguishing AG and BG microbiomes, as detected by both PERMANOVA and FWER-GLMs.
FWER-GLMs that detected significant microbial community differentiation caused by organ type also identified the responsible OTUs (Table S3). These differentially abundant OTUs (N = 120) attributable to organ type accounted for a small portion (8%) of the overall OTUs ( Fig. S3; Table S3). We further assessed the effect size (i.e., fold change [log 2 ], log 2 FC) and sign (i.e., depleted/enriched) of differentially abundant OTUs in leaf or flower microbiomes relative to root microbiomes, controlling for all other factors using edgeR 41 GLMs with false discovery rate control (FDR-GLMs). As a result, 414 OTUs were identified as differentially abundant between leaf and root microbiomes, and 404 OTUs between flower and root microbiomes (Figs S3, S4 and Table S3), more than three times the OTUs identified by FWER-GLMs. This distinction was primarily caused by stringent FWER relative to FDR control.
Host species influence increases from root to leaf to flower microbiomes. As plants harbored distinct AG and BG microbiomes, host species effect (see definition in Introduction) was assessed for root, leaf and flower microbiomes separately, while controlling for abiotic environments and host sex. Host species did not predict microbial α-diversity for any organ-associated microbiomes, nor did abiotic environments (Table 1). Different from α-diversity measures, the extent to which host species overlapped in OTUs changed from BG to AG microbiomes (Fig. 3a-c). In root microbiomes, the three Fragaria species overlapped substantially in OTUs, accounting for 71-82% of individual host species OTUs (Fig. 3a); but this OTU sharing dropped to 20-58% in leaf and 27-45% in flower microbiomes (Fig. 3b,c), suggesting that microbial communities could be more similar among host species in roots, compared to leaves and flowers.
In support of this inference, constrained PCoAs of β-diversity among organ-specific microbial communities (Figs 3d-i and S5) detected diminishing community separation by host species from AG to BG microbiomes, while controlling for all other factors. For flower microbiomes (Figs 3f,i and S5), F. chiloensis (F.chilo) associated microbial communities were segregated with the other two species along the first axis of constrained PCoAs for all three β-diversity metrics, whereas those of F. virginiana ssp. platypetala (F. virg) and the hybrid F. ×ananassa ssp. cuneifolia (F.cunei) were separated by the second axis for compositional β-diversity (Bray-Curtis) metric. In leaf microbiomes, community separation was only seen when using Bray-Curtis metric between the hybrid F.cunei and two parental species along the second axis (Fig. 3e); when using phylogenetic β-diversity metrics, host species overlapped in microbial communities (Figs 3h and S5). This microbiome overlapping among host species was most pronounced in roots.
Consistent with the constrained PCoAs results, PERMANOVA (Table 1) showed that host species explained the largest source of variation in flower microbiomes among all tested predictors (Bray-Curtis, 12.9% of variation, df = 2, F = 1.507, P = 0.049; UniFrac, 20.6%, F = 2.920, P = 0.004; betaMPD, 13.1%, F = 1.575, P = 0.006). Further   Table 1. The marginal effects of individual variables on organ-specific compositional and phylogenetic αand β-diversity of microbial communities. PERMANOVA, permutational multivariate analyses of variance. FWER-GLMs, generalized linear models (GLMs) with negative binomial errors, controlling for family-wise error rate (FWER). PD, total branch lengths of OTUs within a microbial community. MPD, mean pairwise branch lengths of OTUs weighted by abundances; power transformation (with power parameter of 3) being applied to improve normality. Bray-Curtis dissimilarity, a compositional β-diversity metric. Weighted UniFrac is defined as the proportion of branch lengths of OTUs not shared between two communities, weighted by abundances. Weighted betaMPD is defined as mean branch lengths of OTUs between two communities, weighted by abundances. PC1.clim and PC2.clim, the first two axes of the PCA on seven climatic variables and elevation data. Species:Sex, the interaction term between host species and sex. %Var, proportion of variation explained; Dev: Deviance test statistic. Statistical significance is indicated: *P ≤ 0.05; **P ≤ 0.01.
evidence came from FWER-GLMs that host species strongly predicted multivariate abundances of flower microbiomes (deviance = 2120, df = 2, P = 0.030; Table 1). In leaf microbiomes, UniFrac and betaMPD PERMANOVA corroborated the inference from constrained PCoAs that host species did not affect the phylogenetic community structures of leaf microbiota (F = 0.827, df = 2, P = 0.6 and F = 1.001, P = 0.4, respectively; Table 1). Yet, the subtler microbial community separation by host species in leaves relative to flowers (Fig. 3e,f) was captured by FWER-GLMs (deviance = 1771, df = 2, P = 0.031; Table 1), albeit not by less sensitive PERMANOVA (Bray-Curtis, F = 0.966, df = 2, P = 0.5). By contrast, in root microbiomes both FWER-GLMs and Bray-Curtis PERMANOVA (Table 1) supported that host species did not predict multivariate abundances of root microbial communities (deviance = 4657, df = 2, P = 0.056 and F = 1.131, P = 0.3, respectively), as well as the phylogenetic community structures (UniFrac PERMANOVA, F = 0.887, P = 0.5; betaMPD, F = 0.959, P = 0.8). Across root, leaf and flower microbiomes, abiotic environments (PC1.clim and PC2.clim together) explained a similar amount of variation in microbial community structures as did host species (PERMANOVA, Table 1). FWER-GLMs and FDR-GLMs were used to identify the OTUs underlying microbial community differentiation caused by host species, for flower and leaf microbiomes separately. Surprisingly, no OTUs were reported as differentially abundant among host species in flowers or leaves (Tables S4, S5). However, the effect size estimation by FDR-GLMs indicated that many OTUs exhibited relatively large fold changes (log 2 FC in absolute value ≥ 5) among host species, accounting for 37% of flower and 53% of leaf OTUs (Tables S4, S5); this was likely attributable to many host species-specific OTUs and limited OTU overlapping among all three host species (Fig. 3b,c). The presence of non-significant, large-effect OTUs also suggested considerable variation in OTU abundances within host species given relatively small sample sizes, which perhaps limited the power to detecting OTU-level but not yet community-level differences among host species. Sexual dimorphism is present in flower and leaf microbiomes. Sexual phenotype predicted α-diversity of microbial communities in flowers, but not in leaves and roots (Table 1). In flower microbiomes, sex comprised the second largest source of variation in Shannon diversity (16.4% of variation, F = 5.207, df = 1, P = 0.038; Table 1) and MPD (19.7%, F = 5.776, P = 0.030), although not in PD (4.8%, F = 0.943, P = 0.3). Compared to the main sex effect, species-specific sex effect on flower microbial α-diversity was even stronger, explaining the largest source of variation (Shannon diversity, 24.3%, F = 3.857, df = 2, P = 0.045; MPD, 27.5%, F = 4.030, P = 0.040). Specifically, males harbored higher microbial α-diversity than females in F.chilo (Shannon diversity, F = 5.207, df = 1, P = 0.038; MPD, F = 5.776, P = 0.030; Fig. S6) that also has the most pronounced sexual dimorphism (see Discussion); but in the other two species, inter-sexual differences were not significant (Fig. S6).
In contrast, microbial community β-diversity was not influenced by the main effect of sex (Table 1); this pattern was consistent across organ types, β-diversity metrics and statistical models. However, species-specific sex effects on community structures were observed in flower and leaf microbiomes by FWER-GLMs (Table 1), although such interaction effects failed to be captured by less sensitive PERMANOVA.
Because small sample sizes limited the detectability of differentially abundant OTUs, we focused on phylum-level variation in relative abundances between intraspecific male/hermaphrodite and female hosts for flower and leaf microbiomes. In flowers (Fig. 4), males/hermaphrodites harbored proportionately more Bacteroidetes (in all three host species, P < 0.001 in proportion tests) and less Proteobacteria (all P < 0.001) than females, after controlling for FDR (alpha = 0.05). For flower Actinobacteria, the relative abundances were also higher in males/hermaphrodites but only in two species (F.chilo and F.cunei, P < 0.001). In leaves (Fig. S7), sex differences in the three dominant phyla were variable among host species.

Discussion
Our results support the hypothesis that organ type has a prominent role in structuring plant-associated microbiomes across host species 1,4 . In the three wild strawberries in their native environments, organ type not only predicts species and phylogenetic α-diversity of plant-associated microbiomes, but also explains the largest source of variation in compositional and phylogenetic β-diversity, while controlling for the effects of their broad-scale abiotic environments, host species and sex. In other words, the root microbiome of one host species is expected to be more similar to that of a different host species than to its own leaf and/or flower microbiome, potentially owing to niche-specific selection for adapted microbiota in different plant organs 1,4,9 . In line with earlier studies on root microbiomes 11 , Proteobacteria (45%, relative abundance), Actinobacteria (22%) and Bacteroidetes (14%) are the dominant bacterial phyla in Fragaria, with Firmicutes to a lesser extent (6%). However, these previous studies often identified substantial influence of soil type on root microbiomes in both greenhouse and manipulated field settings, stronger than the effects of host genotypes 6,42 and plant species 11 . By contrast, our study observed convergence in root microbiomes, despite distinct soil habitats where the three Fragaria species grow 34 .
The discrepancy with previous studies on the relation of root microbiomes with soil habitats likely has two reasons. First, root microbiomes quantified here comprise microbiota in association with absorptive fine roots (or first-order roots), as compared to those associated, for example, with the primary root 11,42 or the whole root system 6,10 in Arabidopsis and rice. It is possible that metabolically active root parts impose stronger filtering for colonizing microbes and thus cause stronger deviation from source soil microbiota, relative to metabolically inactive parts. Nevertheless, this hypothesis requires further investigation relating root traits to root microbiomes. In fact, several root samples in this study comprising old segments of fine roots, which were perhaps no longer metabolically active, formed a separate cluster different from the other root samples (Fig. 2a). Second, root microbiomes here did not consider low-frequency OTUs, owing to sequencing depth constraint and the use of normalized microbial community matrix and abundance-weighted diversity metrics. Although we cannot rule out the possibility that the low-frequency OTUs, which did not pass bioinformatic filtering or were not retrieved by sequencing, may indeed differ with soil habitats, their influence on microbial community structure and function remains an open question.
Fragaria leaf and flower microbiomes shared most of their OTUs with root microbiomes. This is in line with the idea that the sources of microbial assemblies in phyllosphere 4 involve colonizing microbiota from soil by rain splash, wind and the visits of ground-dwelling herbivores and pollinators, as well as endophytes migrating from root to AG organs 43 . Consistent with grapevine AG and BG microbiomes in managed field settings 12 , we also found that leaf and flower microbiomes shared proportionally more OTUs with root microbiomes than with each other, indicating soil microbiota as a common species pool for microbial assemblies associated with different plant organs. Despite substantial OTU overlap, AG and BG microbiomes differed significantly in community structures, as has been observed in leaf and root microbiomes of Arabidopsis in both controlled and wild settings 9,44 . Such AG-BG microbiota differentiation in Fragaria was attributable to many depleted bacterial taxa in phyllosphere, especially Bradyrhizobium and Steroidobacter. These two bacterial genera were also found enriched in grapevine root microbiome (relative to AG microbiomes), where they probably mediate essential processes including nitrogen fixation in roots in vineyards 12 . The enriched taxa in Fragaria leaves and flowers such as Methylobacterium, Sphingomonas and Pseudomonas have been detected in other plant species, and their abilities to withstand more hostile habitats of phyllosphere have been implicated 3,4 .
Although in wild populations plant species-microbiota association can be influenced by abiotic environments, our study showed that host species differences in microbiomes were not just the by-product of their broad-scale variation in abiotic environments, but that host species explained a substantial amount of variation of organ-specific microbial communities, after accounting for that explained by abiotic environments.
The effect of Fragaria host species was strongest in flower microbiomes followed by leaf, and weakest in root microbiomes. Such enhanced host influence from root to leaf microbiomes has been detected in other plants at both genotype 8 and species levels 17 . In common garden experiments of Boechera stricta 8 , phylogenetic community structure of leaf bacterial microbiota was clustered in accordance with host genotype but not for root microbiota. Similar to our finding of Fragaria species influence on bacterial microbiomes, host species was found as significantly affecting fungal endophytic communities in leaves but not in roots, among three common grass species in their wild populations 17 . If these patterns hold true across plant taxa, we may expect relaxed phylogenetic conservatism of plant microbiome traits (e.g., community structure) from BG to AG, owing to the paralleling phylogenetic signal of host plant functional traits.
Studies examining root and leaf traits 45,46 have revealed stronger phylogenetic conservatism in roots, with root trait variation largely among deeply diverged plant lineages and little variation among species of low taxonomic ranks. The three Fragaria congeners in this study have a relatively short divergence history (originated ~1 Mya) and identical creeping herbaceous life form 33 , and thus they likely have similar fine root traits. This perhaps explains the resemblance of root microbiomes among host species in our study. Compared to root traits, interspecific differentiation in leaf traits (e.g., coriaceousness, leaf thickness) have been seen among the three Fragaria species in the wild (Fig. 1a) and greenhouse 34 ; this may underlie the more heterogeneous leaf microbiomes relative to root microbiomes among host species. Relative to root and leaf traits, floral traits (e.g., color, size, scent and reward) are exceedingly diverse in plants 47,48 , and thus flower microbiomes are predicted to be distinct among host species 3 . Interestingly, although Fragaria flowers are morphologically similar in shape and color ( Fig. 1a; see also ref. 33 ), microbiome divergence among host species was yet strongest in flowers, suggesting that other floral traits such as scent, and pollen and nectar rewards might be critical in shaping species-specific flower microbiome. For plant genera with highly diversified floral traits, host species influence on flower microbiomes may be even stronger than what we observed in Fragaria.
Intriguingly, we found that the degree of sexual dimorphism in microbiomes coincided with the degree of sexual dimorphism in the three host species. In wild strawberries, F. chiloensis has the most complete separation of sexes than other Fragaria species, and has perhaps the greatest sexual differentiation in floral and other traits 49 . In F. chiloensis, male flowers are typically larger in petal size than female flowers 49 , which likely in part explains higher α-diversity of flower microbiomes in males. Although comparable studies on sexual dimorphism in plant microbiomes are lacking, it is noteworthy that culturable nectar-dwelling yeasts appeared to be higher in richness in male than female flowers of Silene latifolia in the wild 24 .
Relative to α-diversity, microbial community structure seems more sensitive to host sex, as many bacterial phyla were found differentially abundant between sexes in flower and leaf microbiomes across species. Sexually dimorphic leaf traits, which may underlie the observed differences in leaf microbiomes between sexes, have been detected in F. chiloensis and F. virginiana ssp. platypetala in common gardens (T-L Ashman and N Wei, unpubl. res.). For both species, males/hermaphrodites possess higher leaf nitrogen content and specific leaf area than females; but the degree of sexual dimorphism in leaf traits is still higher in F. chiloensis. Although similar data are not available for the hybrid F. ×ananassa ssp. cuneifolia, we suspect that it can be similar to F. virginiana ssp. platypetala considering their morphological resemblance 34 . Nevertheless, exactly how these sexually dimorphic traits affect microbiomes is not clear, but deserves additional research.
To conclude, our study provides the first characterization of microbiomes associated with the close wild relatives of the cultivated strawberry. We show, for the first time, enhanced host species influence and sexual dimorphism from root to flower microbiomes in wild populations. While these findings await similar investigations to generalize how plants control microbiota assemblies in the wild, it is important to recognize that such patterns of host species-microbiota association in situ affect plant interactions with AG and BG environments and plant fitness. Moreover, our results of sex-differential microbiota expand the understanding of sexual dimorphism in plants, and also highlight the needs for future research on the underlying mechanisms and on relating these differences to sex-specific fitness. Overall, findings from the wild, like ours here, strengthen those from experimental settings, and together they have broad implications for understanding this extended phenotype of plants.  506°N, 123.285°W). From each population, we randomly selected two female and two male-fertile (male or hermaphrodite) plants that were at least 2 m apart from each other, and collected root, leaf and flower samples from each plant. However, at COR we only sampled roots and leaves from three plants, as this population passed flowering. In total, our collection comprised 78 samples for the three species with organ type of root (N = 27), leaf (N = 27) and flower (N = 24).
For each plant, one flower (~2.5 cm in diameter), and the central leaflet (~2 cm in width) of a healthy trifoliate leaf that showed no evidence of herbivory and pathogen infection, were collected separately using ethanol-rinsed forceps and put directly into 750 µL Xpedition Lysis/Stabilization Solution in a ZR BashingBead Lysis tube (Zymo Research, Irvine, CA). Roots of the same plant were unearthed by the assisting person with ethanol-rinsed gloves, and were shaken vigorously to remove attached soil. Then five segments (~5 cm in length each) of fine roots, including some rhizosphere soil particles, were severed using sterile forceps and stored in the same manner. These samples were transferred to a −20 °C freezer within six hours after field collection, and shipped with dry ice to the University of Pittsburgh for DNA extraction.
Our leaf and flower samples contained both epiphytic and endophytic microbiota. The root samples also included some rhizosphere microbiota, in addition to rhizoplane and endosphere microbiota. For simplicity, we refer to these organ-associated microbiomes as root, leaf and flower microbiomes. OTU profiling and filtering. Paired-end reads were first joined using PEAR v0.9.6 51 with an overlap size of ≥20 bp. The successfully merged reads were used for subsequent open-reference operational taxonomic unit (OTU) picking. Sequence demultiplexing and quality filtering (with Phred quality scores of ≥20) were performed using QIIME v1.9.1 52 . The resulting sequences were clustered into OTUs based on a similarity threshold of ≥97% by PyNAST and assigned with taxonomic identification by RDP classifier based on the Greengenes reference database (13_8 release), as implemented in QIIME. After chimera removal using QIIME ChimeraSlayer, aligned OTU representative sequences were used to build a midpoint-rooted phylogenetic tree of these OTUs using QIIME FastTree.
The QIIME-generated OTU table was further filtered before the conversion into a microbial community matrix. First, we filtered out OTUs belonging to chloroplasts and mitochondria. Second, we removed the singletons as well as low-frequency OTUs that accounted for ≤0.01% of the total observations of the entire OTU table. Third, we removed low-depth samples of <100 observations (N = 6, five leaf and one flower samples). Fourth, we normalized the OTU table to the same observations, which were the product of the median per-sample depth and per-sample OTU proportions (or relative abundances) 37 . The resultant normalized OTU table was used as the microbial community matrix for downstream statistical analyses, because normalization using alternative per-sample depths (e.g., mean or maximum depth) and raw OTU table did not affect the results (data not shown).

Abiotic environments of sampled populations.
To account for abiotic effects on microbial communities, we used seven PRISM climatic variables 35 of the current (1981-2010) conditions at 30-arcsec resolution, and elevation data, for the seven sampled populations. The seven annual climatic variables include temperature (mean, minimum and maximum), mean dewpoint temperature, precipitation and vapor pressure deficit (minimum and maximum). We conducted a principal component analysis (PCA) of these variables, including elevation, using prcomp() in R v3.3.3 53 . The first two principal components (denoted as PC1.clim and PC2.clim) were taken as the abiotic predictors in the following statistical models.
Statistical analyses of microbial community α-diversity. Species and phylogenetic α-diversity metrics considered Shannon diversity, Faith's phylogenetic diversity (PD) and abundance-weighted mean phylogenetic distance (MPD), which were calculated using the R package vegan 54 and picante 55 . These α-diversity metrics were transformed (i.e., log e (PD), MPD 3 ) to improve normality, and used as response variables in general linear mixed models (LMMs) using the package lme4 56 . The fixed effects included PC1.clim + PC2.clim + Species + Sex + Organ + Species:Sex + Species:Organ + Sex:Organ; the random effect included individual plants. We did not include populations in random effects for two reasons: first, models that incorporated nested random effects Scientific RepoRts | (2018) 8:5195 | DOI:10.1038/s41598-018-23518-9 failed to converge given the sample size; second, individuals also captured some of the population variation. For the main effect of each predictor and their interactions, the least-squares means (LS-means) were estimated using the package lmerTest 57 , and the statistical significance was evaluated by Type III sums of squares (SS). When considering organ type separately, we subdivided the microbial community matrix by organ type and re-estimated α-diversity metrics for organ-specific microbial communities. General linear models (LMs) were fitted with PC1. clim + PC2.clim + Species + Sex + Species:Sex, in which the LS-means and Type III SS were estimated using the package phia 58 and car 59 , respectively.
These β-diversity metrics were taken as response variables in permutational multivariate analyses of variance (PERMANOVA) using vegan adonis2(). To assess the statistical significance (i.e., the marginal, instead of sequential, effect) of each main effect (or main term), PERMANOVA included PC1.clim + PC2. clim + Species + Sex + Organ. To assess the marginal effect for each interaction term, PERMANOVA included both the above main effects and their interaction terms (Species:Sex + Species:Organ + Sex:Organ). As a complement to distance-based PERMANOVA, generalized linear models (GLMs) with negative binomial errors were conducted using the package mvabund 39 to assess how community structures changed in response to the main and interaction terms. The marginal effect of each term was assessed by nested model comparison between a full model and a reduced model with the focal term removed using a likelihood ratio test. OTUs that responded significantly to each model term were identified using univariate likelihood ratio tests with P values adjusted by resampling-based multiple testing implemented in mvabund to control for the family-wise error rate (FWER; alpha = 0.05). Here we referred to mvabund GLMs as FWER-GLMs.
PERMANOVA and FWER-GLMs were also conducted to model microbial community β-diversity for each organ type separately. Visualization of organ-specific β-diversity metrics was performed using constrained PCoAs by vegan capscale().
Differentially abundant OTUs among microbial communities. As a complement to the univariate tests of individual OTU abundances in mvabund, we used the package edgeR 41 to estimate the effect size (log 2 fold change) and sign (depleted or enriched) of each differentially abundant OTU attributed to individual predictors and their interactions, as well as between different levels within a predictor. Similar to mvabund, edgeR also uses GLMs with negative binomial errors, but it models individual OTU abundances with false discovery rate control (FDR; alpha = 0.05) for multiple testing. Here we referred to edgeR GLMs as FDR-GLMs. FDR-GLMs allow a design matrix accommodating complex experimental structure. Our design matrix followed PC1.prism + PC2. prism + 0 + Group, in which Group contained all combinations of different levels of predictors (Species, Sex and Organ). The model was fitted using glmQLFit() and specific contrasts were made by glmQLFTest(). The P values were adjusted by the Benjamini-Hochberg (BH) correction using p.adjust(). We also conducted differential analyses for organ-specific microbial communities to detect differentially abundant OTUs between species within each organ type. Data accessibility. Raw reads are available from National Center for Biotechnology Information (PRJNA434446).