Evolutionary consequences of microbiomes for hosts: impacts on host fitness, traits, and heritability

Abstract

An organism’s phenotypes and fitness often depend on the interactive effects of its genome (⁠$G_{host}$⁠), microbiome (⁠$G_{microbe}$⁠), and environment (⁠$E$⁠). These $G$ × $G$⁠, $G$ × $E$⁠, and $G$ × $G$ × $E$ effects fundamentally shape host-microbiome (co)evolution and may be widespread, but are rarely compared within a single experiment. We collected and cultured $Lemna$$minor$ (duckweed) and its associated microbiome from 10 sites across an urban-to-rural ecotone. We factorially manipulated host genotype and microbiome in two environments (low and high zinc, an urban aquatic stressor) in an experiment with 200 treatments: 10 host genotypes × 10 microbiomes × 2 environments. Host genotype explained the most variation in $L.$$minor$ fitness and traits, while microbiome effects often depended on host genotype (⁠$G$ × $G$⁠). Microbiome composition predicted $G$ × $G$ effects: when compared in more similar microbiomes, duckweed genotypes had more similar effects on traits. Further, host fitness increased and microbes grew faster when applied microbiomes more closely matched the host’s field microbiome, suggesting some local adaptation between hosts and microbiota. Finally, selection on and heritability of host traits shifted across microbiomes and zinc exposure. Thus, we found that microbiomes impact host fitness, trait expression, and heritability, with implications for host–microbiome evolution and microbiome breeding.

duckweed, species interactions, adaptation, G × G, G × E, urban pollution

Introduction

The microbial communities making a living on larger hosts are an integral part of their host’s ecology and evolution, affecting trait expression and fitness in hosts across the tree of life (Dittami et al., 2016; Friesen et al., 2011; Gould et al., 2018; Honegger, 1993; Turnbaugh et al., 2008). In turn, host trait variation can affect microbial survival, growth, or transmission (Adler et al., 2021; Buffington et al., 2016). Theoretical models suggest that microbiomes will contain genetic variation that impacts host traits and that this variation will respond to selection on hosts (Henry et al., 2021; Mueller and Linksvayer, 2022; O’Brien et al., 2021; Rebolleda-Gómez et al., 2019). Emerging empirical results support this idea: altering microbiomes inoculated onto hosts can change the additive genetic variance and covariance of traits (O’Brien et al., 2019), and both host traits and host fitness can “evolve” in microbiomes when host genotypes are held constant (Batstone et al., 2020; Lau and Lennon, 2012; Panke-Buisse et al., 2015; Tso et al., 2018). These and similar results have spurred continued calls to engineer microbiomes for agricultural, medical, or other applications (e.g., Epstein et al., 2019; French et al., 2021; Khan et al., 2021).

We expect that much “evolution” of host traits in microbial genomes or communities results from links between microbial and host fitness via reciprocal fitness feedbacks (Sachs et al., 2004). For example, in plants hosting nitrogen-fixing microbes, microbe-derived nitrogen may allow the plant host to grow larger, which in turn may translate into more plant-fixed sugars provided to microbes, and thus a positive fitness feedback between partners (Friesen, 2012). Variation in host fitness caused by expressed host trait values could therefore impact the relative fitness of microbial genotypes or species that have variable effects on those traits. Indeed “evolved” effects of microbes on host phenotypes have been linked to both community composition turnover and genetic changes within strains (Batstone et al., 2020; Lau and Lennon, 2012; Panke-Buisse et al., 2015; Tso et al., 2018), meaning apparent “evolution” in host traits was in fact due to ecological or evolutionary changes in microbes. Relatedly, microbiomes can alter selection on host traits, for example by changing the strength or direction of selection on flowering time (Chaney and Baucom, 2020; Fitzpatrick et al., 2019; Lau and Lennon, 2011; Wagner et al., 2014). As a result, when hosts and microbes share an evolutionary history, ongoing fitness feedbacks may cause local adaptation between them (Houwenhuyse et al., 2021; Rúa et al., 2016).

The dependence of expressed trait values, heritable variation for traits, and nature of selection on microbiome context are not unique to the microbiome environment. Like microbiomes, physical environments can simultaneously change trait expression, the amount of heritable genetic variation, and selection pressures (Wood and Brodie III, 2015, 2016). Furthermore, physical environments may alter the evolutionary forces that act on the variation for host traits and fitness contained in microbiomes (Bertness and Callaway, 1994; Bronstein, 1994; O’Brien et al., 2018). For example, nutrient pollution disrupts positive fitness feedbacks between hosts and microbial symbionts (Shantz et al., 2016), and evolution in microbes in response to nutrient loading can lead to fitness reductions in hosts (Weese et al., 2015).

Before setting out to engineer or evolve host traits in microbiomes, we must ask how much trait variation actually depends on the microbiome versus host genome, and whether microbiome effects are predictable or contingent on host genotypes or physical environments. Previous studies (e.g., Grieneisen et al., 2021) have quantified what is often termed “microbiome heritability,” or the proportion of variation in microbiome composition explained by host genotype. But for microbiome breeding, the converse matters more: we need to measure the “heritability” of host traits within the microbiome, or the proportion of variance in a “host” trait that can be attributed to microbial genes (Henry et al., 2021; Mueller and Linksvayer, 2022; O’Brien et al., 2021). Despite a wealth of knowledge on effects of microbes on host traits, host fitness, and even host fitness across environments (e.g., Friesen et al., 2011; Thrall et al., 2008; Turnbaugh et al., 2008), few studies explicitly consider microbial effects on the trait variation available to selection (e.g., Chaney and Baucom, 2020; Epstein et al., 2023; Fitzpatrick et al., 2019; Heath et al., 2012; O’Brien et al., 2019; Wagner et al., 2014; Wendlandt et al., 2021). The common duckweed Lemna minor and its manipulable, largely mutualistic microbiome (O’Brien et al., 2020, but see Jewell et al., 2023) is a useful model for such experiments because of L. minor’s small size and fast generation time. Here, we carried out a full-factorial $G\times G\times E$ experiment using duckweeds and microbes to directly compare the relative importance of host genotype (⁠$G_{host}$⁠), microbiomes (⁠$G_{microbe}$⁠), physical environment (⁠$E$⁠), and their interactive effects (⁠$G_{host}\times G_{microbe}$⁠, $G_{host}\times E$⁠, $G_{microbe}\times E$⁠, and $G_{host}\times G_{microbe}\times E$⁠) to host phenotypes and fitness.

Duckweed is a floating aquatic plant with a wide distribution across diverse environments (RDSC, 2016) that may alter duckweed interactions with their microbiomes. Urban habitats may be one such environmental factor. Urban-to-rural gradients can be stark ecotones, drive a substantial amount of evolution (Johnson and Munshi-South, 2017), and have been observed to disrupt interactions between a host plant and root-nodulating microbes (Murray-Stoker and Johnson, 2021). Another important environmental characteristic may be zinc exposure. Zinc is a duckweed stressor (Dirilgen and Inel, 1994; Jayasri and Suthindhiran, 2017; O’Brien et al., 2020; Radić et al., 2010) and may vary across duckweed habitats. It is released from car tires and building roofs in urban areas (Göbel et al., 2007), but can also be high in rural areas due to natural geologic variation or use in agricultural pesticides and fertilizers (Ontario Ministry of the Environment, 2011; Zhang et al., 2003). We tested all possible combinations of 10 clonal L. minor lines from an urban-to-rural gradient and the culturable portion of their field-collected microbiomes in two environments: high and low zinc.

