Assembly of the Populus Microbiome Is Temporally Dynamic and Determined by Selective and Stochastic Factors

ABSTRACT Recent work shows that the plant microbiome, particularly the initial assembly of this microbiome, influences plant health, survival, and fitness. Here, we characterize the initial assembly of the Populus microbiome across ten genotypes belonging to two poplar species in a common garden using 16S rRNA gene and ITS2 region amplicon sequencing of the leaf endosphere, leaf surface, root endosphere, and rhizosphere. We sampled these microbiomes three times throughout the first growing season and found that the composition of the microbiome changed dramatically over time across all plant-associated habitats and host genotypes. For archaea and bacteria, these changes were dominated by strong homogenizing selection (accounting for 29 to 62% of pairwise comparisons). However, fungal assembly was generally characterized by multiple ecological assembly processes (i.e., a mix of weak selective and dispersal processes). Interestingly, genotype, while a significant moderator of microbiome composition, generally explained less variation than sample date across plant-associated habitats. We defined a set of core genera that accounted for, on average, 36% of the microbiome. The relative abundance of this core community was consistent over time. Additionally, using source tracking modeling, we determined that new microbial taxa colonize from both aboveground and belowground sources, and combined with our ecological assembly null models, we found that both selective and dispersal processes explained the differences between exo- (i.e., leaf surface and rhizosphere) and endospheric microbiomes. Taken together, our results suggest that the initial assembly of the Populus microbiome is time-, genotype-, and habitat-dependent and is moderated by both selective and stochastic factors. IMPORTANCE The initial assembly of the plant microbiome may establish the trajectory of forthcoming microbiome states, which could determine the overall future health of the plant. However, while much is known about the initial microbiome assembly of grasses and agricultural crops, less is known about the initial microbiome of long-lived trees, such as poplar (Populus spp.). Thus, a greater understanding of initial plant microbiome assembly in an ecologically and economically important plant such as Populus is highly desirable. Here, we show that the initial microbiome community composition and assembly in the first growing season of Populus is temporally dynamic and is determined by a combination of both selective and stochastic factors. Our findings could be used to prescribe ecologically informed microbial inoculations and better predict the composition of the Populus microbiome into the future and to better understand its influence on plant health.

