Large-scale replicated field study of maize rhizosphere identifies heritable microbes.

Soil microbes that colonize plant roots and are responsive to differences in plant genotype remain to be ascertained for agronomically important crops. From a very large-scale longitudinal field study of 27 maize inbred lines planted in three fields, with partial replication 5 y later, we identify root-associated microbiota exhibiting reproducible associations with plant genotype. Analysis of 4,866 samples identified 143 operational taxonomic units (OTUs) whose variation in relative abundances across the samples was significantly regulated by plant genotype, and included five of seven core OTUs present in all samples. Plant genetic effects were significant amid the large effects of plant age on the rhizosphere microbiome, regardless of the specific community of each field, and despite microbiome responses to climate events. Seasonal patterns showed that the plant root microbiome is locally seeded, changes with plant growth, and responds to weather events. However, against this background of variation, specific taxa responded to differences in host genotype. If shown to have beneficial functions, microbes may be considered candidate traits for selective breeding.

A growing number of studies report an influence of plant genotype on various facets of the rhizosphere microbiome, but the majority are conducted under greenhouse and laboratory conditions where environmental variability is controlled (1–5), or are limited in scope (3, 6–10). Studies with small sample numbers (e.g., in the hundreds), or that lack replication over space or time, are limited in their power and may yield spurious results. As a result, they can both overestimate the importance of plant genotype and underestimate genotype interactions with the environment (11). Furthermore, although small-scale studies can assess the effect of plant genotype on abstracted measures of the microbiome, such as within (alpha-diversity) and between (beta-diversity) samples, they lack the power to quantify the effect of plant genotype on abundances of specific taxa comprising the microbiome. Thus, large-scale studies with replication are required to identify heritable taxa, i.e., those whose variance in abundance across samples has a significant plant genotype component (12). Heritable microbes of the rhizosphere (the region of soil surrounding, and chemically influenced by, plant roots) should exhibit, across a plant population, differential abundances significantly influenced by variation in plant genotype (13, 14). Heritable components (e.g., taxa or functions) of the soil microbiome remain to be identified for major crop species.

Maize (Zea mays L. subsp. mays) is a globally important crop with over a billion metric tons of production in 2016 (15) and is used for a variety of food and industrial products (16). Maize possesses exceptional phenotypic diversity that is influenced by genotype, which has the potential to influence the rhizosphere microbiome. For instance, maize lines with mutations affecting their carbon storage patterns have been shown to harbor distinct rhizosphere microbiomes; mutant lines with the sugary endosperm su1 gene have been associated with higher fungal and Gram-positive bacterial biomass compared with lines with the shrunken endosperm sh2 gene, and their overall communities differed (7). Maize lines can also differ in their root structures, with effects on the rhizosphere microbiome (17). Heritable microbiota of the maize rhizosphere remain to be identified.

The largest set of public maize germplasm used to dissect genetically complex traits (i.e., lines developed for nested association mapping, so-called NAM lines) (18) consists of ∼5,000 recombinant inbred lines. Analysis of genotype–phenotype relationships using NAM lines has revealed the genetic architecture underlying a variety of complex quantitative traits, such as flowering time, height, stalk strength, and disease resistance (19). If the parents of the NAM lines, selected for maximal genetic diversity while still flowering in summer in North America, show variation in their microbiomes, then subsequent analysis of a far greater number of NAM inbred lines would be warranted. In a previous study of the maize rhizosphere microbiome of NAM parent lines, we analyzed ∼500 samples derived at flowering time from 27 NAM parental lines planted in a randomized block design replicated in five fields (20). This analysis was sufficient to show that plant genotype drove subtle differences in rhizosphere bacterial diversity. However, as with other small-scale analyses, this study did not have sufficient power to identify the specific taxa that responded to differences in maize plant genotype (20).