We hypothesized (a) that duckweeds and microbes benefit from shared ecological and evolutionary history (i.e., that local adaptation between duckweeds and microbes underlies some $G_{host}\times G_{microbe}$ variation), and that contamination context might expose $G_{host}\times G_{microbe}$ effects (i.e., cause $G_{host}\times G_{microbe}\times E$ effects). As previously observed in other systems (Wagner et al., 2014; Wood and Brodie III, 2016), we also expected that (b) microbiome and environmental context (zinc contamination) would alter selection on host phenotypes. More generally, $G_{host}\times G_{microbiome}$ effects may be detectable only when the difference in microbiomes is sufficiently large, and more diverged microbiomes may be more likely to produce nonadditive effects, in much the same way that more diverged genomes cause more epistasis (Dettman et al., 2007; De Visser et al., 2011). Therefore, we also predicted that (c) more dissimilar microbiomes are more likely to affect duckweed hosts in genotype-specific ways, that is, that $G_{host}\times G_{microbiome}$ effects will be larger when comparing more dissimilar microbial communities.

Methods

Biological materials

In May–August 2017, we collected common duckweed, Lemna minor,¹ and its associated microbes from 10 sites in the Greater Toronto Area (see Supplementary Table S1), spanning a gradient from urban to suburban to rural. From each site, we isolated microbes from 1 to 2 duckweed fronds by pulverizing freshly collected plant tissue (with adhering liquid) and streaking onto an agar plate with yeast-mannitol (YMA) media, which we incubated at 29 ∘C for 5 days before storing at 4 until use in the experiment. These microbes thus comprise the subset of L. minor epiphytic and endophytic microbes that can be maintained in mixed culture on YMA media. We also collected approximately 50 additional L. minor fronds from each site which we froze for 16S rRNA gene profiling of field L. minor microbiomes, where the sample again included plant tissue and adhering liquid. Finally, we isolated one living L. minor frond from each site. We maintained the progeny resulting from individual field fronds in growth media (Krazčič et al., 1995, refreshed every $\approx$2 weeks), in vented glass 500-ml mason jars, in a growth chamber (ENCONAIR AC80, Winnipeg, Canada) with a cycle of 16 hr at 23 ∘C at 150 $\text{μ}$mol/m${}^2$ and 8 hr of 18 ∘C in the dark. Fronds grew from a single individual to hundreds in about 3 months and were maintained at high density until the experiment.

In the field, L. minor produces new fronds mainly clonally (i.e., by budding), but can flower and reproduce sexually via seed. However, genetic evidence shows that this occurs rarely (we never observed flowers in our cultures) and that there is little segregating diversity within populations, meaning that even sexually produced isoparental offspring may be nearly isogenic (Ho, 2017). Furthermore, duckweeds from ponds only a few kilometers apart in this region have been observed to differ in genotype (Ho, 2017). Therefore, we consider our L. minor cultures derived from a single field-collected frond as isogenic genotypes.

Sequencing of microbial communities

To characterize and compare the microbiomes of field-collected L. minor and the cultured communities we used to inoculate experiments, we conducted 16S rRNA amplicon sequencing. We extracted DNA from field-collected frozen duckweed tissue with DNeasy Powersoil kits (Qiagen, Hilden, Germany), using the recommended fresh tissue weight (0.25 g). We extracted DNA from our cultured microbiomes with GenElute Bacterial Genomic DNA kits (Millipore-Sigma, St. Louis, MO, USA). We sent $\approx$10 ng of DNA from each extraction to Genome Québec (McGill University) for PCR amplification of the 16s rRNA gene (V3–V4 region, 341f/805r, Supplementary Table S2, a region expected to be optimal for identifying microbiome composition, Mizrahi-Man et al., 2013; Thijs et al., 2017), normalization, and barcoded, paired-end, 250 base pair format sequencing on an Illumina MiSeq System (San Diego, CA, USA). We received $\approx$962 million base pairs in $\approx$1.9 million demultiplexed reads, ranging from 26,029 to 129,507 reads across samples.

We used QIIME2 software to process reads (Bolyen et al., 2019, version 2020.8.0). We trimmed adapters and primers with “cutadapt” (Martin, 2011), and assigned amplicon sequence variants (ASVs) to paired reads and filtered sequencing errors with tool “DADA2” (Callahan et al., 2016). We then inferred taxonomic affiliation of ASVs using the naive Bayesian “classify-sklearn” method trained on the Greengenes database (Bokulich et al., 2018; McDonald et al., 2012; Pedregosa et al., 2011). Because processing decisions can affect downstream analyses (Bardenhorst et al., 2022; Prodan et al., 2020), we also corrected sequencing errors with “deblur” (Amir et al., 2017) at two different stringency settings: one similar to our DADA2 settings (singleton reads in any sample must occur in $>$1 sample to potentially be retained as ASVs), and one more restrictive (singleton reads in any sample must occur in $\ge$10 samples to be potentially retained as ASVs). When using deblur, reads must be joined (Rognes et al., 2016), quality filtered (Bokulich et al., 2013), and set to equal length (402 bp, the length beyond which most reads did not extend and after which quality declined, identified with QIIME2’s visualization tools). We limited recording taxonomy identifications to those inferred with a confidence of 70% or higher, and we filtered out reads identified as chloroplast (Streptophyta only) or mitochrondria. We inserted reads into the Greengenes reference phylogeny (Janssen et al., 2018), discarding those that could not beplaced.

After all processing steps and using the “denoise” method, 1,026,828 reads remained, ranging from 19,239 to 84,116 across samples (mean 51,341). Field samples had a similar number of reads per sample as cultured samples pre-processing, but fewer postprocessing (average 32,024 for field samples versus 67,995 for cultured samples postprocessing). See Supplementary Table S3 for reads remaining when using deblur methods. Field samples were extracted from plant tissue, and contained numerous Streptophyta chloroplast reads (average of 29,240 per sample), accounting for the postprocessing difference in depth. Sample coverage (the estimated proportion of total sample DNA sequences belonging to sequenced ASVs) ranged from 99.9% to 100%, indicating that missed ASVs likely comprised only a small fraction of DNA samples (calculated with package iNEXT in R, Chao et al., 2014; Hsieh et al., 2020; R Core Team, 2022). We note that while sample coverage estimates are expected to be fairly robust to processing methods, error is an unavoidable part of sequence data analysis (Bardenhorst et al., 2022), and this metric evaluates only whether the sequences in the sample have been thoroughly characterized. It does not evaluate whether the sample contained sufficient DNA to characterize the source. We compared communities with packages vegan and MCMCglmm in R (Hadfield, 2010; Oksanen et al., 2022). Data formatting and plotting used packages qiime2R, tidyverse, and Polychrome (Bisanz, 2018; Coombes et al., 2019; Wickham et al., 2019).

