Time outweighs the effect of host developmental stage on microbial community composition

ABSTRACT

Thousands of microbial taxa in the soil form symbioses with host plants, and due to their contribution to plant performance, these microbes are often considered an extension of the host genome. Given microbial effects on host performance, it is important to understand factors that govern microbial community assembly. Host developmental stage could affect rhizosphere microbial diversity while, alternatively, microbial assemblages could change simply as a consequence of time and the opportunity for microbial succession. Previous studies suggest that rhizosphere microbial assemblages shift across plant developmental stages, but time since germination is confounded with developmental stage. We asked how elapsed time and potential microbial succession relative to host development affected microbial diversity in the rhizosphere using monogenic flowering-time mutants of Arabidopsis thaliana. Under our experimental design, different developmental stages were present among host genotypes after the same amount of time following germination, e.g. at 76 days following germination some host genotypes were flowering while others were fruiting or senescing. We found that elapsed time was a strong predictor of microbial diversity whereas there were few differences among developmental stages. Our results support the idea that time and, likely, microbial succession more strongly affect microbial community assembly than host developmental stage.

rhizosphere, microbes, succession, developmental stage, Arabidopsis thaliana

INTRODUCTION

Soil-dwelling microorganisms are critical for ecosystem function (Berendsen, Pieterse and Bakker 2012; Delgado-Baquerizo et al. 2017), increasing stability of natural ecosystems (Delgado-Baquerizo et al. 2017), and enhancing sustainable agricultural management (Barrios 2007; Bhardwaj et al. 2014; Yadav et al. 2020). Thousands of diverse soil microbial taxa inhabit the region of soil that interacts with plant roots, known as the rhizosphere (e.g. Teixeira et al. 2010; Kamutando et al. 2017; Alawiye and Babalola 2019). These rhizosphere microorganisms can benefit their hosts in many ways; they may influence plant metabolic function (Hubbard et al. 2018), phenological timing (Berendsen, Pieterse and Bakker 2012; Lu et al. 2018), and plant productivity (van der Heijden, Bardgett and van Straalen 2008). Rhizosphere microbe communities can vary substantially over space and time (e.g. Dunfield and Germida 2003; Shi et al. 2015; Kamutando et al. 2017), exhibiting global patterns in composition across coarse-grained climatic gradients (Brockett, Prescott and Grayston 2012; de Vries et al. 2012) and, at finer scales, varying across meters to kilometers (Peiffer et al. 2013; Kamutando et al. 2017; but see Teixeira et al. 2010). Moreover, the diversity, abundance, and activity of rhizobacterial taxa can change seasonally (Shi et al. 2015; Qiao et al. 2017; Li et al. 2020) as well as diurnally (Hubbard et al. 2017; Staley et al. 2017; Baraniya et al. 2018). At the scale of a single host plant, as physiological requirements and environmental sensitivities change over its lifetime (Kalra and Lal 2018), the composition of rhizosphere microbial communities is also likely to differ (e.g. Micallef et al.2009; Chaparro, Badri and Vivanco 2014); the causes of such microbial changes over the host lifespan remain largely unresolved.

A range of host factors influence rhizosphere microbial community diversity (Berg and Smalla 2009; de Vries et al. 2012; Lareen, Burton and Schäfer 2016), including species of host plant (Berg and Smalla 2009), plant genotype (Micallef et al.2009; Micallef, Shiaris and Colón-Carmona 2009; Bulgarelli et al. 2012; Bodenhausen et al. 2014; Edwards et al. 2015), richness of plant community (Essarioui et al. 2017), plant traits such as growth rate and root length (Cantarel et al. 2015), and interactions among microbes that span from facultative to antagonistic (Barea et al. 2005; Philippot et al. 2013; Essarioui et al. 2017). There is also evidence that host plant developmental stage and age contribute to rhizosphere diversity (e.g. Micallef et al. 2009; Chaparro, Badri and Vivanco 2014; but see Dombrowski et al. 2017). Bacterial diversity in the rhizosphere changed with host age in Arabidopsis thaliana, though differently among host genotypes; by the end of the plants’ lives rhizosphere composition converged to similar communities (Micallef et al. 2009). Other examples include a study of peas (Pisum satvium), wheat (Triticum aestivum), and sugar beet (Beta vulgaris) that revealed changes in bacterial and fungal diversity in the rhizosphere microbiome significantly associated with developmental stage of peas and wheat, but not sugar beets (Houlden et al. 2008); in a study of potatoes (Solanum tuberosum), plant developmental stage had a marked effect on the composition of rhizosphere microbes (Rasche et al. 2006). While this growing line of evidence indicates that developmental stage is important for determining rhizosphere diversity, developmental stage is usually confounded with time in experiments (but see Dombrowski et al. 2017) and, therefore, it remains largely unclear whether progression through the host plant's developmental trajectory vs. time since germination and the opportunity for microbial community succession is responsible for change in rhizosphere microbes.