Here, we performed a longitudinal analysis of the 27 NAM parental line genotypes planted in five fields, yielding >4,800 samples collected weekly over the course of the growing season. In addition, we partially replicated the study 5 y later with a subset of the NAM parent genotypes sampled over the growing season in one field. Microbial diversity was assessed with Illumina sequencing of the V4 region of the 16S rRNA gene for all samples, resulting in an extremely large assessment of rhizosphere microbial community membership. We chose 16S rRNA gene sequencing over metagenome sequencing primarily due to the very large size of the sample set and to specifically identify heritable taxa. These data were analyzed with respect to climate variables, field location, plant age, and genotype. We report well-powered identification of heritable taxa of the maize rhizosphere, over and above the prior beta-diversity differences observed in ref. 20.

Results and Discussion

Maize Rhizosphere Core Microbiome.

The rhizosphere microbiome examined here includes bacteria and archaea inhabiting the combined “rhizosphere” and “rhizoplane” (21) of maize. A total of 4,911 samples and ∼627 million 16S rRNA reads were generated for these analyses (). Across all samples, we detected 55 bacterial phyla, which represent a sizable fraction of the current ∼92 named phyla (22). As expected, maize rhizosphere and bulk soil microbiome differed, but all phyla were detected in both sample types. As microbes that are consistently present across samples likely provide critical ecological functions (23), we searched for a core microbiome, defined here using a stringent criterion: the set of operational taxonomic units (OTUs) common to all rhizosphere microbiome samples. Seven OTUs comprised the core maize microbiome in 100% of samples. All seven were taxonomically assigned to the phylum Proteobacteria, including three Alphaproteobacteria (Agrobacterium, Bradyrhizobiaceae, Devosia), two Betaproteobacteria (both Comamonadaceae), and two Gammaproteobacteria (Pseudomonas and Sinobacteraceae; ). The core microbiome for OTUs shared across 95% of samples included 251 OTUs from a wider diversity of phyla. The seven OTUs present in 100% of 2010 samples were also present in 100% of the 2015 samples. Given that we assessed the diversity of a single plant species, a core microbiome might be expected. However, the observation that seven OTUs were present in all maize rhizosphere samples, in five fields across three states, and across a 5-y timespan, implies that these taxa are highly persistent and ubiquitous in agricultural soil.

Plant Age and Field Are the Main Drivers of Rhizosphere Microbiome Diversity.

To identify factors that shape the differences between maize rhizosphere microbiomes (beta diversity), we applied the phylogenetically aware UniFrac metrics, which provide a measure of dissimilarity between any two microbiomes based on shared membership (unweighted UniFrac), and also take into account the abundances of the lineages (weighted UniFrac). Principal coordinates analysis (PCoA) of unweighted UniFrac distances showed that regardless of field, rhizosphere microbiomes showed a strong patterning according to plant age (Fig. 1). Principal coordinate (PC) 1 explains the majority of the variance in the data, and plant age (sampling week) maps well onto this PC. Samples form a gradient from left to right along PC1, ranging from the bulk soil and early samples to the samples taken at later weeks (Fig. 1). The second PC of the unweighted UniFrac PCoA shows a strong patterning of samples by field (Fig. 1). The age patterning, but not field patterning, was stronger when we used the weighted UniFrac metric (Fig. 1).

Fig. 1.

Environment and age strongly structure the rhizosphere microbial community. PCoA of unweighted (A and B) and weighted (C and D) UniFrac distances for maize rhizosphere and bulk soil samples. A and B show the same projection of the data, as do C and D. Symbols represent microbiomes and are colored by plant age (A and C; sampling week; see color gradient) by environment (B and D; by the specific field of origin). The first three PCs are plotted with the percentage of variation explained by each PC.