We constructed pairwise distance matrices of each sequenced microbial community to the other using abundance weighted UniFrac distance (Lozupone and Knight, 2005, a measure of the difference between communities in phylogenetic space, hereafter “UniFrac distance”). Because UniFrac distance captures changes in abundance across taxa; the pairwise difference between sequenced communities in relative abundance of various bacterial families correlates positively with UniFrac distance (Supplementary Figure S1). We plotted UniFrac distances, taxonomic breakdown, and ASV overlap across paired fields and cultured community samples in R (R Core Team, 2022).

Experimental design and data collection

We performed full factorial manipulation of duckweeds from 10 sites, duckweed microbial communities from those same sites, and the presence or absence of an environmental stressor: high concentration of zinc.

We set the high zinc treatment to 3.44 $\text{μ}$M (approximate level in runoff from an urban highway or waste discharge, Göbel et al., 2007; Liskco and Struger, 1996; Ontario Ministry of the Environment, 2011), by adding ZnSO${}_4$ to growth media. We used the unaltered media for the low zinc treatment, in which zinc was present at 0.86 $\text{μ}$M because it is a micronutrient. Zinc is also naturally present at low levels in duckweed habitat (Liskco and Struger, 1996).

To prepare microbial inocula, we stirred a swab from the stored agar plate for each site into separate tubes of liquid YMA, which we then cultured in a shaking incubator (VWR, Radnor, PA, USA) for 3 days at 30 ∘C and 200 rpm. Using lab-cultured inocula increases the relevancy to engineered microbiome applications, which rely on synthetic microbial communities (Kaminsky et al., 2019). However, these methods introduce taxonomic biases and stochastic assembly differences (Jones et al., 2022; Rappé and Giovannoni, 2003). We inoculated to obtain 2,000 cells/$\text{μ}$L, a density within the range we have observed at duckweed field sites (⁠$\approx$850–6400 cells/$\text{μ}$L; C. Carlson & J.R. Laurich, unpublished data). We diluted inocula using a calibration curve for optical density to cell density. Our optical density transformation is predictive of colony-forming units in mixed communities (O’Brien et al., 2020), but it is imperfect: the relationship between optical density and number of live cells can vary both across and within microbial species (Volkmer and Heinemann, 2011).

We replicated the 200 experimental treatments (10 duckweed source sites $\times$ 10 microbe source sites $\times$ 2 zinc treatments) 8 times, and randomly arranged these 1,600 experimental units across 67 twenty-four-well plates. We added 2.5 mL of sterile high- or low-zinc media, one mother-daughter frond cluster of surface-sterilized duckweeds (dipped in 0.5% bleach for 30 s), and 80 $\text{μ}$L of microbial inocula to each well. We sealed plates with BreatheEasy membranes (Millipore-Sigma, Diversified Biotech, Dedham, MA, USA) to prevent microbial contamination and cross-contamination, while still allowing gas exchange. We placed plates in a growth chamber set to the same conditions as for the duckweed cultures (above), and let them grow for 14 days. All microbial or sterile manipulations were performed in a biological safety cabinet (ESCO Micro Pte. Ltd., Labculture, Singapore). We autoclaved all media solutions prior to use.

We photographed plants with a custom-built camera rig (Nikon D3200 with AF-S DX NIKKOR 18-55mm f/3.5-5.6G VR lens, Minato, Tokyo, Japan, Yongnuo YN-300 light, Shenzhen, China) on days 0 and 14, removing plate seals on day 14 for accurate color. We used ImageJ (Schneider et al., 2012, version 1.51) to measure plant growth and traits. As a metric of fitness, we calculated the mm${}^2$ area of live duckweeds in each well at the end of the experiment, which measures duckweed growth from the single added pair, and which is tightly correlated to the number of fronds (each frond is one individual, O’Brien et al., 2020). We averaged three measures across duckweeds in each well: color intensity, aggregation (the ratio of frond area to frond perimeter, see O’Brien et al., 2019), and roundness of separate duckweed outlines (⁠$4\times A/\left(\pi \times d^2\right)$⁠, where $d$ is the major axis length, and $A$ is the area). Aggregation becomes high when daughter fronds are retained on the mother (i.e., do not abscise), a process that varies across duckweed genotypes and in response to stress (Henke et al., 2011; Newton, 1977; O’Brien et al., 2019, 2020). Roundness may capture elements of both aggregation, as fronds that have not abscised have a less round outline, and frond expansion, which duckweed may also vary in response to environmental conditions (Paolacci et al., 2018). Color intensity quantifies darker, or less white pixels in images, and high values may indicate thicker fronds, capturing another axis of frond expansion.

After photographing on day 14, we froze and stored plates at -20 ∘C. We then extracted and defrosted plates individually to measure optical density at 600-nm wavelength (OD600) in 70 $\text{μ}$L from each well using a spectrophotometer (BioTek Synergy HT with Gen5 1.10 software, Winooski, VT, USA). Because freezing alters optical density values, albeit consistently across samples (J.R. Laurich, unpublished data), we do not report results in absolute cell densities. Instead, we use OD600 as a measure of prefreezing total microbial cell density, suspended in liquid at the end of the experiment, including any dead cells. Because wells started at the same cell densities, OD600 captures summed growth across community members (O’Brien et al., 2020), but may unequally sample microbes with different growth habits, for example, those that form biofilms versus swim freely, which may vary among freshwater bacteria isolated from duckweed (Pernthaler, 2017).

Data analysis

Estimating plant genotype, microbiome, and environment effects

We calculated the percent of the variation in duckweed growth explained by each treatment as the treatment sums of squares (the summed squared distances between each point and its respective treatment mean) divided by the total sums of squares (the squared distance between each point and the grand mean). We abbreviate this %SS. For interactive effects, we calculated %SS with means for all unique combinations of the interacting treatments (e.g., each unique host–microbe combination), and subtracted the %SS for each treatment when it was considered alone. Thus, any variation that could be explained by main or interactive effects was attributed only to the main effects. We also calculated these percentages for duckweed frond traits and microbial density.

For each duckweed trait we measured, we calculated its “heritability” in both hosts and microbiomes. We calculated broad-sense heritabilities of traits in duckweeds the standard way: by dividing the phenotypic variance among duckweed genotype means by the total phenotypic variance (i.e., $V_{G_{host}}$/$V_{total}$⁠, Falconer, 1996). Our duckweed lines grow clonally, so broad-sense heritability (including additive genetic, epistatic, dominance, and maternal effects) most accurately describes the heritable variation visible to selection. We calculated duckweed broad-sense heritability for each trait in each unique zinc and microbial treatment combination.