Plant developmental stage or age could affect microbial community composition via, for instance, changes in root exudation patterns that affect microbial nutrient supply or changes in root architecture that affect the physical habitat available for microbes (Saleem et al. 2018; Pervaiz et al. 2020). Alternatively, microbial succession, the sequential assembly of microbial communities in an environment, involves change over time and could account for shifts in community composition independent of host plant stage or age. In nature, disentangling host developmental stage from time and the opportunity for succession is impractical and likely impossible, but this separation is feasible in a manipulative study (e.g. Dombrowski et al. 2017). We designed an experiment to decouple time from developmental stage across multiple genotypes throughout the lifetime of those plants. Our purpose was to identify the contributions of plant developmental stage compared to time since planting on rhizosphere diversity. We performed this experiment in Arabidopsis thaliana (L.) (Brassicaceae), a model plant for which the genetic regulation of developmental stages (e.g. flowering time pathway) has been well studied (e.g. Weinig et al. 2002; Stinchcombe et al. 2004; Salomé et al. 2011) and can be experimentally manipulated through the use of existing knockout lines. In particular, the use of flowering-time mutants with single-gene knock-outs, leading to early- or delayed-flowering after a fixed amount of time, enabled us to characterize the microbial communities of plants at different developmental stages but with the same amount of time since germination.

Here, we present and test two hypotheses regarding changes in microbial community composition associated with host plants of Arabidopsis thaliana (Fig. 1). First, we hypothesized that the simple progression of time and opportunity for succession among microbes determines rhizosphere microbial composition. In this case we would expect to see distinct microbial assemblages at each of the sampling times regardless of the developmental stage of plants within each of those time points. Alternatively, if plant developmental stage determines rhizosphere microbial assemblage, then we would expect to see similarity of the rhizosphere communities from plants in the same developmental stage regardless of time since germination.

Figure 1.

Experimental hypotheses: (A) If time since germination determines microbial assemblage in the rhizosphere, microbial diversity should be distinct at each time point (grouped within dashed-line ovals) throughout the experiment regardless of developmental stage. For example, at T1 rhizosphere microbes (indicated by purple circles) should be the same for all individuals; at T2 microbes (indicated by purple circles and yellow squares) should be the same across developmental stages and distinct from T1, and so on for T3 (microbes indicated by blue triangles, red diamonds, and yellow squares) and T4 (microbes indicated by red diamonds, blue triangles, and pink pentagons). (B) Alternatively, if developmental stage of host plant determines rhizosphere microbial assemblage, plants should have different rhizosphere microbe diversities across developmental stages (grouped within dashed triangle and ovals) regardless of sampling time. In this case, we would expect to observe rhizosphere microbes that are the same among vegetative plants (microbes indicated by purple circles) across all relevant time points, flowering plants with the same microbes (purple circles and yellow squares) across time points, fruiting plants with the same rhizosphere microbes (blue triangle, red diamonds, and yellow squares), and senescing plants distinct from other developmental stages (microbes indicated by red diamonds, blue triangle, and pink pentagon).

MATERIALS AND METHODS

Experimental design

