Microbes can have strong direct effects on the fitness of eukaryotic organisms. Similarly, eukaryotic organisms can strongly affect the composition (Wagner et al. 2016;Gehring et al. 2017;Walters et al. 2018) and evolution of microbial populations (Toft and Andersson 2010;Guidot et al. 2014;Laine et al. 2014). Plantsoil feedback (or "legacy") experiments show that the biotic and abiotic environment at one point in time can alter the growth of plants in future generations (Kardol et al. 2007;Lau and Lennon 2011;Shapira 2016), at least partly through changes in the soil microbial community (Haichar et al. 2006;Kardol et al. 2007;Meisner et al. 2013). These same factors could impose selection and drive adaptation within microbial populations, although the influence of environmental legacies on host-imposed selection has not been previously examined. In the present study, we use a legume-rhizobia system to determine how prior exposure to abiotic and biotic selection affects the relative fitness and diversity of bacterial symbiont associations with host plants.

Most plant-associated microbial species experience temporally and spatially variable selection, because they alternatively live in the soil and in close association with plants (Burghardt 2019;Provorov and Vorobyov 2006). For the microbes, the net selection experienced by the microbial population is determined by the frequency each environment (i.e., soil or host) is encountered, the relative strength of selection, and fitness correlations between environments (Levins 1968;Siepielski et al. 2009). Fluctuations between two or more different selective environments are expected to contribute to the maintenance of greater genetic diversity in a microbial population than what is expected if that population were exposed to only one of the selective environments.

Measuring selection acting within microbial populations in seminatural environments is challenging. Studies that have estimated fitness of symbiotic microbes within hosts have primarily relied on fluorescent tags or a small number of genetic markers and have thus been limited to examining the relative fitness of only a few strains (e.g., Svenning et al. 2001;Heath and Tiffin 2009) or specific mutants (Mellor et al. 1987;Daubech et al. 2017). A handful of studies have examined adaptation to hosts across multiple generations using experimental evolution approaches (Guidot et al. 2014;Marchetti et al. 2017;Meaden and Koskella 2017), but

these have focused on the effects of de novo mutations in a single genetic background. In nature, however, many microbial populations are genetically diverse (Shapiro et al. 2009;Karasov et al. 2018). High-coverage pooled sequencing together with bioinformatic estimation of strain frequencies provide a means to assess selection acting on multistrain populations (Kessner et al. 2013;Burghardt et al. 2018) and thereby provide an opportunity to directly measure the cascading effects of host and soil selection and legacies on microbial fitness and diversity.

The facultative mutualism between leguminous plants and rhizobial bacteria has been intensively studied, primarily because rhizobia, when growing in symbiosis with legumes, convert atmospheric nitrogen into a plant usable form (Dilworth et al. 2008;diCenzo et al. 2019;Poole et al. 2018). Here, we use a genetic and molecular model system of Ensifer bacteria and Medicago hosts. During the host growing season, a very small fraction of the soil Ensifer population will form a symbiosis with the host. A symbiosis occurs when a bacterium enters a compatible host root hair and initiates the formation of a nodule. In nodules with terminal bacteroid differentiation, such as those formed by Medicago, once inside the nodule, the bacterium divides and forms two subpopulations, a nonreproductive population that fixes nitrogen and a reproductive population that is released into the soil when the nodule senesces (reviewed in Sachs et al. 2018).

A single bacterium that forms a nodule can leave hundreds of thousands of offspring and thus plant hosts can impose strong selection on rhizobial populations (Burghardt et al. 2018) by either limiting which strains form nodules or by differentially rewarding or sanctioning strains once they are inside plants (Denison and Kiers 2011;Daubech et al. 2017;Sachs et al. 2018). Several lines of evidence indicate that evolution in the soil is also important: rhizobia maintain genes necessary for free-living growth on diverse substrates (Galibert et al. 2001), there can be high strain diversity and large populations in soil even when hosts have not been present for multiple years (Hirsch 1996;reviewed in Graham 2008), and genes that confer the ability to form a symbiosis with plants move from nonnative bacteria into native strains that are presumably better adapted to the local soil environment (Sullivan et al. 1995;Nandasena et al. 2006;Perrineau et al. 2011;Hollowell et al. 2016).