We also calculated the “heritability” of duckweed traits in the microbiome, or the proportion of total phenotypic variance caused by microbiome genetic variation across hosts (i.e., $V_{G_{microbe}}$/$V_{total}$⁠, Henry et al., 2021; Mueller and Linksvayer, 2022), again separately for each unique zinc and duckweed genotype treatment. Much as the broad-sense heritability of duckweed traits in duckweeds (i.e., $V_{G_{host}}$/$V_{total}$⁠) includes epistatic, dominance, and maternal effects beyond those of additive genetic variation, the “heritability” of duckweed traits in the microbiome (⁠$V_{G_{microbe}}$/$V_{total}$⁠) includes both the additive genetic effects of microbial alleles on duckweed traits and additional effects such as those arising via epistasis between loci of distinct microbial species in the microbiome. Furthermore, while duckweed alleles are predictably transmitted across generations, how faithfully microbial alleles are transmitted from mother to daughter duckweed fronds is unknown under natural or un-manipulated conditions. In our experiment (and, we would argue, most microbiome breeding contexts), however, we enforced transmission by inoculating duckweeds with the same microbes, although there is likely stochasticity in how these microbes subsequently colonize hosts. This makes $V_{G_{microbe}}$/$V_{total}$ somewhat analogous to broad-sense heritability, in the sense that this quantity represents the variation visible to artificial selection during microbiome breeding, and thus we refer to $V_{G_{microbe}}$/$V_{total}$ as the “heritability” of the duckweed trait in the microbiome. Note, however, that this is distinct from the heritability of microbiome composition in hosts, often termed “microbiome heritability,” as in Grieneisen et al. (2021).

We permuted the trait values across our treatments in 10,000 randomizations and recalculated each “heritability” value. We then compared observed “heritability” values to the 95% highest posterior density intervals (HPDIs, Bayesian equivalents of confidence intervals) for the “heritabilities” calculated on the permuted data to determine a detection threshold, below which there was no statistical significance for heritability estimates.

Before quantifying how environmental gradients describing duckweed source, and gradients describing microbe source, affected duckweed fitness and expressed trait values, we first accounted for strong spatial effects. We identified the best model of spatial covariates for each response variable (x and y growth chamber coordinates, row and column plate coordinates, including both linear and squared terms, $f\left(Space\right)$⁠). Linear terms can account for gradients across plates or the chamber, while squared terms can account for edge effects, both of which could alter light or airflow. We removed nonsignificant spatial covariates one at a time and re-fit iteratively until DIC stopped improving (Spiegelhalter et al., 2002, deviance information criterion, evaluates model fit).

Because we $a$$priori$ expected that urban gradients might structure biotic variation in duckweed and microbes, we did not fit unique means for each source, and instead modeled $G_{host}$ and $G_{microbe}$ effects as being structured by this gradient. We characterized the urban-to-rural gradient by estimating the percent of permeable surfaces within a 0.5 km radius (⁠$PS$⁠, using Google Earth 2018 images and ImageJ). Low permeable surface area indicates a high prevalence of roof, road, parking lot, and other developed surfaces. We expect this to be the most relevant metric of “urbanness” for duckweeds, which are exposed to the hydrologic and chemical changes associated with decreases in permeable surface area in urban environments (Masoner et al., 2019; Miles et al., 2021; Walsh et al., 2005). Permeable surface area has the advantage of being quantitative, rather than categorical. We considered two- and three-way interactions between host, microbe, and environment effects, such that the full model was: $\begin{array}{c}Y\sim \alpha +P{S}_{host}+P{S}_{microbe}+Zinc+P{S}_{host}\times Zinc+P{S}_{microbe}\times Zinc+P{S}_{host}\times P{S}_{microbe}+\\ P{S}_{host}\times P{S}_{microbe}\times Zinc+f\left(Space\right)+ϵ\end{array}$⁠, where Y is duckweed frond area, frond color intensity, frond aggregation, frond roundness, or microbial density, and $\alpha$ is the intercept.

We fit these gradient models in R with MCMCglmm (Hadfield, 2010; R Core Team, 2022), using the same methods for these and subsequent linear models. We evaluated the significance of fixed effects by whether the HPDIs overlapped 0 on either side of the interval at pMCMC $<$ 0.05 (Bayesian equivalent of p-values), and if not, whether removing the term worsened model fit. Marginally significant terms (pMCMC $<$0.1, but $>$0.05) were retained in models if they improved fit, but were not considered “significant.” We removed non-significant terms one at a time and re-fit iteratively. All variables were fit with a Gaussian distribution, though we first took the inverse of microbial optical density (Shapiro–Wilk W statistics for transformed optical density and all other variables untransformed ranged 0.95–0.99, where 1 indicates a perfect normal distribution). We used 11,000 MCMC iterations, 1,000 burn-in, and thinning of 10 for model selection, and refit best models with 101,000 iterations.

Functional differences in microbiome effects on hosts might be driven by differences in species or strain composition across cultured communities. Species composition differences can be evaluated with 16s rRNA relative abundance data. To investigate the effects of species composition while avoiding tests for every single taxonomic group, we tested only whether the relative abundance of the two most abundant families of bacteria in cultured communities correlated with effects on duckweeds. We substituted each relative abundance into the above model for $P{S}_{microbe}$⁠, and compared DIC from the best model with the relative abundance to the best model with $P{S}_{microbe}$⁠.

Testing benefits of shared ecological or evolutionary history

We hypothesized that shared history of duckweeds and microbes with each other should have predictable effects. Local adaptation between duckweeds and microbes would be a source of $G_{host}\times G_{microbe}$ effects, where matched hosts and microbes result in higher duckweed frond area or microbial density. If any local adaptation between duckweeds and microbes depends on zinc contamination context, this would be a source of $G_{host}\times G_{microbe}\times E$ effects.

While we can test for local adaptation in this system, we cannot partition any observed local adaptation between the duckweeds and microbes, or among the microbe species. Duckweeds could have evolved to perform best in their local microbial communities or to support more growth of specific local microbes. Likewise, one or more microbes could have evolved to grow best on their local host genotypes, or to better support the growth of their local host genotypes. Further, if only one or a few microbes benefit from local adaptation (grow most on the local host), we will detect this in our microbial density metric only if these microbes contribute a large fraction of the total microbial cells. Our test for benefits to microbes from local adaptation is therefore conservative.

We tested whether there was any evidence of local adaptation between duckweeds and microbiomes (sympatric effects), and whether this was altered by zinc using linear models. We included growth chamber spatial covariates as above. We also included random effects for duckweed line and microbe community (⁠$N_D\left(0\right)$ and $N_M\left(0\right)$⁠, respectively) because genotypes that are always more fit or environments (microbes) that are always better can bias tests of local adaptation (Blanquart et al., 2013). Random effects are assumed to be drawn from normal distributions. To test our hypothesis, we fit the model $Y\sim \alpha +Sympatry+Zinc+Sympatry\times Zinc+N_D\left(0\right)+N_M\left(0\right)+f\left(Space\right)+ϵ$⁠, where $Y$ is either duckweed frond area or microbial density (inverse-transformed optical density at 600 nm).

Because microbiomes vary in composition, non-local microbiomes may differ more or less strongly from local microbiomes. A common finding of studies of local adaptation is that genotypes often perform well in all environments that are close to local conditions (Montalvo and Ellstrand, 2001; Wang et al., 2010; Wilczek et al., 2014). In the same manner, hosts may perform better in microbial communities that are quantitatively more similar to those with which a host shares evolutionary history. Collected microbiomes from the local site may vary in how accurately they represent the local field microbiome, due to the stochastic effects in culturing or assembly (Jones et al., 2022; Kaminsky et al., 2019). For applications of microbiome engineering, it is an open question how closely synthetic microbiomes need to match natural microbiomes. Therefore, we also tested whether duckweed fitness (frond area) or microbial density (inverse-transformed optical density at 600 nm) increased when inoculated microbiomes were more similar in composition to microbiomes observed on duckweed genotypes in the field (pairwise UniFrac distance, $PairwiseD$⁠). Specifically, we fit the model $Y\sim \alpha +PairwiseD+Zinc+PairwiseD\times Zinc+f\left(Space\right)+ϵ$⁠, where Y is either duckweed frond area or microbial density.