Additionally, we analyzed the data for clustering of samples by Adonis testing with unweighted and weighted UniFrac distances, which revealed that plant age is the strongest factor shaping these rhizosphere communities, followed by field, and then plant genetics [R2 of 0.10363, 0.08684, and 0.00770 (unweighted UniFrac); 0.38644, 0.07716, and 0.00724 (weighted UniFrac), for age, field, and inbred maize line, respectively, P value < 0.001 for all tests]. To further assess the relative importance of field, plant age, plant genotype, and their interactions on these patterns, we performed variance decomposition for each principal coordinate of the PCoAs (). This confirmed that the majority of variance exhibited by PC 1 in the unweighted UniFrac PCoA (Fig. 1) is attributable to plant age, and variance in PC2 was attributable to location. For PC1 of the unweighted UniFrac analysis, the second largest component of the variance is the interaction of plant age, genotype, and location. Together, these results imply significant gene-by-environment interaction, implying that plant genotype impacts the community composition differently depending on the plant’s age and location.

Plant Genotype Distance Matrix Versus Microbiome Distance Matrices.

We used a distance matrix for the maize genotypes (24) to ask if the overall genetic dissimilarity of any two maize lines was predictive of the dissimilarity of their rhizosphere microbiomes. Correlations between the maize genotype distance matrix and the UniFrac distance matrices yielded weak (R2 = 0.04 and 0.03) effects with opposite (inverse and positive relationships) directions for unweighted and weighted UniFrac, respectively (). This implies that the overall genetic differences between the maize lines mostly fail to predict the overall diversity differences in their rhizosphere microbiomes (i.e., the most genetically distant NAM lines of maize do not have the largest UniFrac distances between microbiomes). Thus, any plant genetic effects on the diversity of the rhizosphere microbiome will likely manifest at the level of specific taxa.

Identification of Taxa with Differential Abundances Across Plant Genotypes.

To identify microbial taxa responsive to plant genotype, we first compared the abundances of OTUs in pairwise categories of maize groups [six broad functional groups: sweet-corn, popcorn, tropical, stiff-stalk, non-stiff-stalk, and mixed (25)] using a mixed-model approach with Bonferroni correction (). Pairwise comparisons between broad maize groups resulted in 83 instances of differentially abundant OTUs. We observed that these instances involved 48 OTUs, with the same OTUs differentially abundant in several pairwise comparisons ().

Compared with other maize groups, the sweet-corn group had a strong tendency to enrich for taxa that have been associated with nitrogen fixing activities (e.g., Bradyrhizobiaceae, Burkholderia, Rhizobium, Sphingomonas, and Oxalobacteraceae). The sweet-corn inbred lines harbor mutations that cause higher sucrose and glucose concentrations and lower starch production in the endosperm (26), which may lead to different root exudate availability to the microbiome.

We performed similar paired tests of the 27 inbred maize line (), and observed 255 significant instances involving 92 OTUs. Again, maize lines with characteristics likely impacting exudate profiles shaped the rhizosphere microbiome differentially. For instance, members of the Chloracidobacteria and Pedosphaerales orders, and the Opitutaceae family, were enriched on the non-stiff-stalk Oh43 maize line, which has a defective invertase for sucrose uptake in the root system (27). A Paenibacillus sp. was enriched in the tropical CML247 line. This genus includes free-living N2 fixers (28). Many of the OTUs differentially abundant between maize lines belong to the same bacterial families, and they might be expected to show similar patterns of enrichment or depletion. Indeed, we noted that for the 20 cases where OTUs had the same taxonomic classification, all responded in the same manner on the same maize lines. This may indicate that the related OTUs perform similar functions. In the 2015 dataset, no OTUs had significantly different abundance by inbred maize line or maize group after multiple comparison corrections. This is likely due to the reduced power of the smaller dataset.

Microbiome Richness Across Plant Genotypes.

Microbiome richness (alpha diversity) varied between broad maize groups and between inbred lines. Specifically, “mixed” maize lineage had lower diversity than other groups, including the tropical and non-stiff-stalk maize. In accordance, the mixed inbred line Mo18W had lower microbiome richness compared with many other lineages (; this line was not planted in 2015). Differences in richness replicated 5 y later; compared with other maize lines, the non-stiff-stalk Mo17 had higher alpha diversity, and the sweet-corn IL14H had lower diversity in both years ().