Here, we investigate the cumulative effects of multiple bouts of selection imposed by plant hosts and soil on a synthetic community of rhizobia. We ask the following questions: (1) How does the diversity, fitness, and benefits of the strain community in nodules change across host genotypes and growth cohorts? (2) To what extent is the community affected by selective legacies (host presence, host genotype)? (3) What are the putative genetic targets of selection? We accomplish these objectives by exposing a synthetic community of Ensifer meliloti strains to three bouts of selection involving combinations of plant hosts and soil (Fig. 1). By har-vesting hundreds of nodules from each replicate and sequencing the strains inside them, we measured whether prior selection on the rhizobial community has cascading effects on selection and community diversity inside host nodules. Our results highlight the importance of selection acting in the multiple environments experienced by microbial symbionts.

We subjected a synthetic community of 101 E. meliloti strains (see Burghardt et al. 2018 for details on strains) to three bouts of selection. In the first bout, the community was inoculated into 1-L pots of sterile low-nitrogen peat soil (SUN GRO LP5 Sunshine mix) that contained no plant (soil only, 16 pots) or pots with 10-12 surface-sterilized seeds of one of two Medicago genotypes-16 pots each of A17 Jemalong (HM101) or R108 (HM340). Six weeks after inoculation, we removed the vegetative portion of the plants and allowed roots and nodules to decompose in the soil for six weeks; this was the second bout of selection. Together, these 12 weeks were the "legacy" period. To identify the strains with high fitness in host nodules during the legacy plant cohort, we sampled nodule communities from a separate set of 5-weekold A17 (A17 legacy , eight pots) and R108 (R108 legacy , nine pots) plants (these data are described in Burghardt et al. 2018). We then imposed a third bout of selection by planting 12 surface-sterilized seeds of either host A17 or R108 in the soil-only pots (eight pots each, hereafter referred to as soil → R108 and soil → A17), 12 seeds of R108 into seven pots where A17 had grown previously (A17 → R108), and 12 seeds of R108 into the seven pots where R108 had grown previously (R108 → R108) (Fig. 1, due to space and sequencing limitations neither A17 → A17 nor R108 → A17 treatments were included). Five weeks after planting, we harvested nodules from the focal cohort and quantified the strain community (see below). At each stage of the experiment, all pots were randomized and left on a single greenhouse bench. Each week, pots were provided with 150 mL of nitrogen-free fertilizer (Bucciarelli et al. 2006) and sterile water as needed. Plants grown in several uninoculated pots that were randomized with other pots showed little/no evidence of nodulation, indicating that there was little movement of rhizobia between pots or contamination by rhizobia that were not part of the initial inoculum.

We measured strain communities inside host nodules using a select and resequence approach, a variant of evolve and resequence approaches (Burghardt et al. 2018;Burghardt 2019). In brief, we grew each of 101 E. meliloti strains, the genomes of each of which had been previously sequenced, in 3 mL Tryptone yeast media (TY) for three days and then combined an equal volume of each culture to generate a community with approximately equal representation of each strain. One milliliter of this community, approximately 10 9 cells, was used to inoculate each 1-L pot at the beginning of the experiment. When plants were five weeks old, we sampled ß250 nodules from the plants in each pot (Table S1), crushed them, used a series of differential centrifugation steps to isolate reproductive bacteria, and extracted DNA using the Ultr-aClean Microbial DNA Isolation Kit (no. 12224; Mo Bio Laboratories). This DNA, as well as DNA from the initial community, was used to generate NexteraXT libraries that were sequenced using an Illumina HiSeq 2500. We obtained 250 bp paired-end reads for all legacy and half of the initial inoculum libraries and 125 bp paired-end reads for all focal generation libraries and the other initial libraries. Reads were trimmed with TrimGalore (version 0.4.1; bioinformatics.babraham.ac.uk/projects/trim galore/) using default settings except that the quality parameter was set to 30, minimum length was set to 150 bp (100 bp for the shorter reads), and minimum overlap with an adapter sequence required to trigger adapter trimming was set to 3. Variant frequencies were called by aligning to the USDA1106 genome (Nelson et al. 2018) using bwa mem (version 0.7.15) (Li 2013) with default settings.