Testing whether microbiomes and zinc alter selection on duckweed traits

We tested whether there was altered selective pressure on measured duckweed phenotypes in our different zinc and microbial community contexts by fitting models between genotype means for our measure of duckweed fitness (frond area) and genotype means for each trait, as well as between genotype means for frond area and microbial density. We fit linear and univariate selection gradients (⁠$Y\sim \alpha +\beta \times Trait$⁠; Wood and Brodie, 2016) separately for each unique zinc-by-microbial community treatment (10 duckweed genotype means calculated from 80 datapoints per model). We permuted our data set for each unique treatment 1,000 times (randomizing trait and fitness values across each set of 80 points), repeated the analysis, and recorded DIC as a measure of model fit (DIC is comparable because models retain the same y-values). This produced a null distribution for model fit, and we considered a univariate selection gradient significant if the fit was better (DIC was lower) than 95% of the models in permutations.

Note that we used untransformed microbe density (optical density at 600 nm), as the inverse transformation is not appropriate for selection analyses (Hansen and Houle, 2008). We scaled fitness at the full data set level, since this is the appropriate level to scale at when comparing fitness across demes (De Lisle and Svensson, 2017, our duckweed lines are best considered as different populations). The distribution of duckweed trait values was not constant across microbiome treatments (see Results, Figure 2), which means scaling at the full data set level is not appropriate (De Lisle and Svensson, 2017). Therefore, we scaled traits within each microbiome treatment independently.

Testing whether more dissimilar microbiomes increase $G_{host}\times G_{microbe}$

We hypothesized that the effects of inoculated microbiomes might depend on how similar microbiomes were; specifically, that we are more likely to observe host-dependent effects of microbiomes when comparing more dissimilar microbiomes. If host-dependent effects of microbiomes do not exist, then the correlation of trait values expressed by our set of duckweed genotypes in one microbiome to values expressed by the duckweed genotypes in another microbiome should be $\approx$ 1. Conversely, if host-dependent effects of microbiomes are strong, the correlation between trait values expressed by duckweed genotypes in one microbiome and values expressed by duckweed genotypes in another microbiome will be substantially reduced. Correlations could theoretically be as low as -1, which would indicate that the relative trait values of the duckweed genotypes were exactly reversed across the two compared microbiomes.

We calculated the correlation across duckweed genotype means between microbiomes for all pairwise microbiome combinations, and we call this correlation metric “trait similarity.” We measured microbiome similarity as the pairwise UniFrac distance between cultured microbiomes. We then tested whether similar microbiomes produce more similar trait effects with a linear model (⁠$TraitSim\sim \alpha +PairwiseD+ϵ$⁠).

Results

Cultured microbiomes are less diverse than field microbiomes but share taxa

Microbial communities from the field are substantially more diverse than cultured communities: amplicon sequence variants (ASVs) totaled 3,949 across field samples versus 232 across cultured communities. Culturing communities also biased the taxonomic composition of experimental microbiomes. Some families that were common in sequence data from field-collected duckweed, such as Rhodobacteraceae and Pirellulaceae (1–28% and 1–30% of field samples), were very rare in cultured communities (detected in 1 community at less than 0.05%). Cultured communities were instead dominated by Aeromonadaceae and Pseudomonadaceae (Figure 1A), which were sequenced at low abundance in most field communities (9 and 10 field samples, at a max of 0.5% and 0.7%, respectively). We note field communities have a high proportion of unidentified ASVs, which could affect the accuracy of these proportions.

Figure 1

(A) Relative abundance of 16S rRNA gene sequences identified to family. Bars are ordered by field (left 10 bars) or inoculum (right 10 bars) collection site and by percent permeable surface area within a 0.5 km radius (increasing from left to right). “Unidentified” refers to the sum of all groups of sequences that could not be identified to family level. “Other” refers to the sum of all groups that were individually $<$5% of the total sequences. (B) Proportion of sequences in cultured communities identical to sequences (blue) or belonging to the same genus (gray) also captured in a field sample—either any field sample (hatched) or the matched site (solid). (C) Pairwise UniFrac distance (higher indicates communities are less phylogenetically similar, weighted by taxa abundance) between field and paired inoculum plotted against the percent permeable surface area in a 0.5 km radius. Line and shaded region indicate the prediction and 95% HPDI for a linear model of field-inoculum UniFrac distance and percent permeable surface area.

Cultured microbial communities differed substantially from field communities, but shared some similarities. With our analysis settings, ASVs that overlapped between a cultured community and the matching field community sometimes comprised a major proportion of the cultured community (Figure 1B, solid blue bars, average 15.8%), but ASVs identified only in a different field sample also often comprised a large portion (Figure 1B, hatched blue bars, average 52.3%). However, we caution against overinterpretation of these overlap results. Both field and cultured communities represent small samples (⁠$\approx$50 and 2 fronds, respectively) from nonoverlapping individuals, and the amount of strain-level turnover we would expect between them is unclear. Likewise, because sequences are separated at the exact variant level, this may split taxa below the species level. Indeed, the ASV overlap between field and cultured communities was highly sensitive to sequence analysis choices (compare Figure 1C; Supplementary Figures S2C and S3C), and when grouping sequences at a level of genera, the overlap of cultured microbiomes with taxa from matched field microbiomes increased (Figure 1B, solid gray bars, average 52.4%).

Overall community similarity based on weighted UniFrac distance (accounting for both relative abundance and phylogenetic similarity) was highly robust to analysis choices (Supplementary Figure S4) and should be less sensitive tostrain-level turnover across individuals. Matched culture-field pairs from urban sites were closer in UniFrac distance than matched pairs from rural sites (Figure 1C, $\text{pMCMC}<.05$⁠, see also Supplementary Figure S4, note the contrast with ASV or genera overlap in Figure 1B), potentially in part due to higher abundance of Enterobacteriaceae in cultured communities from urban sites, which was also common in one urban field sample (Figure 1A). Members of this family are often culturable and are known signatures of urban water (McLellan et al., 2015). A Mantel test of weighted UniFrac distance matrices between all field communities and weighted UniFrac distance between all cultured communities indicated a nonsignificant positive correspondence (⁠$\text{𝑝}>.1$⁠).

Variation in fitness and traits caused by treatments

Treatments explained 17–32% percent of variation across fitness and traits, explaining the greatest proportion of variation for duckweed fitness and the least for duckweed frond roundness (Figure 2A). Duckweed line (⁠$G_{host}$⁠) was an important source of variation in duckweed growth and trait values (explaining 21%, 18%, 10%, 12%, and 6.6% of the variation in duckweed frond area, total microbial density, frond color intensity, and frond roundness, respectively). Microbial treatment alone (⁠$G_{microbe}$⁠) was a substantial source of variation for color intensity (explaining 11% of variation) and aggregation (5.1%), but not for other measures. We also observed shifts in host contributions to trait values across microbiome contexts. $G_{host}\times G_{microbe}$ interactive effects between host line and microbial treatment consistently explained around 5% of variation across duckweed growth and trait values (4.6%–6.6%, Figure 2A).

