Association analyses of host genetics, root-colonizing microbes, and plant phenotypes under different nitrogen conditions in maize

Curation statements for this article:
  • Curated by eLife

    eLife logo

    Evaluation Summary:

    The authors have produced a rich dataset that illuminates how root-associated bacteria differ among diverse maize lines, under low and high nitrogen treatment. The authors make use of this dataset to describe interesting patterns related to the genetic architecture of factors associated with maize rhizosphere microbiome assembly, although it is not clear yet whether the microbiome is an agent of natural selection in this case, or whether other selective forces shape maize roots in a manner that in turn affect smicrobiome recruitment. In addition to interesting insights reported in the present manuscript, these data are likely to be used for and/or compared to, in many future studies.

    (This preprint has been reviewed by eLife. We include the public reviews from the reviewers here; the authors also receive private feedback with suggested changes to the manuscript. Reviewer #1 agreed to share their name with the authors.)

This article has been Reviewed by the following groups

Read the full article See related articles

Abstract

The root-associated microbiome (rhizobiome) affects plant health, stress tolerance, and nutrient use efficiency. However, it remains unclear to what extent the composition of the rhizobiome is governed by intraspecific variation in host plant genetics in the field and the degree to which host plant selection can reshape the composition of the rhizobiome. Here, we quantify the rhizosphere microbial communities associated with a replicated diversity panel of 230 maize ( Zea mays L .) genotypes grown in agronomically relevant conditions under high N (+N) and low N (-N) treatments. We analyze the maize rhizobiome in terms of 150 abundant and consistently reproducible microbial groups and we show that the abundance of many root-associated microbes is explainable by natural genetic variation in the host plant, with a greater proportion of microbial variance attributable to plant genetic variation in -N conditions. Population genetic approaches identify signatures of purifying selection in the maize genome associated with the abundance of several groups of microbes in the maize rhizobiome. Genome-wide association study was conducted using the abundance of microbial groups as rhizobiome traits, and n=622 plant loci were identified that are linked to the abundance of n=104 microbial groups in the maize rhizosphere. In 62/104 cases, which is more than expected by chance, the abundance of these same microbial groups was correlated with variation in plant vigor indicators derived from high throughput phenotyping of the same field experiment. We provide comprehensive datasets about the three-way interaction of host genetics, microbe abundance, and plant performance under two N treatments to facilitate targeted experiments toward harnessing the full potential of root-associated microbial symbionts in maize production.

Article activity feed

  1. Author Response

    Reviewer #1 (Public Review):

    This manuscript describes a large field experiment in which diverse maize lines were grown under high and low nitrogen conditions. The authors surveyed the root microbiome and performed several different analyses in search of evidence of host control over the microbiome. These analyses include clustering the 16S amplicon sequencing data into 'microbial traits', searching for evidence of heritability and selection, GWAS and transcriptomics. The final section of the paper looks for patterns across the different analyses and highlights a single microbial trait, a corresponding GWAS hit and evidence that one allele correlates with both microbial abundance and plant health.

    The strengths of this manuscript include an impressive dataset, nicely compiled figures and an impressive compilation of distinct analyses to reduce complexity and reveal interactions. This manuscript is interesting on its own and will also serve as a template for future analyses on similar or distinct datasets. In several places, this manuscript could be improved by more clearly articulating the logic behind choices made for the analyses (see specific comments below).

    We thank the reviewer for the comments. We clarified the logic for the analyses as shown below.

    1). How the 'microbe traits' are defined is critical for interpreting the rest of the paper and so needs to be explained more clearly in this manuscript. From the text, I interpret that ASVs were first clustered into genera and then further clustered based on differential abundance between high and low N in both years. If so, this would seem to exclude a bunch of microbes that were specific to the 202 genotypes planted only in the second year. However, Supp Fig 2 shows a handful of 'microbe traits' that are not differentially abundant (grey). Please clarify. Also, the authors should clarify how the 'names' for each microbe trait defined.

    In the revised manuscript, we added a section “Clustering of ASVs into microbial groups” in the materials and methods. Additionally, we added a supplementary figure (or Figure R1) to better explain how microbial groups were clustered. We also added the original plots used to define the 150 microbial groups in Supplementary File 6. See below for detailed responses.

    Added Text:

    Clustering of ASVs into microbial groups

    ASVs were clustered into groups of rhizosphere microbes at the family, genus, and species level using a procedure described previously (Meier et al., 2021). First, the 3,626 ASVs in the present study were grouped at the family levels (the lowest taxonomic rank for which all ASVs were successfully annotated) and the phylogenetic tree derived from 16S V4 alignment was plotted alongside taxonomic annotation at the genus and species level. Because the ASVs are derived from short reads and may constitute variations not covered in the SILVA database, annotation at the genus and species level was often not possible. To close these gaps and form biologically meaningful groups of ASVs at low taxonomic ranks with better confidence, we examined the overall abundance of each ASV as well as the differential abundance in response to the N treatment alongside the sequence-based clustering. The premise here is that ASVs derived from biologically closely related individual microbes are similarly abundant in our dataset and respond similarly to the N treatment imposed on the field, in addition to similar 16S sequences due to common ancestry. An example is given in Supplementary Figure 8 with a subset of ASVs assigned to the Burkholderiaceae family. The plots used to determine all 150 microbial groups in this study are available in Supplementary File 6.

    1. Throughout the manuscript, I expected more consideration of which microbe traits overlap between the different analyses. This comes in only at the final figure and so we never get a sense of how the different analyses overlap. For example, 150 microbe traits are defined and assessed for heritability and evidence of selection. I was not clear on whether the traits that were found to be heritable overlapped with those under selection (as one might expect). Lines 387-395 and Supp Fig 5 attempt to synthesize the different analyses but should be expanded to help the reader understand overlap between the analyses.

    We present evidence of links between microbe abundance, plant genetics, and plant performance through heritability analysis, estimation of selection, GWAS, and correlation of microbe abundance with the canopy coverage phenotype. Summarizing the overlap of these analyses of interest was challenging to do concisely for 150 microbial groups and two distinct N treatments, thus we chose to summarize the results of several analyses for the 62 microbial groups that showed either a positive or negative correlation with canopy coverage.

    Supplementary Figure 5 was revised to include the results of the estimation of selection. Furthermore, the section “Heritable and adaptively selected rhizobiota are associated with plant phenotypes” in the results was revised to better explain overlaps between the separate assays.

    Notably, not all traits that are heritable are expected to be under selection, as traits can be heritable, i.e., transmitted from one generation to the next, without impacting the fitness or performance of offspring individuals under the conditions under which recent natural and/or artificial selection has occurred. This point was added to the discussion.

    Lastly, the complete data with the observed values for all 150 microbial groups and all analyses are made available in Supplementary File 3.

    1. I gather that the authors performed the GWAS separately for each 'microbe trait' and nitrogen condition, but then searched for 'hot spots' where the data from different microbe traits was considered in a pooled manner. The logic behind this decision is not clear to me. Why would we expect different microbe traits to be co-localized in the genome?

    Our logic is that plant morphological (i.e., root architecture) or physiological (root exudation) changes may affect several rhizobiome traits, such that a GWAS signal controlling plant metabolism has an effect on several groups of microbes; therefore, it can be detected as a hotspot. In the revised manuscript we now explicitly describe our reasoning for the expectation of microbial trait hot spots in the section “Genes underlying microbe-associated plant loci are preferentially expressed in root tissue”.

    Line 346-349 indicates that no plant loci were found to associate with microbe 'traits' under both nitrogen treatments and speculate as to why (dynamic interactions or not enough statistical power). The traits were defined based on robust differential abundance between nitrogen treatments and if I understand correctly, the GWAS was run separately for each trait and nitrogen treatment, so it seems logical that this method would only yield microbe trait associations that are differential between nitrogen treatments. If I did understand this correctly, I recommend emphasizing this point as it seems to indicate that the methods are working as expected.

    We are grateful for this suggestion. These points were elaborated in the revised manuscript in the section “Genes underlying microbe-associated plant loci are preferentially expressed in root tissue”.

    To investigate this, we performed GWAS using each of the 150 rhizobiome traits. This analysis was done separately for the -N and +N conditions, as N deficiency induces dramatic changes in plant metabolism, including changes in root gene expression (Choudhury et al, 2022) and root exudation (Zhu et al, 2016), and because N applied to the field directly impacts soil and rhizosphere microbiomes (Meier et. al, 2021). However, it is important to emphasize that microbes which did not exhibit differential abundance in response to nitrogen were indeed included in our analysis. It is clear the explanation for our methodology for grouping ASVs into functionally distinct clades/traits was a significant weakness in the previous version of our manuscript and caused significant confusion. We have revised the relevant sections substantially (see above).

    1. Having access to expression data for 298 genotypes is amazing. It would seem logical to try to more directly connect the MAPLs and microbe traits with this expression data. Do the lines that show association with the microbes also show higher expression of the corresponding gene?

    We agree with the reviewer that it would be logical to connect the MAPLs, gene expression, and microbe traits altogether. Modeling three distinct types of variables in a high dimensional genomics setting, however, is non-trivial. In fact, to address the challenge, another graduate student in our group has developed a genome-wide mediation analysis method with the ultimate goal of establishing a causal chain from genotype to intermediate molecular traits (i.e., gene expression or microbial traits) to plant phenotype. In a recent publication (Z. Yang et al., 2022, Genetics), our results from both simulation and empirical analyses suggested that this model could identify mediating genes with certain power, albeit with some limitations. We are working on improving the model and applying it to the current dataset by considering microbial groups as intermediate traits (mediators).

    We feel conducting and including the results of microbiome mediation analysis will exceed the scope of the current study. In addition, it may not be reliable as the gene expression data published by Kremling et al. does not cover two distinct N treatments. The dataset associating MAPLs with particular genes presented here (Supplementary File 5) was provided to be used alongside new gene expression data to investigate particular associations in greater detail in more targeted experiments.

    The authors generated additional RNAseq data from 2 week old plants from 4 genotypes but the logic for the selection of these lines is missing and I am not sure about the relevance of this since the samples were collected from young plants. Is a nitrogen treatment effect observable at 2 weeks?

    As a small validation experiment, we selected 4 diverse and well characterized maize genotypes from the diversity panel to be grown under controlled greenhouse conditions. N treatments were included for consistency with the field experiment and because application of N fertilizer may have a direct influence on rhizosphere microbiomes independent from the host plant. The relevant section in the manuscript was revised for clarity as below.

    To complement the gene expression data provided by Kremling et. al, we selected 4 diverse and well characterized maize genotypes (K55, W153R, B73, and SD40). Plants were grown in a controlled greenhouse environment under standard N and N deficient conditions and gene expression was analyzed in roots and shoots of two-week old seedlings (for details refer to Xu et al, 2022). In agreement with the dataset provided by Kremling et al, significantly higher expression of 97 MAPL genes was observed in root but not leaf tissue compared to (n = 44,049) other genes available in this dataset (Figure 3C). No strong physiological response to N deficiency was expected for 2-week-old seedlings and no significant differences were observed in the pattern of MAPL gene expression between the two N treatments.

    The authors conclude that the gene expression data is consistent with host control over root microbiome (line 371-373) but, as is, I'm not convinced that this analysis supports that statement. Fig 3C is striking on its own, but based on panel B, I suspect that a similar pattern would be observed for 'third leaf' and 'germinating shoot' so it is harder to make a direct connection with the microbe traits.

    We thank the reviewer for the comment. In Kremling’s experiment (Fig. 3B), leaf tissue samples from adult plants are taken both from the tip of the third leaf and from the base, and similar to roots we do observe higher expression of MAPL genes in the leaf base of adult plants, although other parts of the leaf show the opposite trend. We agree that the current dataset is insufficient to explain this observation, and a direct link between microbiome features and root gene expression can not be conclusively established at the moment. We revised the wording in the relevant section:

    Revised Text: 371-373: Collectively, these data are consistent with the hypothesis that root-associated microbial communities are at least in part genetically controlled by the host plant in a process mediated by plant gene expression.

    1. The authors report that 62 microbe traits associated with canopy coverage, a very exciting result! However, again, this confuses me based on how the microbe traits were defined. To be considered a microbe trait, the microbes had to show differential abundance across the treatments. The logic for how this could manifest in phenotypic changes in both treatments needs to be elaborated.

    How the microbial groups were defined in this study was a major point of confusion, and we included a more thorough explanation of the procedure above. Differences in the pattern of response to N treatments whether positive, negative, or no response), as well as differences in overall abundance were used to separate sister clades of ASVs; however clades of ASVs which did not show differential responses to N treatment were indeed included in our analyses. The purpose of this section was to associate the abundance of microbial groups under either N treatment (not to be confused with the differential abundance between N treatments) with plant performance in the field. More targeted experiments are required to determine the direction of causation and a potential mechanism by which microbe abundance could influence phenotypic changes in the host. This point was added in the discussion.

    1. The final figures summarizes correlations for one microbe trait across the different analyses and looks very promising, especially for noisy field data. The authors are careful to not overstate this finding, perhaps a bit too conservative. They see a significant correlation between microbe abundance and canopy coverage that also correlates with allele frequency. The difference in canopy coverage by allele frequency is not significant, but shows a similar trend and this is not necessarily surprising given all the other factors that will influence this one trait. I expected a comment on gene expression of the genes in the locus and perhaps a peek at the other plant traits to see if any of them also show a similar trend.

    We thank the reviewer for this suggestion. As indicated above, we are hesitant to make strong statements using the gene expression data published by Kremling et. al. because experimental conditions differed from our study.

    As the reader may indeed wonder about the gene expression of the genes in the exemplary locus in Figure 5, we supplied gene expression information of the three relevant genes in an additional supplementary figure (or Supplementary Figure 10).

    We do see elevated gene expression in roots for two out of the three genes, which matches the previously observed trend. A brief literature review of the same genes indicates a possible link to root hair physiology and an altered microbiome. We revised the text as below and incorporated the speculation into the discussion.

    Reviewer #3 (Public Review):

    Strengths:

    The choice of 230 genotypes from a well-known maize diversity panel, with accompanying SNP genotype data, was a good one for this purpose, given the focus on selection during intensive breeding during the 20th century in heavily fertilized conditions.

    The very large dataset (N>3000 with replication of 230 genotypes) is a useful source of information on maize rhizosphere bacterial microbiomes, and the availability of the host genotype SNP data is an especially useful and unusual feature. The authors used a relatively newly developed (2018) Bayesian computational approach to characterize genetic architecture of rhizosphere composition, which is an interesting advance in the microbiome field. The same tool makes inferences about whether each SNP underlying rhizosphere features shows signatures of past selection (inferred as the variation in effect size relative to the minor allele frequency).

    Weaknesses:

    The BayesS results classifying rhizosphere-related SNPs as under positive, negative, or no selection appear to be over-interpreted. First, it is not clear that this method is meant for comparing current patterns of selection in contrasting environments (as in this N+/N- experiment), but rather for detecting signatures of selection in the distant past. Second, it IS clear that this method only reveals signatures of selection on a locus or SNP, and cannot confirm that selection is acting on a particular trait. The Zeng et al. 2018 paper states this quite clearly. The authors of this manuscript did not attempt to rule out that the loci classified as under selection do not have pleiotropic effects on (or are linked to) traits other than rhizosphere microbes. Occam's razor in this case suggests that these loci control root traits that are important for plant survival and also happen to affect microbiome composition. No functional benefit of these microbes was demonstrated beyond correlations with plant phenotypes.

    Thank you for the insightful comments. We agree with the reviewer that the signatures of selection reflect the selection in the distant past during the plant evolution but not the recent crop improvement processes. And that we can’t rule out the possibility that the selection might be acting on the plant fitness and in turn affect microbiome composition.

    Per the reviewer’s comments, we estimated the selection gredients (see detailed responses below) and modified the results of the selection section substantially.

    Additionally, we have clarified the concerns raised by the reviewer in the discussion section as below.

    The BayesS method leverages the relationship between the variance of SNP effects and MAF as a proxy of natural selection in the distant past. This method detects signatures of natural selection on SNPs associated with microbiome traits but is not indicative of selection acting on the particular microbes.

    To further approve the beneficial effects of the microbes on the plant fitness, additional functional analyses (i.e., inoculation experiments) are warranted, and that naturally occurring microbe-plant symbiosis may be harnessed for further crop improvement.

    It is unclear how substantially different the N+/N- treatments were from each other. The entire experiment followed commercial corn, so presumably all plots had been fertilized within the past year. The soil chemical profiles were not subsequently tested, and basic analyses (such as comparison of plant growth in the two treatments) are missing.

    This is a valid concern. The +N plots were applied urea (dry fertilizer) as a source of N at the rate of 120 lbs/acre (approximately 134.5 kg/ha), and the -N received no treatment with the assumption that N had been exhausted by commercial corn planted in the previous year. In addition, we now cite the results of UAV imagery data analysis conducted using the same field experiment (Rodene et al. 2022) comparing plant growth in the two treatments and indicated the different N reponses between +N and -N treatments. We have clarified this in the revised manuscript.

    240-244: Paired-end 16S sequencing of 3,313 rhizosphere samples from 230 replicated genotypes of the maize diversity panel (Flint-Garcia et al., 2005) were collected from field experiments conducted under both +N and -N conditions (Materials and Methods, Supplementary Figure 1). At the time of sampling, visible phenotypic differences were observable between +N and -N plots as measured through aerial imaging (Rodene et. al, 2022). Sequencing produced 216,681,749 raw sequence reads representing 496,738 unique amplicon sequence variants (ASVs) (Materials and Methods).

    Furthermore, the assumption that microbiome differences between fertilization treatments are driven by some activity of the plant host (lines 94-99) is not justified - direct responses of the microbes to N addition would almost certainly be reflected in the rhizosphere, since rhizosphere microbiomes are almost entirely derived from the surrounding soil. No bulk soil samples were collected as controls to rule this out. There is clear evidence that N fertilization directly affects microbial communities in the soil. However, the direct effect of N fertilization as opposed to changes induced indirectly via altered

    This is also a valid concern, and we thank the reviewer for the insight. The issue of a direct effect of N fertilization on microbial communities as opposed to an indirect effect via the plant host was partially addressed in a previous study (Meier et.al, 2021). The focus here was to find candidate associations between plant genetics, microbial groups, and plant performance under two N treatments, and to lay the basis for more targeted experiments in which causality can be inferred.

    We edited the text to clarify this:

    94-99: It was observed previously (Zhu et al., 2016) that soil microbial communities drastically change in response to N fertilization. In bulk soil, this is likely due to a direct effect of N application or lack thereof. In rhizospheres, however, only a subset of the observed changes can be attributed to direct effects of N fertilization, while particular microbial groups may be subject to indirect effects induced by the plant host in response to N availability or deficiency (Meier et al., 2021). A possible explanation for this could be that the vast majority of the interval between maize domestication and the present, beneficial plant-microbe interactions have evolved in low-input agricultural systems characterized by relative scarcity of nutrients, predominantly nitrogen (N) (Brisson et al., 2019).

    The authors report how various patterns differ between the N+ and N- treatments- for example, more rhizobiome features had nonzero heritability in N- than in N+. However there is no statistical support for this apparent difference, i.e. no direct test of heritabilities in the two treatments. Nor does they test for possible differences (or lack thereof) in the magnitude of heritability between treatments. This incomplete style of analysis was repeated several times in the paper, e.g. for comparing patterns of selection between treatments, and for comparing correlations between rhizobiome features and plant traits between features.

    We agree that appropriate statistical tests can strengthen our results. See detailed response below (Reviewer #3 Recommendations for the authors).

    The methods sections contain inconsistencies and omissions that made it difficult to evaluate some of the claims. For instance, lines 147-148 describe collection of rhizosphere samples from entire root crowns, but the appendix (lines 774-778) describe collection of rhizosphere from roots that fit in 50 mL tubes. So it is unclear which part of the root crown was actually used, and whether the focal root type was consistent for all samples. Similarly, the appendix states that the B73xMo17 check genotype was used to correct for small-scale geographic differences (747-748), but no additional detail is provided nor are the results of this process reported. In general, the descriptions of statistical analyses lack important details. For example, by definition a constrained ordination (CAP) analysis requires a formula to be specified, but this was not reported in the paper, making it impossible to interpret the meaning of the constrained axes shown in the figures. Ordinations also require the use of a distance or dissimilarity metric, the choice of which affects interpretation - the metric used in this paper was not provided.

    We thank the reviewer for these comments. The relevant sections in the manuscript were revised below to clarify the sampling procedure and provide missing information about the check plants and the CAP analysis.

    146-148: Eight weeks after planting (2018: July 10 and 11; 2019: July 30, 31 and August 1), plant roots were dug up to a depth of 30 cm and rootstocks were manually shaken to remove and discard loosely adherent bulk soil. For each plant, all roots thus exposed were cut into 5 cm pieces and homogenized, and 20-30 ml randomly selected root material (with adherent rhizosphere soil) was collected to generate the rhizosphere samples (Supplementary Methods).

    159-163: Raw ASV reads were subjected to a series of filters to produce a final ASV table with biologically relevant and reproducible 16S sequences (Supplementary File 1). For the constrained ordination (CAP) analysis performed here, the weighted Unifrac distance metric was used with model “distance ~ year + genotype + nitrogen + block + sp + spb”. Only ASVs that were highly abundant and repeatedly observed in both years of sampling were considered for downstream analysis.

    746-749: In each of 12 split plot blocks per quadrant, at least one subplot was randomly selected and assigned the hybrid genotype (B73xMo17) to be used as a check to test for differences between geographical field locations. Two check genotypes (B73xMo17 and B37xMo17) were used in 2018, and a single check genotype (B73xMo17) was used in 2019. Plant growth across the field was determined uniform within quadrants using the checks as reported in a sister study on the same experimental field (Rodene et. al, 2022).

    774-778: To wash the tightly adherent rhizosphere soil layer off the roots, tubes were filled up to the 40 ml mark with autoclaved PBS buffer (46 mM NaH2PO4, 60 mM Na2HPO4, 0.02% Silwet-77), and shaken horizontally at 8000 rpm for 30s. Rhizosphere suspension was filtered through a 100 μm nylon cell strainer (Celltreat Scientific Products, Pepperell, MA, USA) into a fresh 50 ml tube to capture root debris and large soil particles.

    Finally, many of the analyses throughout the paper take the form of testing 150 different rhizobiome traits, one by one, and then reporting the number of significant results (e.g., differential abundance between N+/N-, significant heritability, selection, correlations with plant traits). This suggests a potentially severe risk of false positives due to repeated multiple testing. After the p-values are corrected for the very large number of statistical tests (using Bonferroni, FDR, or similar) many of the conclusions might change.

    We agree that the large number of tests may lead to false discoveries. To address this, a multiple testing correction method was applied to increase the stringency of the GWAS analysis.

  2. Evaluation Summary:

    The authors have produced a rich dataset that illuminates how root-associated bacteria differ among diverse maize lines, under low and high nitrogen treatment. The authors make use of this dataset to describe interesting patterns related to the genetic architecture of factors associated with maize rhizosphere microbiome assembly, although it is not clear yet whether the microbiome is an agent of natural selection in this case, or whether other selective forces shape maize roots in a manner that in turn affect smicrobiome recruitment. In addition to interesting insights reported in the present manuscript, these data are likely to be used for and/or compared to, in many future studies.

    (This preprint has been reviewed by eLife. We include the public reviews from the reviewers here; the authors also receive private feedback with suggested changes to the manuscript. Reviewer #1 agreed to share their name with the authors.)

  3. Reviewer #1 (Public Review):

    This manuscript describes a large field experiment in which diverse maize lines were grown under high and low nitrogen conditions. The authors surveyed the root microbiome and performed several different analyses in search of evidence of host control over the microbiome. These analyses include clustering the 16S amplicon sequencing data into 'microbial traits', searching for evidence of heritability and selection, GWAS and transcriptomics. The final section of the paper looks for patterns across the different analyses and highlights a single microbial trait, a corresponding GWAS hit and evidence that one allele correlates with both microbial abundance and plant health.

    The strengths of this manuscript include an impressive dataset, nicely compiled figures and an impressive compilation of distinct analyses to reduce complexity and reveal interactions. This manuscript is interesting on its own and will also serve as a template for future analyses on similar or distinct datasets. In several places, this manuscript could be improved by more clearly articulating the logic behind choices made for the analyses (see specific comments below).

    1. How the 'microbe traits' are defined is critical for interpreting the rest of the paper and so needs to be explained more clearly in this manuscript. From the text, I interpret that ASVs were first clustered into genera and then further clustered based on differential abundance between high and low N in both years. If so, this would seem to exclude a bunch of microbes that were specific to the 202 genotypes planted only in the second year. However, Supp Fig 2 shows a handful of 'microbe traits' that are not differentially abundant (grey). Please clarify. Also, the authors should clarify how the 'names' for each microbe trait defined.

    2. Throughout the manuscript, I expected more consideration of which microbe traits overlap between the different analyses. This comes in only at the final figure and so we never get a sense of how the different analyses overlap. For example, 150 microbe traits are defined and assessed for heritability and evidence of selection. I was not clear on whether the traits that were found to be heritable overlapped with those under selection (as one might expect). Lines 387-395 and Supp Fig 5 attempt to synthesize the different analyses but should be expanded to help the reader understand overlap between the analyses.

    3. I gather that the authors performed the GWAS separately for each 'microbe trait' and nitrogen condition, but then searched for 'hot spots' where the data from different microbe traits was considered in a pooled manner. The logic behind this decision is not clear to me. Why would we expect different microbe traits to be co-localized in the genome? Line 346-349 indicates that no plant loci were found to associate with microbe 'traits' under both nitrogen treatments and speculate as to why (dynamic interactions or not enough statistical power). The traits were defined based on robust differential abundance between nitrogen treatments and if I understand correctly, the GWAS was run separately for each trait and nitrogen treatment, so it seems logical that this method would only yield microbe trait associations that are differential between nitrogen treatments. If I did understand this correctly, I recommend emphasizing this point as it seems to indicate that the methods are working as expected.

    4. Having access to expression data for 298 genotypes is amazing. It would seem logical to try to more directly connect the MAPLs and microbe traits with this expression data. Do the lines that show association with the microbes also show higher expression of the corresponding gene? The authors generated additional RNAseq data from 2 week old plants from 4 genotypes but the logic for the selection of these lines is missing and I am not sure about the relevance of this since the samples were collected from young plants. Is a nitrogen treatment effect observable at 2 weeks? The authors conclude that the gene expression data is consistent with host control over root microbiome (line 371-373) but, as is, I'm not convinced that this analysis supports that statement. Fig 3C is striking on its own, but based on panel B, I suspect that a similar pattern would be observed for 'third leaf' and 'germinating shoot' so it is harder to make a direct connection with the microbe traits.

    5. The authors report that 62 microbe traits associated with canopy coverage, a very exciting result! However, again, this confuses me based on how the microbe traits were defined. To be considered a microbe trait, the microbes had to show differential abundance across the treatments. The logic for how this could manifest in phenotypic changes in both treatments needs to be elaborated.

    6. The final figures summarizes correlations for one microbe trait across the different analyses and looks very promising, especially for noisy field data. The authors are careful to not overstate this finding, perhaps a bit too conservative. They see a significant correlation between microbe abundance and canopy coverage that also correlates with allele frequency. The difference in canopy coverage by allele frequency is not significant, but shows a similar trend and this is not necessarily surprising given all the other factors that will influence this one trait. I expected a comment on gene expression of the genes in the locus and perhaps a peek at the other plant traits to see if any of them also show a similar trend.

  4. Reviewer #2 (Public Review):

    This study reports on a field experiment done with 230 maize genotypes under 2 nitrogen fertilization regimens. The abundance of each taxon within the root microbiota of over 3000 samples were used to perform two separate analyses: (i) a genome-wide association study (GWAS), using the abundance of each taxon as a quantitative trait and (ii) a phenotypic study, identifying microbial taxons correlated with different plant phenotypes (mainly canopy coverage).

    This is an impressive manuscript that will help pave the way for a better understanding of heritability in the plant microbiota, which is a major open question in our field. In massive screening studies such as this one, it is often challenging to deliver a concise take-home message or a mechanistic insight. Here, the actual functions of "MAPL" genes are not looked at in detail, besides looking at their expression profiles across the plants (which is an interesting analysis in and of itself). I therefore urge the authors to make their data, in the form of supplementary tables, as accessible as possible for follow-up studies on the loci that they detected.

    This manuscript will provide a valuable and nearly unprecedentedly large dataset on the genetic control maize has on microbiome assembly.

    One methodological aspect that is unclear here, is how were the 150 "rhizobiome traits" defined? This seems like a taxonomic classification, but in this version it not made clear how it was performed.

    A second and related comment, is why are the only "rhizobiome traits" that the authors consider taxon abundances? Other quantities can be extracted from this dataset: alpha diversity, beta diversity (e.g the 1st and 2nd axes of the ordination). Perhaps a measure of absolute abundance can also be extracted from the ratio of bacterial to plastid/mitochondrial reads. These would provide a more holistic picture of how plant genotype influences the microbiome composition.

    There is some disconnect between the GWAS and the phenotypic comparisons. Indeed, the authors show that heritability is correlated with the correlation with canopy coverage, but I would assume that much of the phenotypic differences observed in the first place are a result of the genotypic variability. However, no unified model that accounts for the relationship between plant genotype, microbiome composition and plant phenotype is attempted.
    The relative abundance and the prevalence of ASVs is essentially ignored, beyond the initial screening described in the appendix. It would be important to know for example if heritable taxa are also abundant ones?

  5. Reviewer #3 (Public Review):

    Strengths:

    The choice of 230 genotypes from a well-known maize diversity panel, with accompanying SNP genotype data, was a good one for this purpose, given the focus on selection during intensive breeding during the 20th century in heavily fertilized conditions.

    The very large dataset (N>3000 with replication of 230 genotypes) is a useful source of information on maize rhizosphere bacterial microbiomes, and the availability of the host genotype SNP data is an especially useful and unusual feature. The authors used a relatively newly developed (2018) Bayesian computational approach to characterize genetic architecture of rhizosphere composition, which is an interesting advance in the microbiome field. The same tool makes inferences about whether each SNP underlying rhizosphere features shows signatures of past selection (inferred as the variation in effect size relative to the minor allele frequency).

    Weaknesses:

    The BayesS results classifying rhizosphere-related SNPs as under positive, negative, or no selection appear to be over-interpreted. First, it is not clear that this method is meant for comparing current patterns of selection in contrasting environments (as in this N+/N- experiment), but rather for detecting signatures of selection in the distant past. Second, it IS clear that this method only reveals signatures of selection on a locus or SNP, and cannot confirm that selection is acting on a particular trait. The Zeng et al. 2018 paper states this quite clearly. The authors of this manuscript did not attempt to rule out that the loci classified as under selection do not have pleiotropic effects on (or are linked to) traits other than rhizosphere microbes. Occam's razor in this case suggests that these loci control root traits that are important for plant survival and also happen to affect microbiome composition. No functional benefit of these microbes was demonstrated beyond correlations with plant phenotypes.

    It is unclear how substantially different the N+/N- treatments were from each other. The entire experiment followed commercial corn, so presumably all plots had been fertilized within the past year. The soil chemical profiles were not subsequently tested, and basic analyses (such as comparison of plant growth in the two treatments) are missing. Furthermore, the assumption that microbiome differences between fertilization treatments are driven by some activity of the plant host (lines 94-99) is not justified - direct responses of the microbes to N addition would almost certainly be reflected in the rhizosphere, since rhizosphere microbiomes are almost entirely derived from the surrounding soil. No bulk soil samples were collected as controls to rule this out.

    The authors report how various patterns differ between the N+ and N- treatments- for example, more rhizobiome features had nonzero heritability in N- than in N+. However there is no statistical support for this apparent difference, i.e. no direct test of heritabilities in the two treatments. Nor does they test for possible differences (or lack thereof) in the magnitude of heritability between treatments. This incomplete style of analysis was repeated several times in the paper, e.g. for comparing patterns of selection between treatments, and for comparing correlations between rhizobiome features and plant traits between features.

    The methods sections contain inconsistencies and omissions that made it difficult to evaluate some of the claims. For instance, lines 147-148 describe collection of rhizosphere samples from entire root crowns, but the appendix (lines 774-778) describe collection of rhizosphere from roots that fit in 50 mL tubes. So it is unclear which part of the root crown was actually used, and whether the focal root type was consistent for all samples. Similarly, the appendix states that the B73xMo17 check genotype was used to correct for small-scale geographic differences (747-748), but no additional detail is provided nor are the results of this process reported. In general, the descriptions of statistical analyses lack important details. For example, by definition a constrained ordination (CAP) analysis requires a formula to be specified, but this was not reported in the paper, making it impossible to interpret the meaning of the constrained axes shown in the figures. Ordinations also require the use of a distance or dissimilarity metric, the choice of which affects interpretation - the metric used in this paper was not provided.

    Finally, many of the analyses throughout the paper take the form of testing 150 different rhizobiome traits, one by one, and then reporting the number of significant results (e.g., differential abundance between N+/N-, significant heritability, selection, correlations with plant traits). This suggests a potentially severe risk of false positives due to repeated multiple testing. After the p-values are corrected for the very large number of statistical tests (using Bonferroni, FDR, or similar) many of the conclusions might change.