To our knowledge, no other study has tested whether naturally occurring soil biotas from multiple undisturbed habitats affect flowering phenology and selection on flowering time, nor has any study of the relationship between soil microbes and phenology explicitly accounted for genetic variation among undomesticated plant populations.Sterilized potting soil was saturated with inoculum and sterilized field soils were saturated with buffer for 10 days. Then, 48 accessions of Boechera stricta were transplanted as gnotobiotic seedlings into all treatments . All pots were randomized into blocks of 200 and maintained under controlled greenhouse conditions throughout the experiment, except for 7 weeks of 4°C vernalization in a growth chamber. We measured plant height and leaf number on the date of first flowering. Here we focus on phenology and fecundity of the 51% of plants that flowered successfully. Among those plants that flowered, the experimental design exhibited only modest imbalance, with substantial sample sizes remaining in every subspecies × treatment cell . Factors influencing probability of flowering and other fitness components are beyond the scope of this paper; we found no evidence that any experimental treatments affected flowering probability of surviving plants . Reproductive fitness was estimated as the number of fruits on each individual at 33 weeks of age; in B. stricta, fruit set is strongly and positively correlated with seed production in the greenhouse .Soils were collected from four natural Boechera stricta habitats in central Idaho, USA, separated by ~26 to ~92 km and differing in elevation, temperature, water availability, density and diversity of vegetation, and many soil properties .

Collection locations were named ‘Jackass Meadow’ , ‘Mahogany Valley’ , ‘Parker Meadow’ , and ‘Silver Creek’ . These remote sites have little history of disturbance by humans, are home to endogenous B. stricta populations, hydroponic net pots and function well as common gardens for B. stricta field experiments . Therefore, they are legitimate potential habitats that B. stricta likely encountered during its evolutionary history in this region. Each soil collection comprised five sub-samples from the four corners and approximate center of a ~150 m2 area, at a depth of ~10 to 30 cm. Subsamples were combined and mixed thoroughly, sieved through ~1.25 cm wire mesh to remove rocks and coarse detritus, shipped to Duke University, and stored in plastic bags at 4°C for ~3 months until use. We also collected seven ~1 mL soil vouchers from each site for microbial community analysis: three in August 2011, and at all sites except PAR, four in August 2012. Vouchers were frozen at -20°C until DNA extraction in late 2012. B. stricta seeds were collected from 48 natural populations, including four from the soil collection sites. Their sites of origin span over 1,000 m in elevation and are separated by between ~1 km and ~350 km, with the exception of the “SAD12” genotype from Colorado. Because B. stricta is naturally inbred and exhibits high FIS and FST , each population was represented by a distinct genotype. This diverse collection of genotypes included 24 from each of the ecologically divergent “east” and “west” subspecies . We used seeds from a single mother descended from the original field-collected accession, self-pollinated in the greenhouse for at least one generation, to minimize maternal effects; i.e. individuals within a genotype were self-full sibs .To create four soils that were identical except for their microbial communities, we extracted microbes from field soils into sterile buffer and soaked sterilized potting soil in the resulting suspensions.

We prepared inocula from 75 g sub-samples of each field-collected soil stirred into 1 L of 2.5 mM MES monohydrate in sterile diH2O . After settling for 30 min the suspensions were vacuum filtered to remove particulates. Filtrates were centrifuged 30 min at 3,000×G at room temperature to pellet microorganisms. To remove dissolved nutrients, we discarded the supernatants and re-suspended the microbe-enriched pellets in 1 L sterile 2.5 mM MES. This process mostly eliminated variation in chemical properties that differentiate the field soils . Each rack of 200 pots was bottom-saturated with 400 mL of one of the microbial suspensions, 6 g 20-10-20 fertilizer, and sterile diH2O for a total treatment volume of 4 L. The fertilizer was added to encourage seedling survival and to counteract possible soil impoverishment due to autoclaving . An additional 1 mL of undiluted microbial suspension was pipetted into each pot. Treatments derived by this process are termed “biotic” or “microbial treatments” throughout. It is possible that the filtration and re-colonization processes somewhat altered community structures; however, the differences between our experimental inocula—and their effects on the plants—originate from corresponding differences between real Boechera habitats.We sterilized soils from four natural habitats to create growth substrates with different physical and chemical properties, but without their natural microbiomes. After sub-sampling to extract microbial communities , we sterilized the four field-collected soils via autoclaving . These soils were loosely packed into clean pots and bottomsaturated with 400 mL sterile 2.5 mM MES, 6 g 20-10-20 fertilizer, and sterile diH2O to bring the treatment volume to 4 L.