Figure 2

(A) Percent of variation explained by treatments for duckweed fitness (frond area), each trait, and total microbial density (optical density at 600 nm, a measure of summed microbial growth across community members). See Methods for a description of how variation was partitioned; residual, unexplained variation makes up the difference from the top of the bar to 100%. (B) Broad-sense heritability of duckweed fitness and traits calculated in each unique zinc-by-microbiome treatment (top row) and the ratio of the variance among microbiome means to the total phenotypic variance (“heritability” in microbiomes) in each zinc-by-duckweed genotype treatment (bottom row). The dashed line in each panel indicates the detection threshold (upper bound of the 95% HPDI for the metric on the y-axis calculated from randomly permuted data).

We did not observe large contributions of zinc environment to trait variation. Zinc treatment alone (⁠$E$⁠) was never a substantial source of variation (explained $<$1% of variation across measures, Figure 2A). A low contribution of the main effect of $E$ is not unexpected, since we sampled only two environments, and many other environmental factors vary across duckweed habitats. For example, unmanipulated aspects of the experimental environment contributed to fitness and traits (spatial effects, Supplementary Tables S4–S6). The contributions of host-by-environment (⁠$G_{host}\times E$⁠) and microbe-by-environment (⁠$G_{microbe}\times E$⁠) were relatively small, explaining $<$1% each of duckweed frond area, microbial density, and aggregation variation. There was remaining variation explained by the combination of host-by-microbe-by-environment effects (⁠$G_{host}\times G_{microbe}\times E$⁠), from 3.2% to 4.1%, potentially indicating complex impacts of zinc environments on trait expression mediated by hosts and microbiomes.

Broad-sense heritability for duckweed fitness and traits (the variance in duckweed genotype means divided by total phenotypic variance) shifted across both zinc level and microbiome treatment (Figure 2B, top row). Duckweed fitness (frond area) and total microbial cell density were often estimated to be more heritable than the detection threshold (⁠$\approx$0.23 across traits), depending on conditions. For example, heritability in duckweeds of microbial density reached 0.62 when duckweed genotypes were growing in microbes from the University of Toronto Mississauga at high zinc (Figure 2B). Other traits less often exceeded our detection threshold.

We also report the “heritability” of duckweed fitness and traits in the microbiome: the variance in microbiome treatment means divided by total phenotypic variance, including nongenetic sources of variation (⁠$V_{G_{microbe}}$/$V_{total}$⁠, see Methods). The “heritability” of duckweed fitness and traits in the microbiome shifted across both zinc conditions and microbiome treatments (Figure 2B, bottom row). However, few estimates exceeded the detection threshold, and the estimates were generally smaller than broad sense heritabilities for traits in the host genome (the highest estimate was 0.48, color intensity on duckweeds from Moccasin Trail at low zinc). Which estimates exceeded our detection threshold depended on zinc treatment and duckweed host (Figure 2B).

The urban-to-rural ecotone shapes host, microbe, and interactive effects on traits and fitness

Across microbial communities, duckweeds from urban sites (low % permeable surface area) grew to much larger total frond area (pMCMC $<$ 0.05, Supplementary Figure S6, a $G_{host}$ effect). Likewise, microbial communities from urban sites reached higher cell density when inoculated onto duckweeds from rural sites (a $G_{host}\times G_{microbe}$ effect, pMCMC $<$ 0.05, Supplementary Figure S6). Duckweeds from rural sites inoculated with microbial communities from urban sites also had fronds with less color intensity and aggregation (Supplementary Figure S6, both pMCMC $<$ 0.05, but in high zinc only for aggregation, $G_{microbe}$⁠, $G_{host}\times G_{microbe}$⁠, and $G_{host}\times G_{microbe}\times E$ effects, see also Table S6). Across duckweed lines and microbial communities, zinc negatively affected duckweed fitness (from an average of 5.96 to 5.67 mm${}^2$⁠, $\pm$ a standard error of 0.10 for both), but had no impact on total microbial density (Figure 3A and D).

Figure 3

Duckweed fitness (frond area, A-C) and microbial density (optical density at 600 nm, a measure of summed microbial growth across community members, D-F). A & D: Effect of zinc. Means across all other treatments are shown in points and $\pm$ one standard error of the mean as bars. B & E: Means and standard errors in high (red) and low (black) zinc treatments for sympatric (“S”) or allopatric (“A”) combinations of duckweeds and microbial communities. C & F: Relationship between duckweed frond area (C) or microbial density (F) and the pairwise UniFrac distance between the field microbiome of the duckweed host and the experimentally inoculated microbiome. Predictions from the inverse scale are back-transformed to the data scale for optical density at 600 nm. Significance and model parameters in Supplementary Tables S4 and S5. Note change in y-axis scale across panels.

We did not identify any community composition differences that fully reflected observed functional differences in microbiome effects. The relative abundance of the most common family (Aeromonadaceae) did not explain variation in the effects of cultured microbial communities better than the percent permeable surface area surrounding its site of origin (DIC indicated worse fit). The relative abundance of Pseudomonadaceae, the second most common family, explained variation in microbiome effects on frond roundness only: higher abundance of Pseudomonadaceae in the inoculated community was associated with increased frond roundness for duckweeds from urban sites (pMCMC $<$ 0.05, SupplementaryFigure S7).

Shared history contributes to $G_{host}\times G_{microbe}$ effects on fitness

Microbiomes cultured from the same sites as duckweed hosts, or microbiomes that were more similar to “home” microbiomes, were not equivalent to other microbiomes in their effects on duckweed fitness. Sympatric combinations of duckweeds and microbiomes nonsignificantly increased duckweed frond area (our fitness proxy) at low zinc conditions, but marginally decreased frond area at high zinc conditions (by about 8%, at 90% HPDI, Figure 3B, see also Supplementary Table S4). Similarly, the effect of combining duckweeds with sympatric (“home”/“local”) microbial communities versus allopatric combinations of duckweed and microbes (all nonhome and nonlocal combinations) affected microbial cell density only in high-zinc treatments, where it negatively impacted total microbial density (⁠$\text{pMCMC}<.05$⁠, Figure 3E) with an average decrease in optical density from 0.072 to 0.067 (see also Supplementary Table S5). On the other hand, cultured microbiomes that were more similar to field communities isolated from the duckweed host significantly increased the frond area of duckweeds and microbial density (both pMCMC $<$ 0.05, robust to sequence analysis choices, Figure 3C and F; Supplementary Figure S5). There was about a 30% increase for duckweed frond area from least to most similar microbiomes, and a shift in optical density from 0.067 to 0.084. There was no significant interaction between community similarity and zinc treatment.

Microbiomes and environment altered the strength of selection on duckweed traits

We saw evidence for variable strength of selection on measured phenotypes across different zinc and microbiome contexts, but the direction of selection usually did not change. Depending on the inoculated microbiome and zinc, higher values of aggregation, color intensity, or roundness were significantly associated with increased fitness (Figure 4). Univariate selection gradients were most often significant for aggregation and were only negative for frond roundness (and in only one treatment). A post hoc analysis suggested that microbiome origin along the urban-rural gradient could explain variation in the strength of selection on microbial density across zinc and microbe treatments (Supplementary Figure S9).