Identification of Heritable OTUs in the 2010 Field Study.

To estimate the amount of variance in OTU abundances across all samples that could be attributed to plant genetics, while controlling for environmental factors, we calculated the broad-sense heritability (H2) for OTUs shared by ≥80% of samples (n = 792). We identified 143 OTUs with significantly more heritability than expected by chance based on 5,000 random permutations of the data (empirical P value ≤ 0.001; Fig. 2 and ).

Fig. 2.

Broad-sense heritability of individual OTUs for the 2010 field study. The broad-sense heritability (H2) is shown for the 100 OTUs with highest H2 values. Circles show the actual H2 values for each OTU in decreasing order and blue distributions show the corresponding H2 values from 5,000 permutations of the data. Red circles indicate OTUs with P values ≤0.001. Taxonomies shown are most specific for each OTU.

Overall, the broad-sense heritabilities for these OTUs were relatively low, with a range of 0.15–0.25 (on a scale from 0 to 1, where 1 means that all variance is attributed to genetics). This indicates that environmental factors contribute most of the variance observed in the relative proportion of OTUs across samples, as expected given the strong environmental effects discussed earlier. We observed partial overlap between the heritable OTUs and the 92 OTUs that showed significant differences in abundances in the pairwise comparisons (inbred maize lines and broad maize groups; ).

Mapping the heritable OTUs onto a common phylogeny revealed some clusters of related taxa (Fig. 3). The heritable OTUs are diverse, including 26 Alphaproteobacteria, 9 Betaproteobacteria, 12 Actinobacteria, 6 Verrucomicrobia, and 8 Bacteroidetes. Others belong to WS3, Beta- and Gammaproteobacteria, Planctomycetes, Firmicutes, Chloroflexi, Acidobacteria, Gemmatimonadetes, and interestingly, the Archaeal phylum Crenarchaeota (candidatus Nitrososphaera). Five of the seven core OTUs described previously (Agrobacterium, Devosia, both Comamonadaceae, and the Sinobacteraceae OTUs) were among the heritable taxa. Although the Alpha, Beta, Delta, and Gammaproteobacteria were well represented among the heritable OTUs, the genus Pseudomonas was not (except one OTU, classified as Pseudomonas viridiflava). Of note, only eight OTUs did not match the Greengenes reference database, indicating that a majority of heritable OTUs have previously been observed, at least by sequencing. Indeed, the OTUs with the greatest heritability values (>0.2) could be classified at least to the family level, and many to genus and species.

Fig. 3.

Phylogeny of heritable OTUs. The tree shows the phylogenetic relationships of OTUs found in 80% of all rhizosphere samples. Heritable OTUs have branches colored in red. Major phyla are color-coded in the outer circle, and the deepest named taxonomy is given for each OTU. The tree is a subset of the Aug 2013 Greengenes tree. (Scale bar indicates inferred mutation rate per site.)

The heritability of individual rhizosphere OTUs was lower than for most traditional agronomic traits (0.15–0.25, whereas the heritability of plant yield is around 0.3 and flowering time can be 0.9 or higher). Nonetheless, there is a noticeable effect of plant genetics on the relative abundances of specific taxa. Since plant genotype selects on microbial phenotypes, heritable taxa likely encode functions that are phylogenetically restricted, with function manifested at the taxon level.

Heritable OTUs in the 2015 Field Study.