To eliminate reads that mapped to more than one location, the minimum mapping quality was set to 30. A median of 78.8% of reads from each library (range 57.4-89.9%) mapped to the USDA1106 genome resulting in a median of 5.7 (range 2.7-9.8) million aligned reads per library (Fig. S2; NCBI Accessions SRR6029825-SRR6029912). We estimated the frequency of each strain using HARP (Kessner et al. 2013). HARP calculates the probability that each read came from each strain and then finds the combination of strain frequencies that maximizes the likelihood across all reads. Strain frequencies inferred for initial communities were highly correlated across sequencing runs (the 250 and 125 bp reads, r = 0.99).

To characterize strain diversity in the nodule community, we calculated Renyi entropies (Rényi 1961) on strain frequencies estimated for each pot ("renyi" function in vegan package of R). Renyi entropy is calculated across a continuum of scale parameters (α) that vary in their sensitivity to strain frequencies: at low α, rare and common strains have equally important effects, whereas at higher α, common strains have greater weight than rare strains. Specific parameters produce commonly used diversity measures including richness (0), the exponent of Shannon diversity (1), inverse Simpson (2), and Berger-Parker (infinity).

To estimate the effects of treatments on strain diversity, we ran linear models ("lm" in base R) contrasting pairwise treatments of interest at each of these four commonly used scale parameters.

Results were robust to removal of extremely rare strains (frequency < 0.005) and fairly consistent across diversity parameters (full results in Table S2, Fig. S3).

We calculated the relative fitness of each strain in each replicate pot as the fold change in the frequency of that strain (q x ) in postselection communities relative to mean frequency of that strain across eight sequencing replicates of the initial community (log 2 (q x selected /qx initial )). The median q x initial = 0.0092 (5-95% quantile: 0.0054-0.0154). This transformation both normalizes the strain fitness distribution and controls for slight differences in strain frequency in the initial community. For instance, a strain with q x initial = 0.01 and q x selected = 0.08 would have a relative fitness of 3, whereas a strain with q x selected =0.00125 would have a relative fitness of -3. We used redundancy analysis (RDA, "rda" function in the vegan R package; Oskanen et al. 2017) to collapse the dimensionality of the relative fitness data and visualize changes in relative fitness in response to each treatment. RDA fits a multivariate linear regression to centered and scaled data and then uses principal component analysis (PCA) to decompose the major axes of variation in the fitted parameters. The adjusted r 2 of each RDA model provides an estimate of the proportion of variance in relative fitness explained by the model predictor(s) (i.e., treatment). We permuted the data to determine the probability that fitness differences occurred by chance (anova function, 999 permutations). The specific treatment comparisons made using RDA are presented in Table 1. Estimates of strain fitness variance across replicates and treatments using unconstrained, covariance-based PCA ("PCA" in "FactoMine" package in R) were qualitatively similar to those of the RDA (Fig. S4).

In addition to testing the effects of hosts and legacy environment on strain fitness, we conducted two sets of tests for differences in the strength of the correlations between relative fitnesses. First, we compared the correlation between the same host genotype in different cohorts to the correlation between different host genotypes in the same cohort. Second, we compared the strength of correlations between the R108 legacy cohort and either the focal cohort preceded by soil or the focal cohort preceded by a host.