Figure 4

Estimated univariate selection gradients (fitted slope values) between duckweed genotype fitness means (frond area) and genotype trait means (including microbial density). Each point represents a different microbiome-by-zinc treatment, with points sorted along the x-axis according to the degree of urbanization of the microbiome source site (left = most urban, right = most rural), and point color denoting high (red) and low (black) zinc treatments. Lines indicate 95% HPDI for each estimated univariate selection gradient. Asterisks indicate the significance of the relationship as determined by whether model fit was better (DIC was lower) in the real data than 95% of the permuted data sets. See Supplementary Figure S8 for the full relationships described by each point.

More dissimilar microbiomes drive $G_{host}\times G_{microbe}$ effects

The effect of duckweed line on fitness, trait, and microbial density values remained similar only when comparing duckweeds grown in similar microbial communities. When the correlation of duckweed line means in one microbial community to line means in another is high, this indicates that the relative effects of these microbial communities on duckweeds were similar across all host genotypes, for example, that $G_{host}\times G_{microbe}$ effects are weak for that pair of communities. Conversely, when the correlation of duckweed line means across two microbial communities is low, this suggests that the relative effects of these microbial communities on duckweeds were not similar across duckweed hosts, for example, that $G_{host}\times G_{microbe}$ effects are strong for that pair of communities. The correlation of duckweed line means between pairs of microbial communities (here, trait similarity) was significantly negatively related to the UniFrac distance between those microbial communities (Figure 5, for duckweed frond area, microbial density, aggregation, and color intensity pMCMC is $<$ 0.05, but not for roundness, robust to sequence analysis choices: Supplementary Figure S11). Trait similarity ranged from 0.82 to 0.95 for the most similar microbiomes, but from 0.36 to 0.73 for the least similar microbiomes (across fitness, microbial density, and traits, excluding roundness). This indicates that when two microbiomes are dissimilar, $G_{host}\times G_{microbe}$ effects are stronger. The biggest shift in trait similarity across microbiome similarity was for aggregation, which ranged from 0.90 (95% HPDI 0.78–1) for the most similar microbial communities to 0.36 (95% HPDI 0.25–0.47) for the least similar.

Figure 5

Relationships between pairwise similarity in traits or fitness for the duckweed genotypes across each pair of microbial sources (y-axis) and pairwise microbial community distance of those same sources (x-axis). Plots from left to right represent different fitness and trait measures: duckweed frond area (duckweed fitness), microbial density (optical density at 600 nm, a measure of total microbial growth), duckweed frond color intensity, duckweed area:perimeter ratio (aggregation), and duckweed frond roundness. Lines and shaded regions represent predicted mean and 95% HPDIs for significant relationships (all but frond roundness, where this relationship is n.s., see Supplementary Table S7). Each point’s pairwise trait similarity (y-axis value) is a correlation between duckweed genotype means across one pair of inoculated microbiomes. See also Supplementary Figure S10, where each point here is expanded to a full panel.

Discussion

A focus of recent host–microbiome research has been determining whether microbiome community composition is heritable in hosts (Grieneisen et al., 2021; Peiffer et al., 2013), and whether this heritable variation responds to selection (Henry et al., 2021; Mueller and Linksvayer, 2022; O’Brien et al., 2021). Yet, if selection on hosts directly impacts the relative fitness of different microbial species or strains, responses of microbiomes could contribute to trait change in hosts without any shift in host genotypes: via ecological shifts in microbial communities (Lau and Lennon, 2012), or evolutionary changes within individual microbes (Batstone et al., 2020; Rebolleda-Gómez et al., 2019; Tso et al., 2018). Whether responses to selection on host traits include microbial contributions is expected to hinge only on the amount of host phenotypic variation encoded by microbes that are transmitted to new hosts (Henry et al., 2021; Mueller and Linksvayer, 2022; O’Brien et al., 2021). Going forward, our ability to engineer microbiome variation for desired effects on hosts depends on our ability to predict microbiome effects and the amount of host trait variation in microbiomes. Here, we quantified variation in duckweed fitness and phenotypes that is attributable to the cultured fraction of their microbiome. We found that microbiomes contribute to traits, often via host- or environment-dependent mechanisms.

Host genotype underlies most variation, but effects are only predictable in similar microbiomes

Duckweed line was the largest source of phenotypic variation, though some effects depended on the inoculated microbiome (Figure 2). Other studies have found that most host trait variation is explained by host genotype, or host-specific impacts of microbes and microbial communities (Epstein et al., 2023; Heath et al., 2012; O’Brien et al., 2019; Wagner et al., 2014; Wendlandt et al., 2021). However, few studies have evaluated microbial effects on variation available to selection (see Chaney and Baucom, 2020; Epstein et al., 2023; Fitzpatrick et al., 2019; Heath et al., 2012; O’Brien et al., 2019; Wagner et al., 2014; Wendlandt et al., 2021). That most duckweed trait variation is explained by duckweed genotype is somewhat surprising, given the low genetic variation among duckweeds revealed by sequencing (Ho, 2017). However, host genotype here includes maternal effects such as epigenetic marks, an important source of phenotypic variation in clonal plants (Wilschut et al., 2016).

In contrast, microbiome contributions to host traits and fitness primarily depended on host genotype or environment (Figure 2). Across hosts and environments, identically assembled communities could produce different effects, or communities could assemble differently. The high-throughput, small-volume nature of our experiment prevented us from evaluating community assembly via sequencing. However, previous research indicates that host genotypes often affect the relative growth of taxonomic groups within inoculated communities (Grieneisen et al., 2021; Jackrel et al., 2021; Lebeis et al., 2015; Wippel et al., 2021), altering microbiome assembly via top-down host selection of different microbes (Lebeis et al., 2015), or bottom-up adaptation of microbes to better colonize local hosts (Li et al., 2021). Abiotic environments can also alter microbiome assembly (Wang et al., 2021), as tolerance to stressors including zinc varies across microbes (Davis et al., 2004; Jain et al., 2020), potentially accounting for zinc-dependent effects of microbiomes.

Since we did not verify the elimination of endophytic microbes before inoculations, our estimates of phenotypic variation across host genotypes could include variation contributed by vertically transmitted endophytes. However, this less manipulable fraction of the microbiome should perhaps be considered host-derived variation: tightly vertically transmitted microbes may not segregate separately from host genomes (Moran et al., 1993). Furthermore, because our microbial treatments affected host traits and fitness, some trait variation must reside in the microbes we did manipulate.

Phenotypes resulting from host and microbiome variation were somewhat predictable: more similar inoculated microbiomes resulted in conserved effects of host genotypes on traits, but less similar microbiomes resulted in unpredictable phenotypes across duckweed genotypes (Figure 5). This pattern suggests that for sufficiently similar communities, host-by-microbiome interactive effects (⁠$G_{host}\times G_{microbe}$⁠) may be minimal. When engineering small shifts to microbiome communities, we may be able to model microbiome effects on host traits as largely additive (e.g., Falconer, 1996), simplifying efforts to construct beneficial microbiomes. Reciprocally, this result also suggests that host genotypes that perform poorly in one microbial condition may perform well in a very different microbial condition.