Analysis of the 2015 dataset identified five OTUs that were significantly heritable (P < 0.001; ). One OTU, an Ellin517 taxon, was overlapping between the 2 y, and another was taxonomically related to OTUs heritable in 2010 (Erythrobacteraceae). The OTUs that did not overlap between years were detected in both years; however, they were present in a smaller number of samples (55–70%) than the Ellin517 OTU in 2010 (>95%). These zero counts limit variance in abundance, and hence heritability. Note that the heritability estimates were higher in 2015 than in 2010, likely because of the small 2015 sample size. Moreover, the 2015 data stem from just one field, so that environmental variability is reduced. Indeed, higher H2 estimates emerged when the 2010 data were subsetted to equivalent sample size (). These findings underscore the importance of large datasets and multiple environments in heritability analyses.

Pseudomonas Bloom.

The age gradient allowed a weekly view of rhizosphere microbiome development over the growing season. Averaging the abundances of OTUs across fields and maize inbred lines, we observed that bacteria belonging to the Pseudomonas genus bloomed from week 8 onward (average increased from ∼2.7% during weeks 1–7 to ∼44.5% thereafter; ). This bloom was apparent across all fields. Pseudomonas was not, however, dominant in bulk soil sampled across the season, indicating the bloom was plant-driven (). Of note, Pseudomonas has previously been shown to dominate the maize rhizosphere microbiome (29–31).

After week 8, 43.1% of the total microbial community belonged to just three Pseudomonas OTUs. The top 10 most abundant unique sequences within these OTUs accounted for 78.6% of the 2015 and 73.1% of the 2010 Pseudomonas sequence data. Bray–Curtis dissimilarity clustering revealed that the Pseudomonas community structure differed by field (), and that this field-specific structure remained constant over time. Thus, each field harbored a distinct Pseudomonas population structure that was stable over the growing season.

The three Pseudomonas 97% OTUs that dominated in the latter half of the 2010 sampling season were also present in the 2015 dataset, albeit at lower abundance. Moreover, the three Pseudomonas 100% identity OTUs dominant in 2010 were also the most abundant in 2015, although at lower relative abundance to all other taxa. In 2015, Pseudomonas were slightly more abundant in rhizosphere compared with bulk soil samples (2.0% versus 0.3%; ), but we did not observe the same seasonal bloom. Intriguingly, when we correlated the relative abundances of all unique 100% identity Pseudomonas OTUs between the 2010 and 2015 datasets, we observed a strong relationship (R2 of 0.91 for the 10 most abundant OTUs and R2 of 0.94 for all data), even though the 2010 data included three geographically distinct fields and the 2015 data only one. Thus, the most abundant 100% identity Pseudomonas OTUs in 2010 were also the most abundant in 2015. These data suggest that the field-specific Pseudomonas community structure persists across long time periods, even though the relative abundance of the Pseudomonas as a broader category can vary substantially.

Impact of Climate, 2010 and 2015.

A subset of the rhizosphere microbiome is responsive to weather: their relative abundances were associated with climate variables. We modeled the effects of weather, as recorded from the National Centers for Environmental Information, on the rhizosphere microbiome across the 2010 and 2015 season (). Families of Bacteria and Archaea tended to respond in the same positive or negative pattern for given weather conditions. For instance, Nocardioidaceae, Caulobacteraceae, and Sphingobacteriaceae responded positively (across more than 80% of OTUs) to temperature, whereas members of the Ellin515 and Opitutaceae, Solibacteraceae, and Oxalobacteraceae responded negatively. For same-day precipitation, Sphingobacteriaceae and Cytophagaceae families responded positively 100% of the time, while the Verrucomicrobia Ellin517 and Ellin515, Hyphomicrobiaceae, Chthoniobacteraceae, and Chitinophagaceae responded negatively to precipitation at least 80% of the time. When cumulative precipitation for the prior 2 or 3 d was tested, the OTUs that responded in a positive manner changed. Chitinophagaceae, Hyphomicrobiaceae, Sphingomonadaceae, Comamonadaceae, Gaiellaceae, Nocardioidaceae, Microbacteriaceae, and Opitutaceae responded positively to long-term precipitation for both 2- and 3-d sums. Members of the Pseudomonadaceae responded negatively to long-term precipitation, indicating that wetter and colder conditions disfavor this group. The weather-responsive subset of the rhizosphere microbiome adds variability to the age-gradient trends and underscores the importance of multiple time measurements in field studies.