To deal with the nonindependence of the data in these comparisons, we used Pearson and Filon's z test (Pearson and Filon 1898) and Zou's confidence intervals (Zou 2007) calculated using the "cocor" package in R (Diedenhofen and Musch 2015; other tests implemented in cocor produced similar results). We used a twotailed test of significance to compare correlations between the same host genotype in different cohorts and different host genotypes in the same cohort because we had no a priori expectation of directionality. We used a one-tailed test to compare the strength of the correlations between cohorts because we expected correlations to be stronger when the host was present during the legacy cohort. R code, strain frequency, and phenotype data are available on Dryad (https://doi.org/10.5061/dryad.qf36gd4).

We estimated how shifts in strain communities in nodules would affect plant hosts using single-strain inoculation data from Burghardt et al. (2018). In that experiment, we inoculated each host genotype with a single strain of rhizobia and measured plant dry biomass after 5 weeks of growth. For each replicate, we calculated predicted biomass as the sum of the per strain frequency in nodules multiplied by the dry plant weight from the single-strain inoculation experiment: n x=1 q x × dry plant weight in single strain experiment x n x=1 q x

Because the benefits plants obtain from rhizobia are host genotype specific, for A17 we used A17 single-strain data and for R108 we used R108 single-strain data. We made host-specific predictions comparable by placing them on a scale from zero (complete occupancy by the least beneficial strain) to one (complete occupancy by the most beneficial strain). We tested for differences in predicted biomass between treatments using a linear model. We have single strain data with 100 strains for A17 and 70 strains for R108, but we have no reason to suspect the subsample is biased. A caveat to these results is that the single-strain measurements were performed in an independent greenhouse experiment under different environmental conditions.

To identify the allelic variants showing the strongest statistical signal of change in frequency between treatments, we followed the procedure described in Burghardt et al. (2018). Briefly, each variant site was tested using a generalized linear model with a binomial link and read counts supporting each allele as the response.

The 300 variants (supplemental datasets on Dryad), out of 126,000 tested, with strongest statistical support were then grouped based

Table 1. Data and RDA model predictors used to address questions related to how strain fitness was affected by host genotype, cohort, and legacy type. Question Treatments Predictors in RDA model What is the relative importance of host genotype and growth cohort on strain fitness? A17 legacy R108 legacy soil → R108 soil → A17 Host genotype (A17, R108) Cohort (legacy, focal) What is the relative importance of legacy type versus host genotype on strain fitness? soil → R108 soil → A17 R108 → R108 A17 → R108 Focal host genotype (A17, R108) Legacy type (soil, host) What is the effect of host presence during the legacy period on strain fitness in focal hosts? soil → R108 R108 → R108 Legacy type (soil, R108) What is the effect of host genotype during the legacy period on strain fitness in focal hosts? R108 → R108 A17 → R108 Legacy type (A17, R108) on linkage disequilibrium (r 2 > 0.75). We tested four contrasts: R108 legacy versus soil → R108 and A17 legacy versus soil → A17 (combined effect of growth cohort and extended period of selection in soil), A17 → R108 versus R108 → R108 (effect of host identity during the legacy cohort), and soil → R108 versus R108 → R108 (effect of host presence during legacy cohort).

Hosts formed 189-325 nodules per pot with R108 hosts tending to make fewer, larger nodules than A17 (Table S1). Consistent with previous experiments (Heath and Tiffin 2009;Grillo et al. 2016;Burghardt et al. 2018), strain fitness depended strongly on host genotype and was skewed so that less than one-fifth of the strains increased in relative frequency inside nodules (relative fitness > 0; Fig. 2A). In both growth cohorts, the nodule community in A17, relative to R108, showed greater among-strain variance in relative fitness, fewer strains with relative fitness ࣙ0 (detailed metrics in Table S2), and lower strain diversity (Fig. 2B, legacy F df=1,15 = 25.9, P < 0.001; focal F df=1,14 = 21.4, P = 0.0013, details in Table S3). All of these are indicative of A17 imposing stronger selection than R108. Redundancy analyses (Fig. 3A) revealed that host genotype explained far more variation (37%) in rhizobial fitness than growth cohort (6.7%). Likewise, fitness was more strongly correlated across the two growth cohorts with the same host (A17: r n=101 = 0.95 and R108: r n=101 = 0.88; Fig. 3B and Fig. S1) than between hosts in the same growth cohort (legacy cohort r n=101 = 0.57, focal cohort r n=101 = 0.73; Fig. 3B). The differences in correlation strength were highly significant (P df=98 < 0.001) for each of the four pairwise comparisons (Table S4). Interestingly, strain fitness in the A17 and R108 nodule com-munities were more strongly correlated in the focal cohort (after three months of soil selection) than in the legacy cohort. Strains that had intermediate fitness in R108 legacy nodules and low fitness in A17 legacy nodules had low fitness in both hosts in the focal generation (Fig. S1). This pattern is likely due to selection in the soil environment, but we cannot exclude the potential effect of changes in host growth conditions across generations.

Strain communities inside each host were more beneficial than expected if hosts had formed associations in direct proportion to the representation of each strain in the initial community (Fig. 2C). In accordance with the stronger fitness alignment in A17 than R108 (Burghardt et al. 2018), A17 associated with a more beneficial community than R108 in both cohorts. Although there was no significant change in the predicted benefit of nodule communities across cohorts, R108 tended to associate with a more beneficial community of rhizobia in the focal than the legacy cohort. Taken together, strain fitness inside hosts is genotype specific with A17 consistently enriching a smaller and more beneficial subset of strains than R108.

If symbiosis strongly shapes the evolution of rhizobia populations, then we expect the legacy environment to affect strain frequencies in the focal cohort. After accounting for the 41% of strain fitness variation (P = 0.001) explained by focal host genetic differences (R108 vs. A17), legacy type (soil vs. host) explained 6.4% of the strain fitness variation (P = 0.022; Fig. 3C). Both the presence and identity of a host during the legacy period affected rhizobia relative fitness in the focal cohort. An RDA comparing R108 → R108 to soil → R108 treatments revealed that the presence or absence of a plant host during the legacy period explained 22% of the among-strain variation in fitness (P = 0.003, model r 2 adj. = 0.16, Fig. S5A) in the focal cohort, and comparison of A17 → R108 to R108 → R108 revealed that the genotype of the host during the legacy period accounted for 11.8% of fitness variation during the focal cohort (P = 0.047, model r 2 adj. = 0.043, Fig. S5B). Thus, the legacy environment affected the strain community in the nodules of the focal cohort plants, and the presence of a host during the legacy period was more important than the genotype of that host.

If high host-associated fitness in the legacy period enriches the soil for specific strains, then we expect strain fitness in the legacy nodules to be more strongly correlated with strain fitness in the focal cohort when the host was present during the legacy S4.

generation. Our results are consistent with this expectation (Fig. 3D and 3E); the correlation of relative fitness in the nodule community of R108 legacy hosts and R108 → R108 hosts (r n=101 = 0.93) is significantly stronger than the correlation between strain relative fitness in R108 legacy and soil → R108 nodules (r n=101 = 0.88; Pearson-Filon's z test P df=98 < 0.0001). Similarly, strain representation in A17 legacy nodules is more strongly correlated with strain representation in A17 → R108 nodules (r n=101 = 0.67) than with representation in soil → R108 nodules (r n=101 = 0.63), although the difference is not statistically significant (Pearson-Filon's z test P df=98 = 0.09, Table S4). The stronger correlations are not driven by extremely low fitness strains becoming less common. When strains with fitness less than -5 are assigned a value of -5, the R108 legacy result remains significant (P = 0.0003) and the A17 legacy correlation becomes significant (P = 0.02).

Given that plant hosts impose strong selection on rhizobia, we expected that two rounds of host selection would decrease strain diversity in nodules. Contrary to this expectation, the presence of a host in the legacy period resulted in a more, not a less diverse strain community in the focal cohort (Fig. 2B). When R108, the less selective host, was present during the legacy period, this effect was statistically significant across a range of community diversity metrics (Fig. S3; Table S3). When A17 was present during the legacy period, community diversity was still higher than it was with no host present during the legacy period, but the host effect was smaller. The effect was only marginally significant for richness and Shannon diversity metrics that emphasize the total number of strains present and not significant for metrics that are more sensitive to strain frequencies (Table S3). Interestingly, these shifts in diversity had little effect on the predicted host benefit of the strain communities (Fig. 2C).

The genomic variants most strongly affected by treatments are candidates for targets of environment-specific selection (Table S5). A comparison of allele frequency in the nodule communities of legacy versus focal cohorts (R108 legacy vs. soil → R108, A17 legacy vs. soil → A17) revealed that the effects of soil selection were largely independent of host identity; five of the top 10 variants with A17 were among the top 10 with R108 and all 10 were among the top 400 variants with R108. Allele frequencies of top variants changed dramatically across legacy and focal cohorts in both A17 (mean change 0.43) and R108 (0.40) hosts. Four of the variants most affected by soil selection were nonsynonymous variants that became fixed during soil selection. These variants were found in genes involved in resource use and include a ribose import protein (RbsA), a putative amino acid (D-alanine) exporter (Algl), a nitrate reductase involved in the first step of the denitrification pathway (NapC; Torres et al. 2014Torres et al. , 2018)), and a siderophore receptor involved in iron scavenging (FoxA; Cuív et al. 2008). Allele frequency differences between top candidates due to past host and soil legacies were, on average, weaker than direct shifts across cohorts (Table S5). The mean difference between top variants when R108 hosts were preceded by different host genotypes was 0.18 (A17 → R108 vs. R108 → R108) and by host presence versus soil was 0.19 (R108 → R108 vs. soil → R108). There was little overlap between the variants most strongly affected by legacy environments, suggesting that alleles enriched due to legacies of host genetic variation differ from those affected by a legacy of host presence.

