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.

Introduction

The plant-associated microbiome is considered part of the extended plant phenotype[1,2]. Microbiota living on the surface and inside of host plants mediate processes vital to plant fitness, ranging from nutrient acquisition and stress responses, to pollination[1,3-5]. Recent quantitative characterizations of plant-associated microbiomes have focused on microbiota harbored by leaves and roots in greenhouse and common garden settings[6-11], and have identified important drivers underlying microbiota assemblies in plants aboveground (AG) and belowground (BG).

These drivers include environment-related factors such as soil source[6] and experimental site[7,8] that may vary in climatic conditions, and host-related factors including organ type[9], within-organ compartment[10], host species[11] and genotype[8]. However, there has been a lack of comparative data that broaden AG and BG microbiomes to also encompass those associated with plant reproductive organs (e.g., flowers) and sexual phenotype, and generalize relevant findings from a single to multiple host species, especially in the wild. These are, nevertheless, an important first step towards understanding how plants recruit microbiota that may in turn affect host fitness in wild populations, and for determining whether individual drivers identified under experimental settings are also of paramount importance in the wild.

While consensus from studies in experimental or managed field settings is growing that organ type (root, leaf, flower)[12] and host plant species[11] influence microbial communities, their relative roles have not been explicitly quantified. This is in part because empirical studies profiling plant AG (leaf and flower) and BG (root) microbiomes are often not directly comparable, and vary in terms of environmental settings, microbial community type and anatomical details of plant organs. For example, comparative investigations of leaf microbiomes tend to focus on epiphytic microbiota, and have often involved tens to hundreds of host species in semi-natural or natural habitats[13,14]. By contrast, those on root microbiomes have typically evaluated microbiota of a few host species[11], or different genotypes and root compartments within a single host species[6,7,10], often in experimental and/or managed field sites. Studies of flower microbiomes mostly focus on compartment differences, and/or in comparison to other organs, within a single host species in the wild or managed field sites[12,15,16]. While these studies begin to suggest a defining role of organ type in structuring bacterial microbiota in plants[1], the differences among studies make the quantitative inference of respective effects difficult. Quantifying microbiota across multiple organ types (root, leaf, flower) and host plant species in situ will provide a direct evaluation of the extent to which organ type exceeds host species in shaping microbial communities in the wild.

The effect of host species on microbiomes has been examined in controlled settings (in roots)[11] and in the wild (in roots and/or leaves)[13,17]. In wild populations, host species–microbiota association can be influenced by range differences among plant species[18]. This is because the regional species pool of microbiota that potentially colonize host plants (e.g., those harbored by soil) can vary with broad-scale abiotic environments of plant species[19]. Additionally, even when ranges overlap and plants occur in sympatry, species can vary in their local habitats, where they also alter the local species pool of microbiota[20], by enriching and depleting certain microbes. As a result, host species effect on microbiomes in the wild could be viewed as (residual) host effect, after accounting for that attributable to broad-scale abiotic environments (e.g., temperature, precipitation)[21]. Such (residual) host species effect represents species–microbiota association attributable to host phenotype (i.e., the outcome of plant genotype and local environment), and local biotic conditions (microbiota) that are modified by plant species over time; this is, by definition, the extended phenotype of host plants. Thus, we may expect host species effect on microbial communities to be stronger in the wild, owing to long-term accumulated feedbacks between plants and local microbiota, than in controlled settings. This knowledge is essential for understanding the variation of plant–microbiota association that can influence plant fitness in nature.

Microbiomes can be also influenced by host sex, as demonstrated in animal systems including humans[22,23]. Such sexual dimorphism, however, has rarely been explored in plants (but see ref.[24]), despite the fact that sexual phenotype (female, male or hermaphrodite) is known to influence floral and functional traits, and several ecophysiological processes[25,26]. As a result, plant sexes may differ in the principles governing microbiota assemblies (i.e., dispersal, habitat filtering and niche partitioning). First, the species pool of colonizing microbes can vary between female and male plants, owing to differential visitation of pollinators that may carry and disperse microbes, as an outcome of differences in floral rewards and attractive traits between sexes[27,28]. Second, plant sex can influence the niche space available for microbiota. Sexual dimorphism in flower size and longevity[26] likely affects the size and dynamics of microbial habitats. Likewise, sexual dimorphism in leaf traits (e.g., trichomes and leaf toughness)[29] can potentially define the living environments of leaf microbiota[4,30]. Finally, sex-differential susceptibility and/or allocation to defense[25,31,32] could alter resident microbial communities via microbe–microbe interactions[1]. Comparisons across AG and BG organs would inform broadly on the potential for sexual differences in microbiomes.