Similar patterns of responsiveness to weather were replicated 5 y later (). Since there is limited overlap between the specific 97% OTUs detected in 2010 and 2015, and because we had observed that OTUs belonging to the same bacterial family usually had similar responses to weather events in the 2010 data, we compared the 2010 and 2015 weather responses by examining if OTUs of the same family responded to weather similarly in both years. Of note, members of the family Hyphomicrobiaceae responded negatively to temperature in both years. In addition, several taxa showed a positive response to precipitation in both years (2- and 3-d precipitation sums), including Saprospiraceae, Hyphomicrobiaceae, Verrucomicrobia, A4b (Anaerolineae), Ellin6075 (Chloracidobacteria), Sphingomonadaceae, and Chitinophagaceae.

Interestingly, while many of the heritable OTUs belonged to the same families, these also included taxa responsive to weather events (e.g., members of the Oxalobacteraceae, Nocardioidaceae, Sphingomonadaceae, and Chitinophagaceae). The observation that the heritable set of OTUs and the weather-responsive set are related suggests that many of the heritable taxa are r-selected: they likely react to short-term changes in carbon availability that are modulated by the plant (32).

Prospectus.

This well-powered study revealed the heritable components of the maize rhizosphere under field conditions. Despite strong environmental patterning, we identified close to 150 OTUs with significant heritability. These were highly diverse, including Archaea. Many are related to taxa with putative beneficial functions. While our 2015 replication study was limited in scope, it nevertheless partially recovered the 2010 findings concerning heritability and weather responses. The maize inbred lines studied here were selected because they are the parents of the recombinant inbred lines (RILs) that comprise the largest set of public maize germplasm used to dissect genetically complex traits (18). The present work was conducted in part to ascertain if the NAM population might exhibit enough variance in the rhizosphere microbiome to warrant an expanded study using the NAM RILs, aimed at mapping plant genes that drive aspects of the rhizosphere microbiome. Our results suggest that there is indeed a variable selection of NAM RILs upon specific microbial abundances.

A next step is to discern the plant genes that associate with these taxa and to better characterize their functions. GWAS in mammals highlight an effect of immune-related genes on the microbiome, and similarly, studies in Arabidopsis point to effects of plant immune-related genes on microbiome composition (4). It remains to be seen if immune-related phenotypes (33) are as important in shaping the maize rhizosphere microbiome.

Future studies may address other aspects of the microbiome that are more directly related to function, using approaches such as metagenomics and metabolomics. The rhizosphere microbiome phenotypes could also be defined closer to the root surface (e.g., rhizoplane or endosphere), but then the impact of the plant on the soil may be lost. Another possibility is to incorporate the rhizosphere microbiome as a factor in determining desirable traits in maize, such as yield or drought tolerance, and to assess how plant genotype interacts with the rhizosphere microbiome to shape host phenotype.

Materials and Methods

In 2010, we sampled replicate NAM parental lines in three fields in New York State (USA) weekly, and two more fields (Illinois and Missouri) were sampled once at flowering time. Combined with bulk soil from each field, this yielded 4,866 samples. From these samples, we generated ∼627,000,000 sequences of the V4 region of the 16S rRNA gene using Illumina MiSeq sequencing. We later partially replicated the field sampling in 2015 in one location (New York) with a subset of the maize lines and sampling times. From the replicate field study, we obtained 45 rhizosphere samples, which we processed similarly to the 2010 set. Beta diversity, Adonis testing, core microbiome, and raw alpha-diversity values were calculated with the QIIME (34) software package. Differential abundance, alpha-diversity significance tests, microbial response to weather, and heritability results were generated using the lme4 (35) mixed-model package. More specific details for sampling handling and data analyses are described in .