Many microbial symbiont populations, including the vast majority of plant-associated symbionts, experience differences in selection as they alternate between living in soil and host environments (Cao and Goodrich-Blair 2017;Martiny and Walters 2018). Because rhizobia are horizontally transmitted, fluctuating selection can also occur through association with different host genotypes, assuming strain fitness depends on host genotype-as we found here and has been found in other legume-rhizobia systems (Pahua et al. 2018). Our results provide several lines of evidence, suggesting that both soil and host environments not only impose direct selection on rhizobia populations, but also have legacy effects across host generations. First, one host (A17) selected a less diverse and more helpful strain community than R108 regardless of whether the strain community experienced three months of soil selection. However, after soil selection, strain communities were far more similar in the two hosts with many of the same variants showing large changes in allele frequency. Second, we found that the presence of a host during the legacy period was more important than host identity during the legacy period. This result is consistent with experiments finding stronger rhizobial fitness correlations across different plant hosts than between plants and soil (Burghardt et al. 2018) and suggests that rhizobial strain fitness might be more affected by the abiotic environment than by genetic differences among hosts. Lastly, we found that strains with high host fitness in the legacy cohort tend to be enriched in the focal cohort, suggesting the potential for hosts to enrich microbial communities for specific strains. Despite this potential, two rounds of host selection resulted in increased, rather than decreased, strain diversity in nodule communities-perhaps because these communities were exposed to three selective environments (legacy host → soil → focal host) rather than two selective environments (soil → focal host). Further, strain communities after two rounds of R108 selection were not more beneficial to hosts, thus demonstrating that strain community shifts do not always feedback to influence host phenotypes. Of course, it is possible that if we had used the more selective host (A17) instead of R108, we would have observed a different outcome.