numerous distinct genotypes and their genomes have been described (29). These genotypes vary not only phenotypically (e.g., resistance to fungal pathogens differs among genotypes of the same poplar species [23]), but also in their microbiome composition (7,30,31). Genotypic influences on the plant microbiome are common across species due to differences in metabolites, elemental concentrations, and root and leaf traits (11), and these influences are generally stronger in aboveground, endospheric microbiomes (15). However, the degree to which underlying assembly mechanisms influence these microbiome compositional differences across poplar genotypes is currently unknown.
A common goal of plant microbiome research is to define and characterize the core microbial community (12,14,17,(32)(33)(34)(35). Core plant microbiomes generally represent a small proportion of microbial richness but are consistently present among samples at high relative abundances (33). While it is debated whether the core microbiome confers benefits to the host (12,33,34), reducing the complexity of the microbiome presents a greater opportunity to investigate the impact of individual microbial members on host fitness (32,36). Hence, developing an understanding of the core microbiome of poplar trees could lead to future targeted microbiome interventions (i.e., synthetic communities) that promote plant health and productivity. Indeed, the addition of a synthetic community consisting of Burkholderia and Pseudomonas strains isolated from Populus hosts increased lateral root formation and root hair production in Arabidopsis plate assays, and these organisms are predicted to carry out different functions related to growth and plant growth promotion in Populus (37).
To increase our understanding of the initial assembly of the Populus microbiome across genotypes, we planted ten Populus genotypes consisting of P. deltoides and P. trichocarpa in a common garden and sampled the leaf endosphere, leaf surface, root endosphere, and rhizosphere microbiomes three times throughout the first growing season-directly before field propagation (T0, 16 May 2017), 29 June 2017, and 18 September 2017. We used 16S rRNA gene and ITS2 region amplicon sequencing to characterize the archaeal and bacterial and the fungal communities, respectively. We hypothesized that microbiome community composition as well as the underlying assembly mechanisms (e.g., deterministic versus stochastic assembly) would change significantly during this time. Specifically, we hypothesized that deterministic selection would increase over time, as has been demonstrated in other successional contexts (38). Furthermore, we expected greater deterministic selection in plant endosphere microbiomes compared to the rhizosphere and leaf surface microbiomes due to the selective pressures of plant defense compounds (39). Overall, we expected the development of a dominant core community during the first growing season. We also attempted to resolve the colonization source of the root and leaf endosphere. We hypothesized that the rhizosphere would be the greatest source of root endosphere microorganisms, and the leaf surface would be the greatest source of leaf endosphere microorganisms due to spatial proximity. Our overall goals were to characterize the initial microbiome assembly across plant-associated habitats (i.e., plant tissues and immediate environment) and to determine the colonization source of endospheric microorganisms.

RESULTS
a-Diversity changes over time in multiple plant-associated habitats. We detected a significant three-way interaction in the effect of plant-associated habitat, genotype, and sample date on archaeal and bacterial a-diversity expressed as Hill numbers (40) at both q (defined in Materials and Methods) = 0 (analogous to richness) and q = 1 (analogous to Shannon's entropy) (analysis of variance [ANOVA]q = 0: F 42,343 = 1.921, P = 0.001; q = 1: F 42,343 = 2.622, P , 0.001). When analyzed separately among plant-associated habitats, we generally found that sample date only affected a-diversity in the exospheric microbiomes (i.e., leaf surface, rhizosphere) for both q = 0 and q = 1 ( Fig. 1; see Fig. S1 at https://doi.org/10.6084/m9.figshare.14251463.v1). While the leaf surface decreased in a-diversity over time, the rhizosphere increased over time. Alternatively, endospheric microbiome a-diversity was relatively stable over time. At q = 0, there was a significant effect of genotype on both the leaf endosphere and rhizosphere a-diversity, but at q = 1, the effect of genotype was only significant in the leaf endosphere (see Table S1 and Fig. S2 at https://doi.org/10.6084/m9.figshare.14251463.v1). However, for both q = 0 and q = 1, the effect of genotype interacted with sample date for leaf endosphere a-diversity such that there was only a genotype effect during the latest September sample date (see Table S1 and Fig. S2 at https://doi.org/10.6084/m9 .figshare.14251463.v1).
Taxonomic succession, for the most part, varied among plant-associated habitats. For instance, in the leaf endosphere, there was a relative increase in Proteobacteria over time (particularly Gammaproteobacteria), while in the rhizosphere, the relative abundance of Proteobacteria declined (see Fig. S4A at https://doi.org/10.6084/m9 .figshare.14251463.v1). Similarly, at the order level, the leaf surface and rhizosphere followed diverging patterns: the relative abundance of Sphingomonadales increased in Genotype was also a significant moderator of archaeal and bacterial community composition in all but the leaf surface microbiome (Table 1). Genotype effects were  strongest in the two endosphere microbiomes, and in these microbiomes, the effect of genotype was stronger than that of sample date (Table 1). However, in the leaf endosphere, there was a significant sample date-genotype interaction (PERMANOVA, P , 0.001), where there was only a genotype effect during the June and September sample dates (T0, P = 0.261; June, P = 0.021, R 2 = 0.407; September, P = 0.003, R 2 = 0.569). Interestingly, genotypic differences in the belowground microbiomes were intraspecific, as host species was not a significant moderator of community composition (root endosphere, P = 0.440; rhizosphere, P = 1.000; Fig. 2A). However, in the latter two sample dates for the leaf endosphere microbiome (when genotype was a significant moderator of community composition), there was a significant difference between communities in P. deltoides and P. trichocarpa (P = 0.026, R 2 = 0.031). This difference in the microbiome composition was noted at the genus level with a greater relative abundance of Streptococcus spp. in the P. deltoides leaf endosphere (see Fig.  S6A at https://doi.org/10.6084/m9.figshare.14251463.v1). Upon further multiple-comparison testing among genotypes, we failed to detect significant differences between genotypes within the same plant-associated habitat across all sample dates (P , 0.05; see Table S2 at https://doi.org/10.6084/m9.figshare.14251463.v1). Plant-associated habitat explained 24.1% of the variation in fungal community composition (PERMANOVA, P , 0.001, R 2 = 0.241). Similar to the archaeal/bacterial results described above, because there was a three-way interaction among plant-associated habitat, genotype, and sample date (P , 0.001), we analyzed each plant-associated habitat separately. When separated by plant-associated habitat, sample date was consistently a significant moderator of fungal community composition and, on average, explained 22% of the variation ( Fig. 2B; full statistics in Table 1). Taxonomically, this was represented as a gradual increase in the relative abundance of Basidiomycota over time across all plant-associated habitats (see Functionally, there was a shift in the root endosphere fungal community away from pathogenicity at T0 toward mycorrhization during the June and September sample dates (Fig. 3). Accordingly, at T0, fungal pathogens constituted over 60% of the fungal reads in the root endosphere, while mycorrhizal fungi were almost nonexistent (,0.02%). By the September sampling date, the relative abundance of pathogens decreased to 16% of the root endosphere fungal reads, respectively, while ectomycorrhizal (EM) fungi increased to 7%. Arbuscular mycorrhizal (AM) fungi also increased dramatically over the sampling period, although relative abundances were significantly smaller than those of both pathogens and ectomycorrhizal fungi (Fig. 3).
Genotype was also a significant moderator of fungal community composition in all but the root endosphere microbiome (Table 1). However, in the leaf endosphere and surface microbiomes, there were significant sample date-genotype interactions (PERMANOVAleaf endosphere, P = 0.032; leaf surface, P , 0.001), where the effect of genotype was strongest in the later sampling dates (although not significant for the leaf endosphere; see Table S3 at https://doi.org/10.6084/m9.figshare.14251463.v1). Genotypic differences could only be attributed to host species in the leaf surface microbiome (P = 0.030, R 2 = 0.019), where P. trichocarpa had greater relative abundances of the plant pathogen Alternaria spp. (see Fig. S6B at https://doi.org/10.6084/m9 .figshare.14251463.v1).
Microbial assembly processes differ by microbial domain. Using a null modeling approach, we assessed ecological assembly processes of the microbial communities following Stegen et al. (41). Assembly processes were broadly categorized as (i) variable selection, whereby selective processes lead to disparate microbial communities, (ii) homogenous selection, whereby selective processes lead to similar microbial communities, (iii) dispersal limitation, whereby limitations to dispersal allow ecological drift to lead to disparate microbial communities, and (iv) homogenizing dispersal, whereby high rates of dispersal lead to similar microbial communities (see Materials and Methods for complete statistical criteria).
Across all plant-associated habitats, no single assembly process dominated archaeal and bacterial community assembly (Fig. 4A). However, when assessed for each sample date and plant-associated habitat individually (the two strongest moderators of community composition), homogenous selection was the primary driver of between-community shifts in composition for 29 to 62% of pairwise comparisons within each plant-associated habitat-sample date combination (Fig. 4B). The relative dominance of homogenous selection as the primary assembly process varied with both plant-associated habitat and sample date (Fig. 4B). For example, homogenous selection decreased during the growing season in aboveground microbiomes but remained fairly constant in the belowground archaeal and bacterial microbiomes (Fig. 4B). Other assembly processes, such as dispersal limitation, also varied by both plant-associated habitat and sample date. While dispersal limitation increased from T0 to June and decreased from June to September in the leaf endosphere and rhizosphere, it consistently increased in the root endosphere and leaf surface archaeal and bacterial microbiome (Fig. 4B).
Interestingly, variable selection was also a relatively important moderating factor in the root endosphere throughout the first growing season, which would appear to contradict the dominance of homogenous selection. However, pairwise comparisons for the same genotype within each plant-associated habitat-sample date combination were often 100% dominated by homogenous selection (see Table S4 at https://doi .org/10.6084/m9.figshare.14251463.v1). Thus, significant variable selection among genotypes occurred in the root endosphere, particularly among P. trichocarpa genotypes (see Table S5 at https://doi.org/10.6084/m9.figshare.14251463.v1). Differences between the leaf endosphere and leaf surface and the root endosphere and rhizosphere were characterized by a combination of dispersal limitation and variable selection, while similarities were mostly characterized by homogenous selection, and these patterns were fairly consistent over time (see Fig. S7A and B at https://doi.org/10.6084/m9.figshare.14251463.v1).
Across all plant-associated habitats, sample dates, and genotypes, the fungal community was characterized by both dispersal limitation and undominated assembly processes, with selection playing a relatively minor role (Fig. 4C). When assessed for each sample date and plant-associated habitat individually, fungal community assembly was, for the most part, undominated by selective or dispersal processes (Fig. 4D). However, similar to the archaeal and bacterial communities, the dominance of assembly processes differed by sample date. In the leaf surface fungal microbiome, there was a shift from homogenizing dispersal to homogenous selection during the first growing season (Fig. 4D). Alternatively, dispersal processes, at the expense of homogenous selection, became relatively more dominant in the rhizosphere for fungi during the growing season. Unlike for archaea and bacteria, pairwise comparisons within plant-associated habitat-sample date-genotype combinations for fungi were most often characterized by assembly processes undominated by selection or dispersal (see Table S6 at https://doi.org/10.6084/m9.figshare.14251463.v1), a pattern that was consistent within host species as well (see Table S7 at https://doi.org/10.6084/m9 .figshare.14251463.v1). Differences between the leaf endosphere and leaf surface were mostly controlled by weak selection and dispersal factors (i.e., "undominated"), while differences between the root endosphere and rhizosphere were characterized by a combination of undominated assembly processes and dispersal limitation (see Fig. S7C and D at https://doi.org/10.6084/m9.figshare.14251463.v1). These patterns were consistent over the growing season.
Core taxa disproportionately represent the Populus microbiome. Using occupancy-abundance distributions of microbial genera across all time points, we detected 23 distinct bacterial (Table 2) and 12 distinct fungal (Table 3) core genera across the four plant-associated habitats based on .95% occupancy and .1% average relative abundance (see Fig. S8   detected no core archaeal and bacterial genera in the leaf endosphere, on average, archaeal and bacterial core genera constituted 51%, 24%, and 49% of reads in the leaf surface, root endosphere, and rhizosphere microbiome, respectively. The core genera accounted for 1.4%, 0.8%, and 2.9% of archaeal and bacterial genera in the leaf surface, root endosphere, and rhizosphere microbiome, respectively. Across the growing season, the relative abundance of the core archaeal and bacterial community increased in the leaf surface (ANOVA, F 2,102 = 65.74, P , 0.001), decreased in the rhizosphere (F 2,103 = 13.11, P , 0.001), and did not significantly change in the root endosphere (F 2,102 = 1.04, P = 0.357; Fig. 5A).
For fungi, core genera constituted, on average, 44%, 44%, 7%, and 67% of reads in the leaf endosphere, leaf surface, root endosphere, and rhizosphere fungal microbiome, respectively. The core genera accounted for 0.7%, 1.4% 0.4%, and 1.8% of fungal genera in the leaf endosphere, leaf surface, root endosphere, and rhizosphere microbiome, respectively. Similar to archaea and bacteria, the relative abundance of fungal core genera across the growing season differed by plant-associated habitat. During the first growing season, the relative abundance of fungal core genera increased in the leaf endosphere (ANOVA, F 1,28 = 14.06, P , 0.001), remained the same in the leaf surface (F 2,103 = 0.38, P = 0.685), and decreased in the root endosphere (F 2,58 = 9.46, P , 0.001) and rhizosphere (F 2,103 = 105.1, P , 0.001; Fig. 5B). Half of the core fungal genera were classified as potential plant pathogens, while the other half were classified as saprotrophs. Interestingly, many core genera were shared among plant-associated habitats. For instance, the bacterium Sphingomonas and the fungal plant pathogen Fusarium were core genera in the leaf surface, root endosphere, and rhizosphere. Similarly, in addition to the aforementioned genera, Rhodoplanes and Streptomyces were core genera in both the root endosphere and rhizosphere. Additionally, the two core genera in the leaf endosphere, the fungal plant pathogens Alternaria and Boeremia, were also leaf surface core genera.
Endospheric microbiomes have both above-and belowground sources. Across the sink microbiomes (e.g., leaf and root endosphere), we were able to determine the source for 6.9% and 3.5%, on average, of new archaeal and bacterial and fungal taxa, respectively (see Table S8 at https://doi.org/10.6084/m9.figshare.14251463.v1) using SourceTracker (42). Generally, source microbiome contributions to the sink microbiome were comparable, with each source contributing about half of the source-attributable  new taxa within each sink microbiome (ANOVA, F 1,124 = 8.165, P = 0.004). The singular exception was for the leaf endosphere archaeal and bacterial microbiome. Here, while model predictions estimated that the root endosphere contributed almost 12% of new taxa to the archaeal and bacterial leaf endosphere community, the leaf surface was estimated to contribute only 2.6%. Source contributions were also consistent over time for both the leaf and root endosphere across both amplicons (16S-leaf, P = 0.096; ITSleaf, P = 0.440; 16S-root, P = 0. 213; ITS-root, P = 0.159).

DISCUSSION
The initial assembly of the plant microbiome may establish the trajectory for forthcoming microbiome states (via priority effects) and could determine the overall future health of the plant (12,43,44). Thus, a greater understanding of initial plant microbiome assembly in an ecologically and economically important plant such as Populus is highly desirable (25,26). Here, we show that the initial microbiome composition and assembly in the first growing season of Populus are temporally dynamic and are determined by a combination of both selective and stochastic factors.
Sample date was consistently a strong predictor of the microbial community composition when assessed within each individual plant-associated habitat, which highlights how changing plant, soil, and climactic factors over time influence microbial community composition. First of all, increasing root and leaf biomass during the first growing season may affect the diversity of the microbial community through speciesarea relationships (45,46). Additionally, seasonal internal redistribution of sugar and nitrogen among plant tissues (47,48) would affect substrates necessary for microbial growth leading to changes in microbial community composition. Our finding that Gammaproteobacteria gradually replace Alphaproteobacteria in aboveground sample types over time contradicts previous studies in grass microbiomes (14,49), highlighting differences between tree and grass microbiomes. Additionally, as plants grow during the first growing season, exudation of plant photosynthates into the rhizosphere from roots likely increases (50). This increase in exudation has been shown to favor Proteobacteria and Bacteroidetes at the expense of Actinobacteria in Arabidopsis (50), Avena (51), and Citrus (52). Interestingly, our results showed opposite patterns (see Fig.  S4 at https://doi.org/10.6084/m9.figshare.14251463.v1), which could reflect other environmental changes that would favor Actinobacteria. Increased temperatures during the summer months could positively impact Actinobacteria (53). Furthermore, increased activity during summer months (54) may enhance nutrient availability (55,56) and shift the microbial community composition. These results contribute to the growing literature suggesting that microbial communities are temporally dynamic (57). However, we build upon this to show that, across habitat types (e.g., leaf endosphere, rhizosphere, etc.), the initial community assembly of the plant microbiome cannot be captured through examination of a single time point and that stochastic processes as well as deterministic ones shape the plant microbiome.
The effect of sample date was stronger in exospheric microbiomes (i.e., the leaf surface and rhizosphere) than in endospheric microbiomes (i.e., the leaf and root endosphere) for archaea and bacteria. This is corroborated by multiple studies that show that the host effects dominate the endosphere microbiome, while environmental factors influence the exospheric microbiomes to a greater extent (12,15,17). However, host effects likely also changed during the initial growing season as leaf area increased and patterns of nutrient and sugar distribution changed (47,48). This suggests that environmental changes over the first growing season are stronger than changes to the host. The muted community change with season in endospheric microbiomes could also be due to dispersal limitation, where microbes adapted to the changing exospheric conditions are prevented from colonizing the plant endosphere. Indeed, we found evidence for dispersal limitation to be a major factor in structuring differences between endo-and exospheric microbiomes (see Fig. S7 at https://doi.org/10.6084/m9 .figshare.14251463.v1) and within the endospheric microbiomes (Fig. 5). However, this finding that exospheric archaeal and bacterial microbiomes were influenced by sample date to a greater degree than in endospheric microbiomes was not followed by fungi (Table 1). This difference could be attributed to the overall greater stochasticity or weak selective/dispersal factors found in fungal community assembly (Fig. 5). Such differences highlight the need to study the plant microbiome holistically, as different microbial domains and plant-associated habitats can follow distinct patterns.
An increased understanding of fundamental community assembly processes describes the elementary units that underpin microbial community composition. Consistent with our hypothesis, the relative dominance of assembly patterns (e.g., variable selection, homogenous selection, dispersal limitation, and homogenizing dispersal) changed over time. However, a general pattern with time across plant-associated habitats for both archaea and bacteria and fungi did not emerge. This is surprising because in other successional contexts (e.g., soils after wildfire disturbances and a salt marsh soil chronosequence), selective assembly processes tend to dominate over time (58,59). While disturbances may lead to historical contingencies that could affect assembly mechanisms (9), the environmental pressures on the microbiome of a growing plant could also lead to large changes in assembly. Indeed, during the first 4 months of sorghum (Sorghum bicolor) growth, the relative dominance of community assembly processes of the fungal phyllosphere microbiome shifted from stochasticity (i.e., weak selective and dispersal processes) to homogenous selection and homogenizing dispersal (60). Furthermore, evidence of deterministic selection increased during successive microbial passages in the Solanum phyllosphere (61). These conflicting results among studies suggest that temporal patterns of microbial assembly may be idiosyncratic and depend on the host species. This demonstrates the need to understand plant microbiome assembly in multiple species aside from "model" grass systems.
Our finding that homogenous selection was the dominant assembly process for archaea and bacteria at each sample date across plant-associated habitats suggests that the initial Populus archaeal and bacterial microbiome assembly is repeatable and predictable even among different Populus genotypes. This is remarkable given the fact that this consistency was noted even when compared between two different Populus species. This is not to say that there were no differences among genotypes; however, the impact of genotype on microbial community composition across microbial domains was relatively similar or less than that of sample date. In the leaf and root endosphere, the archaeal and bacterial microbiome was, in fact, more structured by genotype than sample date. These genotypes vary in production of numerous metabolites, such as salicylic acid (7,62,63), which can correlate with changes in the microbial community composition (64). Plant phenotypic differences are even more pronounced across Populus species such as P. trichocarpa and P. deltoides, which have differential resistances to plant pathogens such as Sphaerulina musiva and Marssonina brunnea (31,65). Differences in microbiome composition among genotypes may grow stronger, especially within the endosphere, as the plants mature and differences in pathogen infectivity begin to emerge. For example, the effects of the leaf spot and stem canker fungal pathogens Marssonina brunnea and Sphaerulina musiva are known to dramatically increase over the initial years after establishment in plantation settings (66), and similar patterns have been observed at our site in P. trichocarpa trees in the subsequent years after this study (M. Cregger and W. Muchero, unpublished data). However, we did not find differences in pathogen relative abundance among different genotypes or Populus species ( Fig. 3; Fig. S6 at https://doi.org/10.6084/m9.figshare.14251463.v1), which contradicts early findings in mature Populus trees (31), further supporting out hypothesis that genotypic differences likely emerge as trees age.
With a few exceptions, our hypothesis that the dominance of a microbial core community would increase over time was unsupported by the data. Only the archaeal and bacterial leaf surface microbiome and the fungal leaf endosphere core microbiomes increased in dominance over time. However, across all time points, these 23 distinct bacterial and 12 distinct fungal core genera constituted a substantial proportion of the plant microbiome (7 to 69%). This supports earlier findings that overall beta diversity in the plant microbiome is controlled by the variation in a relatively small number of taxa (4,14,32,67). This is corroborated by the fact that most of the fungal core genera in our study were considered potential plant pathogens, which can disproportionately affect the plant microbiome (36). Interestingly, although Populus plants form both EM and AM symbioses (68), none of the core fungal genera were classified as such. This is likely because at T0, mycorrhizal associations had not yet established. Indeed, the proportion of both EM and AM reads increased during the first growing season in the root endosphere and rhizosphere, which suggests that mycorrhizal colonization in the field can occur relatively rapidly (,1 month), and there is a strong preference of such a relationship with time.
The bacterial core microbiome consisted of several taxa that have also been classified as "core" in numerous other studies across vastly different plant types. For example, core bacterial genera in the leaf surface in our study included Methylobacterium, Sphingomonas, Pseudomonas, and Hymenobacter, which have been classified as core members of the Panicum virgatum, Miscanthus x giganteus, Arabidopsis thaliana, and Medicago truncatula phyllospheres (14,69,70). Additionally, in the rhizosphere, core Alphaproteobacteria genera in our study, such as Bradyrhizobium, Burkholderia, Mesorhizobium, and Rhodoplanes are also considered core rhizosphere bacteria across numerous plant phyla (35). This suggests that there is a consistent consortium of bacteria that frequently colonize different plant-associated habitats regardless of plant type or species. In other words, these core genera may be selected based on the plantassociated habitat (e.g., phyllosphere versus rhizosphere), and differences in plantassociated habitats among plant species may be less important. It is possible then that these bacteria perform important functions in the plant. However, in contrast to the fungal core genera, the functions of the bacterial core genera still remain unclear because metagenome prediction based on amplicon data has yet to be verified by large-scale genomic sequencing of isolates or shotgun metagenomics in plant tissues. Until then, the functional significance of the core microbiome is still debated (33,34); however, these core members represent taxa that should be prioritized for isolation and subsequent genomic sequencing to understand their functional significance.
Consistent with previous studies, we found that the colonization of the root and leaf endosphere likely occurs from multiple sources (71)(72)(73). However, our hypothesis that the adjacent exospheric microbiome would be the dominant source of the endospheric microbiome (e.g., the rhizosphere for the root endosphere) was largely unsupported. Our results suggest that under most circumstances new taxa colonize from a combination of aboveground and belowground source environments. This finding somewhat contradicts previous studies that suggest that the soil microbiome is the dominant reservoir of plant endosphere taxa (14,70). While it is true that highly diverse soil microbial communities, in general, share many microbial members with those in the plant endosphere, a mechanistic model of microbial colonization cannot necessarily be inferred from such observations. For example, aeolian-derived microbes from different environments may be deposited on leaf as well the soil surface (74,75). Such microbes have been shown to enter the leaf endosphere through stomatal openings (73), so even though there is shared membership between the leaf endosphere and the soil, this does not necessarily mean that the soil is the source of said microbe. Our source tracking models attempt to remedy this by including multiple aboveground and belowground sources. Nevertheless, the lack of source attribution (,10% in many cases) limits the conclusions that can be made from these models. However, we show that it is likely that multiple sources contribute to plant endosphere colonization.
Extending these source tracking predictions to our community assembly models, we show that colonization of the plant endosphere by archaea and bacteria is influenced by both variable selection and dispersal limitation. This suggests that differences in the community composition between endo-and exo-spheric microbiomes is due to both selective pressures and physical limitations of colonization. While this has been previously hypothesized (4,73), we are the first to demonstrate with community assembly models that this is indeed the case. Furthermore, we show that, for fungi, stochastic factors likely play a major role in influencing the colonization of the plant endosphere from its adjacent exosphere, and in the root endosphere, dispersal limitation is also a factor. The mode of microbial colonization is an important consideration of plant microbiome assembly, and these models of microbial colonization (i.e., source tracking and assembly processes) should be paired with future experiments that empirically test modeled results.
Our results comprehensively characterize the early assembly of the temporally dynamic Populus microbiome among multiple plant genotypes and species and across different plant-associated habitats. While the fundamental assembly processes and colonization sources remained largely similar throughout the first growing season, the microbiome community composition shifted dramatically. This suggests that changing seasonal factors as well as changes associated with plant growth are primarily responsible for these shifts. Future work, such as long-term, seasonal plant microbiome studies, may be able to partition the influence of season and plant growth. It is hypothesized that the initial assembly of the microbiome "sets the stage" for future microbiome states through priority effects (76). Hence, where differences exist (e.g., the modest differences among plant genotypes), our findings could be used to predict the composition of the Populus microbiome into the future and better understand its influence on plant health.

MATERIALS AND METHODS
Plant collections. During December of 2016, eight P. trichocarpa genotypes were collected from Corvallis and Clatskanie genome-wide association study (GWAS) populations (65), and two P. deltoides genotypes were collected from a previous University of Tennessee Institute of Agriculture (UTIA) Sun Grant experiment in Blount County, TN, in the same field as the plots described below (31,77). Cuttings were kept on ice, shipped overnight, and maintained at 4°C until propagation (March 2017). Cuttings were rinsed in a 1% Zerotol 2.0 solution for surface sterilization, rooting powder was placed on sterile cutting surfaces (0.1% indole-3-butyric acid), and plants were placed in sterile potting soil (Fafard 52 Mix, Sun Gro Horticulture, Massachusetts, USA). Once significant root growth took place and leaf buds opened, cuttings were transferred to an experimental plot in Blount County, TN, managed by the UTIA-East Tennessee Research and Education Center (ETRED) (see Cregger et al. [31] for site details).
Experimental design, harvest, and sample processing. Soils were augured to transplant cuttings and minimize the effects of soil compaction on belowground root establishment. We planted six rows of 80 individuals spaced 1 m apart. These rows were broken up into three blocks made up of 4 to 26 replicates of each genotype. Sampling was conducted on three sampling dates in the initial growing season-day of transplant T0 (in May of 2017, representing the microbial community that was present in the cuttings and that colonized the plants in the greenhouse), June 2017, and September 2017 to assess initial microbiome community assembly across four plant-associated habitats (leaf endosphere, leaf surface, root endosphere, and rhizosphere). At each time point, three random replicate plants (one from each block) were destructively harvested per genotype, and samples were transferred to the laboratory. Roots, and the attached soil (operationally defined as rhizosphere soil), were stored at 280°C until root washing, sterilization, and genomic DNA (gDNA) extractions took place (within 2 months of sampling). Briefly, fine roots (,2 mm diameter) were sorted and surface-sterilized by sequential washing with bleach (3.125%) and then ethanol (70%) and rinsing of the roots with autoclaved water (4 times) as previously described (31); the rhizosphere was collected from an initial rinse with sterile water. Leaves were stored at 4°C after harvest and processed within 3 days. Leaves were rinsed (rinsate was collected as leaf surface samples) and surface-sterilized as previously described (31) by washing leaves with bleach and rinsing the leaves with autoclaved water (4 times). Postprocessing, leaf samples were also stored at -80°C until DNA extractions.
DNA extractions and Illumina MiSeq preparation and sequencing. Prior to extraction, root and leaf tissues were cut into fine pieces (;5 mm or less), and leaf and root rinsates (rhizosphere) were centrifuged at 10,000 Â g, and the supernatant was removed. These pelleted soil rhizosphere samples were then extracted using the PowerSoil DNA kit (Qiagen, Venlo, The Netherlands) following the standard protocol except that a Precellys tissue homogenizer was used to bead-beat extractions (30 sec of 5,500 Â g beadbeating with a 30 sec rest in triplicate). Root and leaf endosphere samples were extracted using the Qiagen PowerPlant Pro DNA kit following the standard protocol except that prior to extraction, 50 mg of tissue per extraction was bead-beaten for 1 min in liquid nitrogen blocks with one sterile steel bead two times. Extractions were quantified on a Nanodrop 1000 spectrophotometer (NanoDrop Products, Wilmington, DE, USA). Root and leaf endosphere extractions were purified and concentrated using a DNA Clean and Concentrator-5 kit (Zymo Research Corporation, Irvine, CA, USA) prior to PCR amplification.
A two-step PCR approach was used with barcode-tagged and frameshifting nucleotide primers targeting the 16S rRNA gene for archaea and bacteria and the ITS2 region for fungi (78) using pooled primer sets to increase coverage of archaeal, bacterial, and fungal taxa (31) (see Table S9 at https://doi .org/10.6084/m9.figshare.14251463.v1). The first step of PCR included 2.5 mM peptide nucleotide acid (PNA) blockers for 16S rRNA amplifications, and 2.5 mM PNA targeting plant nuclear rRNA genes for ITS2 region was used to reduce the amplification of plant material (see Table S9 at https://doi.org/10.6084/ m9.figshare.14251463.v1). Thermal cycler conditions for the primary PCRs for soils were 5 cycles of 95°C for 1 min, 50°C for 2 min, and 72°C for 1 min. Primary PCR conditions for plant tissues were 5 cycles of 95°C for 1 min, 78°C for 5 s, 50°C for 2 min, and 72°C for 1 min. Primary PCR products were cleaned with 17 ml of Agencourt AMPure beads and eluted in 21 ml of nuclease-free water. Secondary PCRs had purified DNA tagged with barcoded reverse primers and forward primers (see Table S9 at https://doi.org/10 .6084/m9.figshare.14251463.v1) in the 50-ml reaction mixture, except with 20 ml of purified DNA from primary PCRs. Thermal cycler conditions for secondary soil PCRs consisted of denaturation at 95°C for 45 s, followed by 32 cycles of 95°C for 15 s, annealing at 60°C for 30 s, 72°C for 30 s, and final extension at 72°C for 30 s. Secondary PCRs for plant tissue consisted of denaturation at 95°C for 45 s, followed by 32 cycles of 95°C for 15 s and 78°C for 5 s, with remaining cycle parameters the same as with soil secondary PCRs. After PCRs, all experimental units were pooled based on band intensity and purified with Agencourt AMPure XP beads (0.7:1 bead to DNA ratio; Beckman Coulter, Inc., Pasadena, CA, USA). Illumina MiSeq sequencing was carried out using a 9-pM amplicon concentration with a 15% PhiX spike (2 Â 300 cycles).
Bioinformatics. Paired-end .fastq sequence files had primers and adaptors removed individually using the Cutadapt program (79). Trimmed .fastq files were then joined using QIIME 1 and imported into the QIIME 2 environment (80). Joined sequences were demultiplexed, and median Phred quality scores were visualized. Due to poor quality on the 39 end for 16S data, they were truncated to 190 bp. The ITS reads were not further trimmed or truncated. Both 16S and ITS2 data sets were denoised and delineated into amplicon sequence variants (ASVs) using the DADA2 algorithm in QIIME 2 (81). Representative sequences were then assigned a taxonomic classification using the naive Bayes classifier through the sklearn Python package for 16S rRNA sequences with the SILVA database (82), while ITS2 gene representative sequences were assigned a taxonomic classification in QIIME 1 using BLAST and the UNITE reference database (83). Contaminants (unassigned reads, mitochondria, chloroplasts for 16S; Protista, Chromista, Animalia, and Plantae reads for ITS2) were removed. Contaminants were a low percentage of total reads for bacteria (;9%) and fungi (;15%). Fungi were further classified into functional guilds using the FUNGuild database (84).
Quantifying microbial assembly processes (see below) relies on phylogenetic turnover (41). The phylogenetic tree for 16S reads was created in the QIIME2 environment using FastTree 2 (85) and was midpoint rooted. Because ITS genes do not align well across large phylogenetic distances, we used a "ghost tree" approach to create a phylogenetic tree for ITS2 reads (86). Briefly, the 18S SILVA database (82) was used to create a foundational tree, and class-level ITS branches from the UNITE database (83) were grafted onto the ends of the foundational tree. The ITS2 ASVs were then clustered into 99% operational taxonomic units (OTUs) by closed reference OTU picking using VSEARCH (87). These ITS2 OTUs were only used to quantify microbial assembly processes; other analyses were conducted using ASVs.
Statistical analysis. All statistical analyses were conducted in R (88) using the Picante (89), phyloseq (90), and vegan (91) packages. Significance was determined at the a = 0.05 level for all statistical tests.
Differences in a-diversity over time and among genotypes were compared by means of Hill numbers (40) of samples rarified to 900 reads (for both amplicons) at orders of q = 0 and q = 1. The parameter q determines the relative weighting of rare species. At q = 0, all species are weighted equally (richness); at q = 1, species are weighted proportionally to their relative abundance (analogous to Shannon's index). Differences in means of Hill numbers among plant-associated habitats, sample dates, and genotypes were assessed by ANOVA. Where independent variables were significant, we assessed multiple comparisons by Tukey's test of honest significant differences. We used Q-Q plots and scale-location plots to inspect normality and homoscedasticity, respectively. Where these assumptions were unmet, independent variables were log-transformed and reanalyzed satisfying the aforementioned assumptions.
Differences in the microbial community composition among plant-associated habitats, sample dates, and genotypes were assessed by PERMANOVA (92), using Bray-Curtis distances applied to proportionally normalized data. We visualized differences in community composition using principal-coordinate analysis (PCoA).
Assembly processes were assessed using a null modeling approach following Stegen et al. (41). Assembly processes are broadly categorized as (i) variable selection, whereby selective processes lead to disparate microbial communities, (ii) homogenous selection, whereby selective processes lead to similar microbial communities, (iii) dispersal limitation, whereby limitations to dispersal allow ecological drift to lead to disparate microbial communities, and (iv) homogenizing dispersal, whereby high rates of dispersal lead to similar microbial communities. This approach uses both phylogenetic and taxonomic turnover to classify the dominant assembly process in pairwise sample comparisons. To characterize phylogenetic turnover, we use the between-community version of the (abundance-weighted) b-mean-nearest taxon distance (bMNTD) (93). Observed bMNTDs were then compared to a null model distribution of bMNTDs generated from 999 null model expectations (i.e., randomized taxa reshuffling among phylogenetic tree tips). The difference between observed bMNTD and the mean of the null distribution was measured in units of standard deviation (of the null distribution) and is referred to as the b-nearest taxon index (bNTI). Values of bNTI of $2 signify variable selection as the dominant assembly process, and bNTI values of #-2 signify homogenous selection as the dominant assembly process following Stegen et al. (41).
Taxonomic turnover was used to define the dominant assembly processes of jbNTI valuesj , 2. For taxonomic turnover, we use the Raup Crick index (94) modified to include species abundances (41). We calculated separate null models for archaea and bacteria and fungi using probability-based randomization. Null model microbial compositions were assembled for each sample by randomly sampling from the total ASV (or OTU for fungi, see above) pool in proportion to occupancy to define ASVs and then in proportion to abundance to define abundances into those selected ASVs, thereby maintaining similar levels of a-diversity. The community composition for each sample was probabilistically generated 999 times. For each iteration, the Bray-Curtis dissimilarity index between samples was calculated, and the proportion of iterations in which the index was smaller than or equal to the observed Bray-Curtis dissimilarity index between those pairs of samples was our resulting metric. We standardized this metric (RC Bray ) to range from 21 to 1 by subtracting 0.5 and multiplying by 2 (94). Values of jbNTI j , 2 and RC Bray $ 0.95 signify dispersal limitation, and values of jbNTI j , 2 and RC Bray # -0.95 signify homogenizing dispersal following Stegen et al. (41). Values of jbNTI j , 2 and jRC Bray j , 0.95 signify weak selection and moderate levels of dispersal such that no process dominates assembly and are therefore classified as "undominated" (95).
The core microbiome of each plant-associated habitat was characterized using occupancy-abundance distributions of microbial genera across all time points (14). Microbial genera were considered "core" when their occupancy (the proportion of presence in all samples) was at least 95% and their relative abundance was at least 1% averaged across all sampling dates within a plant-associated habitat.
We used SourceTracker (42) to determine the relative contribution of the rhizosphere and leaf surface as microbial sources of new microbial taxa (i.e., not present in the previous sample dates) in endospheric microbiomes for the June and August sampling dates. Differences in the relative contributions of sources among sample dates and genotypes were also assessed by ANOVA. These ANOVAs were analyzed using the same approach as for a-diversity.
Data availability. The data sets presented in this study can be found in the Sequence Read Archive under BioProject number 685817. The R code for all statistics and figures as well as the final ASV, taxonomy, and sample data tables used in this analysis can be found at https://github.com/nicholascdove/ Populus_microbiome_assembly.