In this study, we aim to quantify root, leaf and flower microbiomes in three closely related wild strawberries in the broadly sympatric portion of their native ranges. These perennial, sexually polymorphic Fragaria[33] include the dioecious F. chiloensis, which is a coastal specialist growing in front dunes, and the subdioecious (hermaphrodite, male and female) F. virginiana ssp. platypetala growing in fertile mesic forest edge, as well as their natural hybrid, subdioecious F. ×ananassa ssp. cuneifolia growing in intermediate habitats[34]. All are the wild relatives of the cultivated strawberry (F. ×ananassa ssp. ananassa). Specifically, we aim to address three key questions concerning the effects of host species, organ type and sexual phenotype in these wild populations: (1) What is the relative importance of host species and organ type in structuring microbial communities? (2) Does the magnitude of the host-species effect on microbiome vary among root, leaf and flower? (3) Do microbiomes differ between host plants of different sexual phenotype?

Results

Root, leaf and flower microbiomes of Fragaria

These Fragaria plants and associated microbiota experienced distinct abiotic environments in Oregon where they are in broad sympatry (Fig. 1a,b), as revealed by the principal component analysis (PCA) of seven climatic variables[35] and elevation data of the seven populations. The first two principal components (PC1 and PC2, denoted as PC1.clim and PC2.clim) accounted for 92% of the abiotic variation. Amplicon sequencing of the V4 region of 16S small subunit ribosomal RNA (rRNA) gene, with peptide nucleic acid (PNA) clamps[36], generated substantial chloroplast and mitochondrial sequences (averagely 63% of reads per sample). After removing the plant contaminants, low-frequency OTUs and low-coverage samples, OTU observations averaged 7896 per sample (median = 2796). Individual microbial communities were normalized to the same observations for each sample using the median per-sample depth (of 2796), while keeping per-sample OTU relative abundances unchanged[37]. The resultant microbial community matrix consisted of 1577 OTUs from 27 root, 22 leaf and 23 flower samples of the three host species, representing the core set of microbiota (comprising both epiphytic and endophytic taxa) in our samples.

Figure 1

Collection and summary of Fragaria microbiomes. (a) Field collection of root, leaf and flower microbiomes of F. chiloensis (F. chilo), F. virginiana ssp. platypetala (F. virg) and F. ×ananassa ssp. cuneifolia (F.cunei) from seven wild populations (solid circles) in Oregon, USA (map generated using QGIS v2.18.10, https://www.qgis.org/; Fragaria photo credit: N. Wei). (b) PCA of climatic variables and elevation data of the sampled populations. 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.

The root, leaf and flower microbiomes contained predominantly bacterial taxa (with one Archaea OTU), among which Proteobacteria, Bacteroidetes and Actinobacteria were the most abundant phyla (Fig. 1c–e). Root microbiomes possessed more OTUs than leaf and flower microbiomes (Figs 1f and S1). The majority of leaf (94%) and flower (98%) OTUs were also found in roots (Fig. 1f).

Organ type structures plant-associated microbiota in the wild

Organ type predicted species α-diversity of microbial communities (Shannon diversity, F = 70.33, df = 2, P < 0.001; Fig. 2c; Table S1), after controlling for abiotic environments (PC1.clim and PC2.clim), host species and sex (see Methods). Shannon diversity (in terms of the least-squares mean) was highest in root microbiomes, and significantly decreased in leaf (t = 9.20, df = 51, P < 0.001) and flower microbiomes (t = 10.91, df = 51, P < 0.001; Fig. 2c). Likewise, root microbiomes harbored significantly higher phylogenetic α-diversity (Fig. 2d,e; Table S1), using the metrics[38] that both scale positively with (e.g., phylogenetic diversity, PD; Fig. 2d) and are insensitive to species richness (e.g., abundance-weighted mean phylogenetic distance, MPD; Fig. 2e). Within AG microbiomes, microbial α-diversity was comparable between leaves and flowers (Fig. 2c–e).

Figure 2

Plants harbor distinct above- and belowground microbiomes. (a) A heat map of the top 100 most abundant OTUs across host species and organ type. The colored scale bar on the left indicates (loge) OTU abundances. Hierarchical clustering was performed using the complete linkage method of Euclidean dissimilarity among microbial communities. (b) NMDS of Bray–Curtis dissimilarity revealed that microbial communities were primarily separated according to organ type rather than host species. The ellipses based on 2 s.d. are indicative of the spread of microbial communities within each organ type. (c–e) The least-squares means of Shannon diversity, log-transformed Faith’s phylogenetic diversity (logePD) and power-transformed abundance-weighted mean phylogenetic distance (MPD3) are plotted for each organ type, after controlling for all other factors. Error bars represent the 95% confidence intervals. Statistical significance is indicated: ***P ≤ 0.001; n.s., not significant.

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 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 [log2], log2FC) 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.