At first glance, our results showing that selection by plant hosts only weakly enriches strains with high fitness and actually promotes diversity in rhizobial communities appear to be in contrast to those of Daubech et al. (2017) who showed that plant hosts can cause rapid increases in the frequency of rhizobia that can provide plant benefit via nitrogen fixation (Fix +) when competed against strains that cannot (Fix -). However, Daubech et al.'s (2017) study differs from our study in two important ways. First, Daubech et al. asked whether a nitrogen-fixing strain could outcompete a mutant strain that differed only in the ability to benefit the host. By contrast, we investigated whether host-imposed selection could rapidly enrich a community of naturally occurring strains for those with the highest host-associated relative fitness, a metric that may not be strongly associated with plant benefit (Burghardt et al. 2018). Second, Daubech et al. allowed only the bacteria that formed nodules with a single host genotype to form the next generation of microbes, in effect, removing spatial environmental variability from the experiment. By contrast, we specifically incorporated the potential for selection on resident soil (and rhizosphere) populations as well as host genotype to shape community diversity (Marschner et al. 2002).

In general, we expected that repeated exposure to the same environment would reduce strain diversity and exposure to multiple environments that select for different strains to contribute to the maintenance of strain diversity. By contrast, we found the highest strain diversity in nodules when the same less selective host genotype (R108) was used in two generations of selection. Although we cannot pinpoint the mechanism responsible for this higher diversity, it might be due to soil selection that occurred between the legacy and focal plant cohorts, or due to selection by plant hosts that occur outside of the nodules, such as through root exudates (Haichar et al. 2006;Dennis et al. 2010) or niche construction (Dinnage et al. 2018). Regardless of the mechanism, our results suggest that bouts of selection outside hosts might be an underappreciated contributor to the maintenance of strain variation in rhizobial communities.