An additional 1 mL sterile 2.5 mM MES was pipetted into each pot. Treatments derived by this process are termed “abiotic” or “sterilized field soils”. Although it is likely that autoclaving these soils changed their fertility, they appear to have retained at least some of their natural chemical variation .Surface-sterilized seeds of 48 genotypes were stratified on auto claved filter paper at 4°C in the dark for two weeks, then placed in a growth chamber to germinate for one week . Four germinated seedlings per genotype were transplanted into each of the eight experimental soils described above, one seedling per pot. Eight pots per treatment were left unplanted as controls. All pots were immediately rearranged into randomized blocks and maintained in controlled greenhouse conditions for the duration of the experiment. Plants were top-watered as needed with RO water, and received an additional 4 mL 20-10-20 fertilizer via pipet when one month old. Two-month-old plants were transferred to a 4°C vernalization treatment, where they remained for seven weeks. After vernalization, plants were returned to the greenhouse, checked three times weekly for flowers and allowed to set fruit. Flowering was defined as sufficient separation of the corollasuch that four distinct petals could be identified. The number of days between end of vernalization and first flowering is termed interchangeably “flowering time” and “flowering phenology” throughout. On the day of first flower for each plant we measured the individual’s height and number of leaves. The last census was done eight weeks post-vernalization; the 749 plants that had not flowered by this date were excluded from all future analyses. The experiment ended two months after the final flowering census, when almost all fruits had matured and dehisced. These cutoffs for flowering time and fruit production are realistic given the short growing season observed in the field. At this time we counted the number of fruits produced by each individual.Due to current methodological limitations, in this study we focus on the prokaryotic component of the soil microbiome . We extracted DNA from field collected soil vouchers using the MoBioTM PowerSoil DNA Isolation Kit and amplified variable region 4 of the bacterial 16S rRNA gene using established primer pairs 515F and 806R and PNA PCR clamps to reduce plastid and mitochondrial contamination . Paired-end 2×250bp sequencing of barcoded amplicons was performed on a MiSeq machine running v2 chemistry at the Joint Genome Institute . The primer sequences were trimmed from the paired-end sequences, blueberry grow pot which were then overlapped and merged using FLASH . Merged sequences were grouped into operational taxonomic units based on 97% sequence identity, and chimeric sequences were removed, using the UPARSE pipeline . Taxonomies were assigned as in Lundberg et al. . Unclassifiable OTUs at the kingdom level, OTUs matching Viridiplantae, mitochondrial, or plastid sequences were excluded by using BLAST to compare them to a custom database of contaminant sequences . Unclassifiable OTUs at the kingdom level and rare or non-reproducible OTUs were also excluded as in Lundberg et al. , resulting in 7,844 OTUs. To control for unequal sequencing effort, we normalized data by rarefaction to 40,000 reads/sample using QIIME-1.7.0 . Diversity analyses were performed after correcting data for 16S gene copy number variation using scripts provided in Kembel et al. ; OTUs without taxonomic information were assigned the mean copy number . Linear regressions were performed before this correction but the resulting parameter estimates were adjusted as necessary.To test the hypothesis that soil properties affect flowering time, we used restricted maximum likelihood linear mixed models with treatment, subspecies, and treatment × subspecies as fixed factors; block, genotype nested within subspecies, and treatment × genotype as random factors; and elongation rate , height at first flowering , and leaves per mm stem as covariates .

To test for treatment effects on overall plant size we used MANOVA with treatment, subspecies, and treatment × subspecies as fixed factors. We analyzed the parallel “biotic” and “abiotic” experiments separately, i.e. one model tested only for effects of soil microbiomes on flowering time, and another identical model tested only for effects of physical soil differences; we did not directly compare the effects of these two types of soil variation. These models were run using JMP® Pro version 10.0.0 . Statistical significance of random effects was determined by REML likelihood ratio test and results were graphed using ggplot2 in R version 3.0.2 . To test the hypothesis that soil properties alter selection on flowering time, we used a REML linear mixed model with the same terms as above, plus flowering time and flowering time × treatment as additional fixed effects. The response variable was number of fruits. Thus, the flowering time term describes the change in fecundity attributed to a change in flowering time, i.e. selection. We analyzed the parallel “biotic” and “abiotic” experiments separately as above. We performed these models both with and without including elongation rate, height at first flowering, and leaves per mm stem as covariates; the former model describes the selection gradient on flowering time , and the latter model describes the selection differential . Introducing a quadratic term did not improve fit, so our model considers only linear effects . Sample sizes were slightly smaller than for the flowering time models because 17 individuals were accidentally discarded after flowering. Main models were performed in JMP; selection differentials and selection gradients were calculated in R. For some models, non-uniformly distributed residuals might have influenced judgments about significance. In general, results with and without covariates were similar, suggesting that significance was robust to heteroscedasticity in the standard model. For the sake of caution, for all of our major results, we performed permutation tests to verify the results of the standard ANOVA. Rarefied microbial communities were analyzed in the R package ‘vegan’ . Principal coordinates of Bray-Curtis pairwise dissimilarities were identified using the vegan function ‘capscale’. Similarity of samples within vs. among sites was tested using the non-parametric permutation test ADONIS with 9,999 permutations constrained by collection year. To ask which components of microbial communities affect flowering phenology and selection, we regressed mean flowering time in each biotic Treatment onto the mean PCo score from the corresponding Site. To ask which OTUs underlie the observed phenotypic effect, we identified the ten OTUs most highly correlated with the PCo axes and regressed the same flowering time residuals on the OTUs’ mean abundances at each site. We used the Wilcoxon rank-sum test to compare relative abundances of common taxa between groups of samples associated with extreme phenotypes. P-values were adjusted using Benjamini-Hochberg false discovery rate. Our method of searching for microbial community members that underlie our phenotype of interest is described in more detail in Appendix S3.Selection on flowering time depended on soil microbiome. This result held both for selection gradients and selection differentials . The most extreme change was between the PAR and SIL soil biotas , with selection differentials of +0.034 and -0.043 fruits/day , respectively. The magnitude of this difference in selection is 1.2 times the selection differential measured in a nearby field site .