We tested our hypotheses experimentally with five genotypes of A. thaliana in the Col background: the wild type, two early-flowering knock-out mutants (svp-32 and tflp2-2), and two late-flowering knockout mutants (ap1-10, ft-10). The use of these genotypes allowed us to sample different developmental stages (early-, wildtype, and late-flowering) on the same day at progressive time points throughout the experiment (e.g. at time point 3, the ap-1 and ft-10 mutant genotypes would still be vegetative and not flowering, wildtype genotypes would be flowering, and the svp-32 and tflp2-2 genotypes would be fruiting, Fig. 1). Because the mutants differ at a single locus from each other and from the wildtype, genotypic effects on microbial composition should be limited, allowing for a direct test of plant developmental stage vs. time on rhizosphere microbe community composition. Further, we used multiple early- and late-flowering mutants to test and control for potential pleiotropic effects of the single-gene flowering-time mutations on microbial community composition.

To standardize growing conditions across treatments, we used a soil mix consisting of 95% autoclaved potting soil and a 5% inoculum of soil (v/v) collected in southeastern Wyoming (Rd234; 41.325912, -106.466076). Before planting we surface-sterilized seeds by washing them in 70% EtOH for one minute, 10% bleach for 10 min, and then rinsing with autoclaved RO water four times (Lundberg et al. 2012). We stored the clean seeds in sterile RO in the dark at 4°C for four days.

We grew plants in four growth chambers under 10 h: 14 h light-dark cycles with associated temperatures of 21.1°C and 15.5°C, respectively. We arranged plants using a randomized complete block design. Before germination, we watered pots from the bottom. After germination we removed pots from shared water to prevent exchange of microbes and misted them every day until they were large enough to be watered from above as needed.

Sample collection and processing

We collected rhizosphere samples at four different times: first when all plants were vegetative and then again when each of the genotypes flowered, which we identified as ≥ 50% of each genotyping bolting. This sampling approach allowed us to collect plants at fruiting and senescing stages as well. The first time point for collection (T1) occurred a month (30 days) after planting (30 June 2017). The second collection (T2) was at 51 days (21 July 2017) when the early-flowering mutants were in bloom. The third collection (T3) occurred on day 76 (15 August 2017), and the final collection (T4) took place on day 84 (23 August 2017). We collected a total of 21 rhizospheres at T1, with between 3–6 replicates per genotype. At time points T2–T4, the total number of samples that we collected per genotype was between 6−10 (Table S1, Supporting Information).

To collect rhizosphere samples, we removed plants from their pots, shook off the bulk soil, carefully removed remaining soil clumps greater than two millimeters in diameter, and placed the roots in sterile 50 mL centrifuge tubes. To each tube we added 40 mL of 1x PBS (with 0.02% Silwet L-77) and then vortexed the tubes for 15 min on the highest setting. To reduce bias introduced by root endophytic microbes, we separated the roots from the rhizosphere by passing each sample through a Steriflip filter of 100 μm nylon mesh (Millipore Corporation, Billerica, MA, USA) and retaining the filtered solution. We centrifuged tubes at 3200 RCF for 15 min. Following this spin, which created a pellet of the rhizosphere soil, we discarded the supernatant, added 1 mL of the PBS buffer, vortexed for five minutes, and transferred to a clean 2 mL microcentrifuge tube. We again centrifuged the 2 mL tubes at 3200 RCF for 15 minutes, poured out the supernatant, and flash-froze the tubes with liquid nitrogen. We stored these samples at −80°C until we extracted DNA.

We performed DNA extractions using the DNeasy Powerlyzer Powersoil kit according to manufacturer instructions (Qiagen 2017). For amplification and sequencing of the V4–V5 region of the 16S rRNA gene common to all microbes, we sent samples to the Marine Biological Laboratory (MBL; Woods Hole, MA, USA) for amplicon library preparation and paired-end sequencing (2 × 250 bp) on Illumina MiSeq (Illumina, San Diego, California, USA).

Data analysis

Sequence data processing