Theory predicts that temporal variation alone will not maintain genetic variation in populations at equilibrium (Felsenstein 1976;Mustonen and Lässig 2009). However, our experimental work is not examining an equilibrium situation, and the rhizobia experience spatially as well as temporally variable selection. Indeed, it is possible that few naturally occurring microbial symbiont relationships are at a selective equilibrium. If this is so, then to predict the potential importance of fluctuating selection on strain diversity we need to examine models in nonequilibrium states. Existing models suggest that the influence of fluctuations on genetic diversity depends on the length and predictability of selective environments as well as population sizes, the shape of fitness tradeoffs, and genetic architectures responsible for fitness variation (Cvijović et al. 2015;Saether and Engen 2015;Dean et al. 2017;Novak and Barton 2017;Wittmann et al. 2017). Unfortunately, little is known about these parameters in natural microbial populations, and fluctuating environments have, until recently, been neglected even in laboratory experimental evolution studies (reviewed in Collins 2011; but see Spencer et al. 2007;Turkarslan et al. 2017). Given that hosts can impose strong selection (Guidot et al. 2014;Meaden and Koskella 2017), studies examining how the frequency of exposure to hosts (e.g., Marchetti et al. 2017) as well as explicit consideration of host and microbial life histories (Barrett et al. 2008) are likely to be important moving forward.

In sum, our experiment highlights the importance of considering the multiple selective environments that plant-associated microbes experience when trying to understand adaptation of microbial populations. Clearly, this study is just a first step. Our understanding of how fluctuating selective environments shape rhizobial, and other bacterial, populations will be advanced by incorporating greater ecological realism and microbiome complexity in experiments (e.g., Wood et al. 2018); by designing experiments to tease apart the separate effects of soil, rhizosphere, and nodule selection on rhizobial populations; by conducting experiments for longer than two host generations; and by expanding the genomic diversity of host genotypes. Such experiments would not only advance our understanding of how selection affects microbial populations, but also may inform efforts to leverage beneficial plant-microbe interactions to increase agricultural yields and aid in ecosystem restoration (Gopalakrishnan et al. 2014;Toju et al. 2018).

LTB and PT designed and conducted the experiment, BE and LTB analyzed the data, LTB took the lead in writing the manuscript, and PT and BE provided extensive input on content and writing throughout.