Any benefits of shared history between duckweeds and microbiomes disrupted by zinc

In other systems, plants often perform better with microbial communities from their local site, suggesting local adaptation between plants and microbiomes is common (Rúa et al., 2016). Furthermore, experimental evolution has documented the adaptation of microbes to increase both plant and microbe fitness in local combinations (Batstone et al., 2020). Here, sympatric and allopatric combinations of hosts and microbiomes revealed only local maladaptation; in high-zinc conditions, combining duckweeds and microbial communities from the same site reduced duckweed fitness and microbial density (Figure 3B and D).

Benefits of locally adapted mutualisms can depend on environmental conditions matching the conditions under which local adaptation occurred (Bronstein, 1994; O’Brien et al., 2018). For example, local adaptation between a plant and its mycorrhizal fungus appeared greater when they grew in soils from the home site (Johnson et al., 2010). Zinc may be one such condition: different responses of duckweed lines and microbial communities to zinc (Supplementary Figure S12), could reflect adaptation to local zinc conditions, which are anticipated to vary (Ontario Ministry of the Environment, 2011), and have elicited adaptation in other plants (Babst-Kostecka et al., 2016; Schvartzman et al., 2018) and microbes (Davis et al., 2004). Urban environments shape evolution in many species (Johnson and Munshi-South, 2017) and may therefore be another such condition. Like our duckweeds, urban populations of various species grow faster or sooner than rural populations (Brans & De Meester, 2018; Gorton et al., 2018; Santangelo et al., 2020, see Supplementary Figure S6). Additionally, culturing-induced shifts in microbiomes may cause us to underestimate the benefits of local microbiomes if other microbiomes are more similar to the field microbiome than the “home” microbiome (Supplementary Figure S4, Figure 3, next section).

A closer “match” to the natural microbiome matters for both hosts and microbes

We saw that microbiomes more similar to the community present on duckweed lines in the field consistently supported greater duckweed fitness and increased microbial cell density (Figure 3C and F). By this metric, we would expect that local adaptation increases duckweed fitness and microbial growth in the field, in contrast with our results from the allopatric-sympatric test. If this is a general pattern, studies evaluating single microbes outside of their community context or using a culturing or storage step may underestimate the benefits of local microbes (e.g., in estimates of the frequency of local adaptation in plant–microbe interactions, Rúa et al., 2016). Indeed, some cultured and field communities deviate substantially (Figure 1). More widespread use of holistic distance metrics for experimental and home microbiome “environments” (developed for studying adaptation to abiotic environments, Montalvo and Ellstrand, 2001; Wang et al., 2010; Wilczek et al., 2014) may help quantify local adaptation to complex microbial communities that are difficult to perfectly replicate.

Inocula more similar to field microbiomes may be more beneficial because the culture step may select for bacteria that grow well outside the host environment, whereas the most beneficial bacteria likely grow better in hosts than in culture. Panke-Buisse et al. (2017) found that conditions in an intervening culture or storage phase drastically influenced microbiome effects on host traits, and Burghardt et al. (2018) found that microbial fitness in the host is only weakly correlated to microbial fitness in cultures. Aquatic bacteria vary in swimming and biofilm-forming ability, including abundant taxonomic groups in our field and cultured microbiomes (e.g., Rhodobacteraceae, Aeromonas, Pseudomonas, Bartling et al., 2018; Pernthaler, 2017; Rossi et al., 2018; Talagrand-Reboul et al., 2017), which may affect growth in lab cultures versus on duckweeds.

Other studies of cultured “synthetic” microbiomes may similarly miss testing the most beneficial microbiomes. Synthetic microbiome methods have taxonomic biases (Rappé and Giovannoni, 2003) or involve re-assembly from low abundance (Jones et al., 2022), which may introduce deterministic shifts in benefits. Many beneficial effects of microbiomes may be lost during sampling and culture, suggesting that microbiome scientists should sample inocula more similar to field microbiomes (e.g., Christian et al., 2017; Zahn & Amend, 2017). Yet, amplification, storage, and inoculation seem to be necessary for applying engineered microbiomes, and overcoming the shifts they introduce remains a major hurdle (Kaminsky et al., 2019), as confirmed by our results. Hopefully, novel methods reducing sampling and culturing bias (Thrash, 2021) will eventually scale to production level.

Microbiomes alter trait heritability and estimated strength of selection on duckweed traits

That microbes affect trait values expressed by their hosts is well known (Friesen et al., 2011; Gould et al., 2018; Honegger, 1993; Turnbaugh et al., 2008). Only a handful of studies have tested whether microbes alter estimated trait heritability in host genomes or selection on host traits, but all have found they do (Fitzpatrick et al., 2019; Horn et al., 2013; Lau and Lennon, 2011; O’Brien et al., 2018; Wagner et al., 2014, Figures 2 and 4). While most work to date considers plant hosts, shifting microbiome composition in fruit flies produced signals of selection in host genomes, implying that microbiomes change selection on or heritability of fly traits (Rudman et al., 2019).

Our selection analysis results suggested that duckweed traits were positively genetically linked to duckweed growth in only some microbial communities and zinc contexts. Both biotic and abiotic contexts can therefore shift the strength of selection. However, the field contains additional biotic and abiotic contexts, and selection may differ in direction and strength, or the true targets of selection could be other, genetically correlated traits.

Conclusions

While plant genotype harbors the most variation for plant traits, microbiomes do alter traits, fitnesses, and their heritability in host- and environment-dependent ways. Yet, releasing engineered microbiomes into nature could have unintended effects (Jack et al., 2021). We will need to understand the potential ecological and evolutionary consequences of engineered microbiomes to avoid negative outcomes. Our results here suggest that not all synthetic microbiomes will be equally beneficial: those that differ greatly from what a host naturally encounters may reduce host fitness. Furthermore, we found that the reduction of $G_{host}\times G_{microbiome}$ effects for similar microbiomes may increase the predictability of synthetic microbiome effects on hosts. Therefore, synthetic microbiota that better match natural microbial communities could reduce risks while increasing benefits and trait predictability.

Data availability

Data, scripts, and select intermediate files are accessible via GitHub at https://github.com/amob/duckweedgxgxe, an archived repository. Data are also archived on figshare at DOI 10.6084/m9.figshare.24059196.

Author contributions

All authors contributed substantially to the design of the study, provisioning of materials, or revising of the manuscript. A.M.O., M.E.F., J.R.L. proposed the study; M.E.F. and J.R.L. provisioned materials; A.M.O. collected the data and performed analyses; A.M.O. provided the first draft of the manuscript, and all authors revised.

Funding

The project was funded by the Gordon and Betty Moore Foundation through Grant GBMF9356 to M.E.F. (https://doi.org/10.37807/GBMF9356) and a Natural Sciences and Engineering Research Council of Canada Discovery Grant to M.E.F.

Conflict of interest

The authors declare no conflict of interest.

Acknowledgments

The authors thank E. Lash for help collecting duckweeds, students and volunteers who have contributed to maintaining duckweed and microbe cultures in the lab, and members of the Frederickson lab for useful discussion.

Footnotes

1

Accumulating recent evidence suggests that Greater Toronto Area Lemna may in fact be more appropriately referred to as Lemna$\times$japonica (Braglia et al., 2021a, b)