Depletion effects dominated AG microbiomes relative to BG microbiomes (Fig. S3). In flower microbiota, FDR-GLMs identified 22 OTUs as significantly enriched and 382 OTUs as significantly depleted (Fig. S3a). The enriched OTUs concentrated in the three dominant phyla. Thirteen out of the 22 enriched OTUs in flower microbiomes, with large effect size (typically of ≥5 log2FC), overlapped those identified by FWER-GLMs as responsive to organ type (Fig. S3a). These large-effect, enriched OTUs were primarily from Sphingomonas, Hymenobacter, Janthinobacterium, Pseudomonas, Methylobacterium and Salinibacterium (Table S3). The 382 OTUs that were significantly depleted from flower microbiomes ranged across diverse phyla (Fig. S3a; Table S3), among which 95 OTUs were also identified by FWER-GLMs as significantly influenced by organ type (Fig. S3a). These 95 OTUs had (log2) fold changes averaged −4.50 (Table S3), among which OTUs from Streptomyces, Bradyrhizobium and Steroidobacter showed the largest effect sizes (in absolute values). Flower and leaf microbiomes were similar in enriched and depleted OTUs (Figs S3, S4 and Table S3).

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.

Table 1

The marginal effects of individual variables on organ-specific compositional and phylogenetic α- and β-diversity of microbial communities.

	df	α-diversity						β-diversity
		Shannon diversity		Phylogenetic diversity (PD)		Weighted mean pairwise distance (MPD)		PERMANOVA						FWER-GLMs
								Bray–Curtis		Weighted UniFrac		Weighted betaMPD
		%Var	F	%Var	F	%Var	F	%Var	F	%Var	F	%Var	F	Dev
Flower
PC1.clim	1	4.1%	1.291	3.3%	0.643	0.2%	0.066	5.6%	1.437	12.0%	4.037**	7.2%	1.913**	484 *
PC2.clim	1	2.6%	0.809	3.4%	0.668	0.6%	0.207	5.2%	1.375	3.9%	1.559	5.1%	1.376	484 **
Species	2	5.3%	0.842	6.4%	0.627	0.8%	0.126	12.9%	1.507*	20.6%	2.920**	13.1%	1.575**	2120 *
Sex	1	16.4%	5.207*	4.8%	0.943	19.7%	5.776*	3.8%	0.884	3.5%	0.987	4.2%	1.014	1501
Species:Sex	2	24.3%	3.857*	6.1%	0.600	27.5%	4.030*	11.8%	1.192	12.9%	1.523	11.1%	1.139	1080 **
Leaf
PC1.clim	1	6.3%	1.340	0.4%	0.098	13.0%	2.772	5.6%	1.025	5.0%	0.895	4.7%	0.889	209
PC2.clim	1	0.2%	0.053	3.8%	1.013	0.0%	0.011	5.3%	0.992	5.6%	1.035	5.9%	1.101	209
Species	2	0.3%	0.033	9.6%	1.281	0.3%	0.035	10.0%	0.966	8.6%	0.827	10.4%	1.001	1771 *
Sex	1	17.6%	3.770	17.2%	4.602	14.9%	3.202	6.4%	1.231	8.3%	1.602	6.4%	1.230	1332
Species:Sex	2	19.7%	2.110	24.1%	3.217	15.6%	1.673	10.1%	0.787	10.3%	0.802	12.3%	0.978	627 *
Root
PC1.clim	1	5.4%	1.151	6.8%	1.372	0.9%	0.197	5.3%	1.187	5.2%	1.061	4.3%	0.982	1624
PC2.clim	1	1.4%	0.289	2.1%	0.415	1.8%	0.441	5.3%	1.198	3.9%	0.831	4.2%	0.954	1624
Species	2	6.6%	0.706	4.6%	0.458	5.4%	0.612	9.6%	1.131	7.7%	0.887	8.4%	0.959	4657
Sex	1	7.4%	1.575	4.5%	0.897	0.0%	0.012	3.8%	0.890	4.6%	1.052	4.4%	1.000	2079
Species:Sex	2	3.8%	0.409	2.4%	0.240	21.9%	2.486	7.1%	0.698	4.8%	0.448	9.6%	0.953	3371

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.