We used Trimmomatic version 0.36 (Bolger, Lohse and Usadel 2014) to remove adapters and primer sequences. We processed amplicon data in R Eggshell igloo (v3.52) (R Core Team 2020) using the DADA2 pipeline version 1.14.1 (Callahan et al. 2016). To filter reads we set the maximum number of expected errors to two (MaxEE = 2,2) and the minimum quality score to 11 (truncQ = 11), and removed all ambigious base calls (maxN = 0). Errors were learned using 1 × 10E8 reads. We then visually checked the error rates before merging forward and reverse reads, with a minimum overlap of 12 base pairs and 0 mismatches. The output of this process was tabulated, producing a sequence table. We then removed chimeras using the ‘consensus’ method and assigned taxonomy to the species level, where possible, using SILVA (SILVA 132, Quast et al. 2012). The final taxon table was then imported into phyloseq package version 1.26 (McMurdie and Holmes 2013). All reads that did not assign to the kingdom Bacteria reads removed prior to statistical analyses.

Statistical analyses

To test our hypotheses and identify relationships among microbes within and across sampling times, we looked for differences in both alpha (α) and beta (β) diversities. We performed comparisons of α- and β-diversity metrics (Richness, Shannon (H`), and Simpson (D)) on rarefied sequence data using Kruskal–Wallis tests followed by Dunn's tests in the stats and FSA packages respectively (Ogle, Wheeler and Dinno 2020; R Core Team 2020). Samples were rarefied to an even depth of 77213 reads. Multiple pairwise comparisons were accounted for using a Bonferroni correction. To examine α-diversity metrics visually, we used ggplot2 (Wickham 2016). Next, we quantified differences in β-diversity among samples using Bray–Curtis pairwise dissimilarities and plotted differences via non-metric multi-dimensional scaling (NMDS). Differences in β-diversity were assessed using within sample proportional abundances, calculated by dividing the counts of a single taxon by the total counts within a sample. To identify significant differences among genotypes, developmental stages, and time, we evaluated each variable individually using pairwise PERMANOVA testing with Bonferroni corrections. To perform these tests, we used the packages vegan (Oksanen et al. 2019) and pairwiseAdonis2 (Arbizu 2020).

First, we tested for differences among all experimental genotypes when at the vegetative stage to identify if flowering-time mutations could pleiotropically affect microbial community composition. Given that genotypic differences in flowering time arose from single-gene mutations, we anticipated little difference between, for instance, the two early-flowering genotypes (vp-32 and tflp2-2) or among any of the genotypes when all plants were in the same vegetative stage of development. If there was no evidence of pleiotropic effects of the single-gene mutations on the microbiome, then observed differences in microbial community could be attributed to either developmental stage or time point rather than artefacts of experimental host genotypes. Next, to test our hypothesis that microbial diversity changes over time potentially due to microbial community succession and regardless of flowering status (Fig. 1), we looked for differences in community composition and indicator taxa among time points. To test our second hypothesis, that plant developmental stage determines microbial assemblage in the rhizosphere (Fig. 1), we compared composition of rhizosphere microbial taxa among early-, wild type, and late-flowering phenotypes within specific time points. We categorized individuals by their stage at each sampling point (vegetative, flowering, fruiting, and senescing) and compared rhizosphere microbiomes within each time point. In addition to splitting samples by time since germination, we used nested PERMANOVA that allowed us to treat time as a blocking factor and control for it in our statistical analyses, that is, developmental stage was nested within time point.

To identify individual microbial taxa whose abundances changed, we conducted an indicator species analysis using DESeq2 (α = 0.01, Bonferroni corrected) (Love, Huber and Anders 2014) on genotypes, time points and developmental stages at the family level using rarefied data (Weiss et al. 2017). Again, samples were rarefied to 77213 reads. We made all possible comparisons and retained all taxa deemed an indicator in any of the comparisons. We then pruned the list of indicators to only those taxa appearing as an indicator for an individual group (e.g. indicators of plants at T2 or of flowering plants). Subsequently, we visualized shifts in the abundance of indicator taxa at the phylum, order, and class levels using the plot_bars function in phyloseq.

RESULTS

Microbial diversity in the rhizosphere of A. thaliana plants was similar on average across genotypes at time point one, with no significant differences in either α-diversity (P > 0.34, df = 4, Fig. 2A) or β-diversity (F_4,16 = 0.9601, P > 0.6634, R² = 0.1935, Fig. 2B). Using nested PERMANOVA analyses, we show that genotype was not a significant predictor of β-diversity when we controlled for the effects of developmental stage and sampling time (F_4,139 = 0.844, P = 0.297, R² = 0.0237). These results indicate that the single-gene mutations affecting flowering time did not pleiotropically affect microbial communities, and that any differences in microbial community composition arose from factors other than host genotype.

Figure 2.

Diversity across genotypes at T1, when all plants were vegetative. A minimum of three replicates per genotype were included in these analyses, for a total of 21 samples altogether. (A) Richness, Shannon, and Simpson diversities of all A. thaliana genotypes included in the experiment. Top and bottom lines of boxes indicate 25th and 75th quartiles, respectively. The middle line represents the median value for each α-diversity metric. Whiskers represent 1.5 x the interquartile range (IQR). There were no significant differences among A. thaliana genotypes at α = 0.05.(B) Nonmetric multidimensional scaling (NMDS) of Bray–Curtis dissimilarities of the bacterial community among genotypes. Genotype was non-significant in explaining variation in bacterial community structure as per adonis testing (p = 0.79). Ellipses represent 95% confidence intervals for the mean of each A. thaliana genotype. Points are colored by A. thaliana genotype.

Diversity changed over time (Fig. 3A and B), such that the first two time points (T1 and T2) were significantly less diverse than the later time points (T3 and T4) (p < 0.01, df = 3, Fig. 3A). Pairwise comparisons of bacterial α-diversity in the rhizosphere microbiome showed that T1 and T2 were not significantly different from one another nor were T3 and T4. For β-diversity, rhizosphere bacterial assemblages were significantly different from each other among all time points (F_3,140 = 17.36, p < 0.001, R² = 0.27, Fig. 3B). We found the same result in the nested PERMANOVA, with time since germination remaining statistically significant (F_3,140 = 17.36, p < 0.001, R² = 0.27).

Figure 3.

Diversity across time points. These analyses include a minimum of three replicates per genotype per timepoint, for a total of 144 samples. (A) Richness, Shannon, and Simpson diversities of individual sampling time points of the experiment. Top and bottom lines of boxes indicate 25th and 75th quartiles, respectively. The middle line represents the median value for each α-diversity metric. Whiskers represent 1.5 x the interquartile range (IQR). Bars above groups indicate no significant difference between those groups and letters indicate significant differences among sampling time points at α = 0.05.(B) Nonmetric multidimensional scaling (NMDS) of Bray–Curtis dissimilarities of the bacterial community among time points. Time since planting was significant in explaining variation in bacterial community structure as per adonis testing (P < 0.001). Ellipses represent 95% confidence intervals for the mean of each sampling time points. Points are colored by sampling time, and point shapes are based on A. thaliana genotype.

When subset by time point or analyzing with developmental stage nested within time point to specifically test the effect of plant developmental stage on α-diversity of the rhizobacterial community, we observed few differences. The effect of developmental stage on α-diversity was non-significant at T2 when plants of genotypes svp-32 and tfl2-2 were flowering, but the wildtype and genotypes ap1-10 and ft-10 were still vegetative (P > 0.3, df = 1). At T3, when the early-flowering genotypes svp-32 and tfl2-2 were fruiting, the wildtype was flowering, and the late-flowering genotypes ap1-10 and ft-10 were vegetative, rhizosphere Shannon diversity of flowering plants was somewhat lower than vegetative plants (p = 0.089, df = 2), but we observed no significant differences at α = 0.05. At T4, when genotypes svp-32 and tfl2-2 were senescing, the wildtype was fruiting, and the ap1-10 and ft-10 were flowering, there were no differences in α-diversity among developmental stages, but there was a significant effect of developmental stage on β-diversity (F_2,35 = 1.586, P < 0.05, R² = 0.083). Pairwise Adonis testing showed that the flowering and senescing stages were different at sampling time four (F_1,33 = 1.91, P < 0.01, R² = 0.056). In the related nested analysis of developmental stage within time point, inclusion of time since germination as a blocking factor resulted in a non-significant effect of developmental stage (F_3,140 = 4.757, P = 0.13, R² = 0.092).

The indicator analyses at the family level revealed zero taxa indicative of genotype, again consistent with the hypothesis that monogenic mutations affecting flowering time did not affect rhizosphere microbes. About 80 taxa were indicative of developmental stage, and 129 taxa indicative of time point. Of the 80 taxa indicative of phenological status, one was a unique indicator of the vegetative stage (phylum Bacteroidetes; Fig. S1, Supporting Information), five of the flowering stage (all phylum Proteobacteria; Fig. S2, Supporting Information), six of the fruiting stage (Proteobacteria, 5; Bacteroidetes, 1; Fig. S3, Supporting Information), and five of plants that were senescing (Proteobacteria, 2; Bacteroidetes, 3; Fig. S4, Supporting Information). The other 62 were shared by at least two different developmental stages (Table S1, Supporting Information). When plotting common indicator taxa (those with greater than 0.5% total relative abundance across all samples) by developmental stage, a more diverse array of indicator taxa were found in the later vs. earlier developmental stages (Fig. 4). Of the 129 taxa indicative of time point (Figs S5–S8, Supporting Information), 11 were uniquely indicators of T1, one was uniquely an indicator of T2, seven were uniquely indicators of T3 (Fig. S7, Supporting Information), and seven were uniquely indicators of T4 (Fig. S8 and Table S1, Supporting Information). Of the 11 indicator taxa for T1, eight were Actinobacteria, one Bacteroidetes, one Deinococcus-Thermus and one Proteobacteria. The lone indicator of T2 was of the phylum Proteobacteria. Of the seven indicators for T3, three were Actinobacteria, one Armatimonadetes, one Firmicutes and two Proteobacteria. Finally, of the indicator taxa unique to T4, one was Actinobacteria, three Bacteroidetes, one Chloroflexi, and two Proteobacteria. The other 103 taxa were shared between at least two time points. Plotting of common indicator taxa by time showed a shift between each time point (Fig. 5).

Figure 4.

Common indicator taxa, as determined through DESeq2 analysis (α = 0.05), with relative abundance > 0.05%, are plotted at the family level and split by developmental stage. All 144 samples were included in this analysis.

Figure 5.

Common indicator taxa, as determined through DESeq2 analysis (α = 0.05), with relative abundance > 0.05%, are plotted at the family level and split by sampling time. All 144 samples were included in this analysis.

DISCUSSION

Our objective was to identify whether variation in microbial diversity in the rhizosphere was best explained by developmental stage of host plants, or by the progression of time and the potential for microbial community succession since inoculation of the soil and emergence of host plants. We designed our experiment to disentangle these factors by selecting A. thaliana knockout mutants that have subtle genomic differences (i.e. at a single locus) but differ in flowering time, allowing us to observe plants in different developmental stages after a fixed amount of time since germination. Microbial diversity did not significantly differ among genotypes when all plants were in a common, vegetative stage. Nor did diversity differ among genotypes when we controlled for the effects of developmental stage and time since planting. Therefore, we conclude that observed differences were due to factors other than potential pleiotropy of flowering time mutations on the rhizosphere microbiome. We found strong support for the hypothesis that time since germination explains differences in microbial diversity, with diversity increasing over time regardless of host developmental status. In contrast, we did not find strong support for the hypothesis that developmental stage determines microbial diversity.

Microbial communities change over time in diverse ecosystems ranging from freezing freshwater lakes, which may shift in assemblage over months (Butler et al. 2019), to tropical air in which observable changes can occur in as few as two hours (Gusareva et al. 2019). Likewise, the rhizosphere microbiome has been shown to change with exudate chemistry (Guyonnet et al. 2018; Hu et al. 2018; Sasse, Martinoia and Northen 2018). While we did not sample root exudates in our experiment, changes in root exudate chemistry occur over the lifespan in A. thaliana (Chaparro et al. 2013) (either due to developmental stage or increasing size of plants with time), and such shifts were likely important in determining rhizosphere microbial diversity in the current study. Considering the existing literature, observing changes in rhizosphere diversity over the lifetime of a host plant is hardly surprising (e.g. Micallef et al. 2009; Chaparro, Badri and Vivanco 2014; Na et al. 2019). One study also used a genetic manipulation of phenology decoupled flowering time from time since germination; Dombrowski et al. (2017) and found that the rhizosphere microbiome of Arabis alpina changed with soil residence time but was unaffected by the onset of flowering. These results indicate that developmental stage may not be the key to changes in a rhizosphere microbiome. Observations in our experiment are consistent with Dombrowski et al. (2017), and expand those results by looking at a greater number of life-history transitions, namely onset of flowering and fruiting through senescence. We likewise find a greater effect of time rather than any developmental transition.

Separating the contribution of host plant developmental stage from time since germination was central to identifying which factor is a key determinant of rhizosphere microbial diversity. When we looked among developmental stages without sub-setting by time (data not shown), developmental stage appeared to be an important and significant predictor of rhizosphere microbial diversity. However, this result is potentially misleading because developmental stage was confounded by time; even though the timing of flowering differed across genotypes, the sequence of stages was necessarily the same for all, with vegetative growth followed by flowering which was followed by fruiting and then senescence. Thus, by examining the contribution of developmental stage within separate time points, our experimental design allowed us to assess the contributions of developmental stage independent of time since germination. In our models, when we controlled for the effect of time since germination, the influence of developmental stage on microbial beta-diversity became a non-significant predictor. This result, coupled with the fact that time since germination remained a significant predictor when we controlled for developmental stage, indicates that the major driver of microbial β-diversity in the rhizosphere of Arabidopsis thaliana is time since germination rather than host developmental stage.

Critically, the differences we observed were driven by time, and the few small differences attributable to developmental stage varied with diversity metric. While time is not separated from developmental stage in other studies (e.g. Rasche et al. 2006; Houlden et al. 2008; Chaparro, Badri and Vivanco 2014), such separation is challenging if not impossible if a study focuses on plants in a natural setting or involves the great majority of species for which relevant phenology mutants do not exist (e.g. Rasche et al. 2006; Houlden et al. 2008; Na et al. 2019). Our results do not dispute the value of other work studying changes in rhizosphere microbial communities, but do suggest that the seeming role of host developmental stage may well be attributable to time since host germination and opportunity for succession within the rhizosphere.

Biological interactions among microbes are important in steering microbial assemblages over time, as microbes depend on each other for a variety of metabolic functions (Zelezniak et al. 2015; Sanchez-Gorostiaga et al. 2019) and often form complex trophic cascades that result in an accumulation of diverse metabolic pathways and substrates (Gralka et al. 2020). Competition for resources among microbes can shift diversity or presence of specific taxa (Foster and Bell 2012; Gielda and DiRita 2012); in biofilms (comprising various bacteria, genetic strains, nutrient availability, and physical properties) interactions among cells affect competitive ability and fitness and thus microbial abundances (Nadell, Drescher and Foster 2016). We were interested in potential interactions among microbes in the rhizosphere over the course of the experiment, and in which microbes were present across developmental stages and sampling events. Our efforts to identify unique indicator taxa of either developmental stage or time only begin to tease apart these complex interactions.

The number of unique indictor taxa present at each time point followed the same pattern as diversity, with fewer taxa present at T1 and T2 than at T3 and T4. These parallel patterns show that as richness and diversity increased over time, no one or two taxa dominated the rhizosphere microbiome. A variety of taxa were present in similar proportions across all time points while other bacteria were unique indicators of each time point. These unique indicators could have key functions or ecological importance that correspond to activities among microbes present at those points. For example, at T1 when all plants were vegetative, unique bacteria included a member of Burkholderiaceae, a family that exhibits exceptional metabolic diversity (Staley 2007) and may be suppressive to fungal root pathogens (Carrión et al. 2018). Also an indicator of T1, Sphingomonadaceae cope well with stress and act as reservoirs of antibiotic resistance (Vaz-Moreira, Nunes and Manaia 2011). These characteristics could be particularly important for newly germinated plants. We could interpret the presence of these taxa, and general characteristics or functions with which they are associated, as explanations for their indicator roles at each time point; despite these possibilities, however, these hypothesized functional rolls would have to be supported in future experiments.

While rhizosphere microbial diversities did not change with developmental stage, abundance of individual microbial taxa did. We saw different indicator taxa across developmental stages, but because our identification of indicator taxa by developmental stage also included time (unlike our diversity analyses), we cannot say for certain whether time or developmental stage is responsible for the shifts in abundance in these taxa. The changing identities of indicator taxa with developmental stage makes sense as root exudates can change with developmental stage (Chaparro et al. 2013). It is possible that sampling all experimental plants at the latest plant developmental stage could reveal distinct indicators of stage independent of time. For example, only the two early-flowering genotypes reached the point of senescence, while the late-flowering genotype did not progress past the flowering stage. The most significant difference that we observed among developmental stages was between fruiting and senescing plants at the final sampling point. If we had carried this experiment to the point of senescence for all genotypes, perhaps we would have observed convergence in rhizosphere microbial diversities across all genotypes, as has been the case in other studies (Micallef et al. 2009). Even if the experiment had run long enough for communities to converge by plant senescence, if developmental stage were as strong a determinant of microbial indicators and diversity as was time we should have observed strong differences across developmental stages at each time point.

Microbial succession, the sequential and predictable assembly of microbial communities in an environment (Fierer et al. 2010), inherently means change. Primary succession in particular, wherein microbial colonization occurs in a sterile environment (or one close to sterile) such as glacial till (Nemergut et al. 2007), new lava beds (Kelly et al. 2014), or a developing human baby (Cong et al. 2016; Robertson et al. 2019), is the most extreme starting point for assembling microbial communities (Fierer et al. 2010). Succession in microbial communities may also occur following disturbance events that result in a shift in the dominant microbial taxa; extreme examples include after application of herbicides (Itoh et al. 2014; Schlatter, Burke and Paulitz 2018), following wildfire (e.g. Ferrenberg et al. 2013), or after bark beetle infestation (Mikkelson et al. 2017; Custer, van Diepen and Stump 2019). Even in environments with an established microbiome, fine-scale changes in microbial communities may occur in predictable patterns over time, such as changes in microbes during fermentation processes (Jeong et al. 2013; Sulaiman et al. 2014; Walsh et al. 2016). In our experiment, the microbial community experienced a coarse-grained disturbance via soil autoclaving, followed by establishment of a plant overstory. Succession may well account for shifts in microbial community composition that we observe.

Future studies should consider our findings when attempting to describe the importance of host developmental stage on rhizosphere diversity. By including a temporal control in future work, researchers will be able to assess the contribution of host developmental stage without over emphasizing its importance. Additionally, by concurrently testing for predictable changes over time, researchers may be able to identify patterns of succession among their study subjects.

DATA AVAILABILITY

Metadata, fastq files of sequence reads, and R code to process reads are publicly available at Pathfinder: https://doi.org/10.15786/14736912.v1. Metadata and R code to complete statistical analyses are publicly available on github: https://github.com/gcuster1991/Sucssion_Flowering_Time.

ACKNOWLEDGEMENTS

This research was supported by the Ronald E. McNair Post-Baccalaureate Achievement Program, National Science Foundation Grants [IOS-1444571], and the Microbial Ecology Collaborative with funding from National Science Foundation award [#EPS-1655726]. Particular thanks to Liz Nysson, Susan Stoddard, and Mallory Lai.

Conflicts of interests

None declared.