Figure 3

Host plant species influences aboveground but not belowground microbiomes. (a–c) OTU overlap among the three host species for root, leaf and flower microbiomes, respectively. (d–f) Constrained PCoAs of Bray–Curtis dissimilarity indicated enhanced microbial community separation by host plant species from root (d) to leaf (e) and flower (f) microbiomes, controlling for abiotic environments (PC1.clim and PC2.clim) and sex (female and male/hermaphrodite). (g–i) Constrained PCoAs of abundance-weighted UniFrac distance for root, leaf and flower microbiomes, respectively.

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 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 (log2FC 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.

Figure 4

Sexual dimorphism in relative abundances of dominant bacterial phyla in flower microbiomes. Statistical significance was assessed using proportion tests with false discovery rate control for multiple comparisons (alpha = 0.05): **P ≤ 0.01; ***P ≤ 0.001. Bars with dark color represent females; bars with light color represent males or hermaphrodites.

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.

Methods

Fragaria microbiome collection

The three Fragaria species (F. chiloensis, F. virginiana ssp. platypetala and F. ×ananassa ssp. cuneifolia) are widely distributed in western North America[50]. In Oregon, where they occur in sympatry but not in the same microhabitats[34], we collected microbiota samples from seven wild populations over a 6-day period in May 2016: two populations of F. chiloensis (Salishan, ‘SAL’: 44.919°N, 124.027°W; Strawberry Hill, ‘SH’: 44.254°N, 124.112°W); two of F. virginiana ssp. platypetala (Willamette National Forest, ‘WNF’: 44.638°N, 121.941°W; Fisherman’s Bend Recreation, ‘FBR’: 44.755°N, 122.515°W); three of the natural hybrid F. ×ananassa ssp. cuneifolia (Marys Peak, ‘MP1’: 44.497°N, 123.546°W; ‘MP2’: 44.507°N, 123.569–123.579°W; Corvallis, ‘COR’: 44.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.

DNA extraction and 16S rRNA amplicon sequencing

Samples were homogenized using a TissueLyser II (QIAGEN, Germantown, MD), and genomic DNA was extracted using Xpedition Fungal/Bacterial DNA MiniPrep (Zymo Research) under a sterile laminar flow hood. The same extraction procedure was conducted on three negative controls without plant samples. The amplification and sequencing of 16S rRNA (the V4 region) were performed at the Environmental Sample Preparation and Sequencing Facility at Argonne National Laboratory. In brief, the V4 region was amplified using 515F-806R primer pair following the Earth Microbiome Project protocol (http://www.earthmicrobiome.org/protocols-and-standards/16s/) with 12-bp barcodes. Peptide nucleic acid (PNA) clamps designed from Arabidopsis thaliana[36] were added in amplification to reduce Fragaria plastid contamination. The three negative controls failed in PCRs and generated primarily primer dimers. The 16S rRNA amplicons of the 78 samples were sequenced using a 1/5 lane of 2 × 151 bp on an Illumina MiSeq instrument.

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., loge(PD), MPD3) 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 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.

Statistical analyses of microbial community β-diversity

Compositional and phylogenetic β-diversity metrics considered Bray–Curtis dissimilarity (in vegan), inter-community MPD (betaMPD in picante) and abundance-weighted UniFrac distance in the package GUniFrac[60]. Visualization of β-diversity metrics used the nonmetric multidimensional scaling (NMDS) in vegan for Bray–Curtis dissimilarity, and principal coordinates analyses (PCoAs) by cmdscale() for UniFrac distance and betaMPD.

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 (log2 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).

Figures S1–S7

Tables S1–S5