Quantitative Genetics of the Maize Leaf Microbiome

A major question in microbiome science is the degree to which a host organism can affect its associated microbial communities. In animals, several host-associated microbiomes are significantly affected by host genetics, especially in the gut (Demmitt et al. 2017; Gomez et al. 2017; Goodrich et al. 2014, 2016; Xie et al. 2016). In some cases specific genes have been implicated in these associations (Belheouane et al. 2017; Bonder et al. 2016; Goodrich et al. 2016; Lynch et al. 2017). In plants, host genotype has been implicated in structuring the microbiomes of temperate and tropical forests (Kembel et al. 2014; Kim et al. 2012; Laforest-Lapointe et al. 2016; Redford et al. 2010), Arabidopsis thaliana (Bulgarelli et al. 2012; Horton et al. 2014; Knief et al. 2010; Lundberg et al. 2012), maize (Peiffer et al. 2013), grapevine (Zarraonaindia et al. 2015), tobacco (Knief et al. 2010), wild mustard (Wagner et al. 2016), and others. These studies have generally found significant effects of plant genotype on the associated microbial communities, though these genetic effects are usually overshadowed by environmental effects (Finkel et al. 2011; Müller et al. 2016; Rastogi et al. 2012).

Whereas soil is arguably the most complex microbial habitat on earth (Daniel 2005) and plant rhizosphere communities are only slightly less complex (Mendes et al. 2013), the leaf microbiome is comparatively simple (Bodenhausen et al. 2013; Delmotte et al. 2009). The leaf environment is actually very harsh, since plant surfaces have little water and nutrients available, and both surface and interior are exposed to environmental stresses such as temperature fluctuations and ultraviolet radiation (Vorholt 2012). Most leaf communities are dominated by a small number of bacterial clades capable of living in this environment, especially the Alphaproteobacteria, Gammaproteobacteria, Bacteroidetes, and Actinobacteria (Lindow and Brandl 2003; Vorholt 2012). These organisms subsist off a variety of carbon compounds derived from their plant host, including carbohydrates, amino acids, organic acids, sugar alcohols, isoprenes, and short-chain carbon compounds (Lindow and Brandl 2003; Vorholt 2012). Methanol is a particularly important carbon source for bacteria in the family Methylobacteriaceae, which contains many facultative methylotrophs present on leaves (Knief et al. 2012; Vorholt 2012).

Investigations into how plant genotype can shape leaf communities have focused mostly on individual gene mutants. Work in Arabidopsis indicates that the leaf community can be influenced by the jasmonic acid defense pathway (Kniskern et al. 2007), cuticle composition (Bodenhausen et al. 2014), and ethylene signaling (Bodenhausen et al. 2014). Attempts to map natural genetic variation that affects the leaf community have been more limited, but have identified the gamma-aminobutyric acid pathway in maize (Balint-Kurti et al. 2010) and both defense and cell wall integrity in Arabidopsis (Horton et al. 2014). More generally, leaf microbial communities have been shown to relate to disease resistance (Agler et al. 2016; Balint-Kurti et al. 2010; Manching et al. 2014), and to respond to drought stress (Methe et al. 2017), ultraviolet radiation (Jacobs and Sundin 2001; Kadivar and Stapleton 2003), and nitrogen fertilization (Manching et al. 2014).

Although several studies have identified significant impacts of plant genotype on the leaf microbial community (Bálint et al. 2013; Lambais et al. 2014; Sapkota et al. 2015), most of these have been comparisons of a small number of individuals and/or been measured across multiple species instead of within a species (Arabidopsis being the chief exception [Agler et al. 2016; Horton et al. 2014]). A better understanding of how within-species genetics affects microbial communities is needed to understand how these communities assemble and how they are impacted by their host.

One measure of how genetics impacts a trait is heritability. “Heritability” in this sense is simply the fraction of total phenotypic variation that is due to genetics (as opposed to the environment, chance, etc.). Heritability ranges from 0 to 1, with 0 meaning no genetic control and 1 meaning complete genetic control. Heritability can be further broken into two types: broad-sense heritability (H²) estimates all genetic contributions to a phenotype, while narrow-sense heritability (h²) estimates only the ones that are due to additive genetic effects (as opposed to dominance, epistasis, etc.). Breeders tend to focus on narrow-sense heritability because it determines the response to artificial selection. Heritability can be estimated through several methods, including using genome-wide marker sets to estimate heritability among unrelated individuals (Yang et al. 2017a).

To investigate the role of host genetics in the composition of leaf microbial communities, we sampled the leaf microbiome across 300 diverse maize lines growing in a common environment. Our primary goals in this research were to (i) identify the extent to which plant genetics impacts common measures of the microbiome (including diversity metrics, individual bacteria, and metabolic functions), (ii) determine which of these measures is most responsive to plant genetic variation, and (iii) identify individual maize genes or loci responsible for this influence. We chose maize both because of its high genetic diversity (Chia et al. 2012) and because of its relevance to global agriculture. Over 1 billion metric tons of maize are harvested each year (37% of global cereal production) (Food and Agriculture Organization of the United Nations 2017), and it covers an estimated 188 million hectares each year—an area nearly the size of Mexico and representing 13% of global arable land (Food and Agriculture Organization of the United Nations 2017). Given that many serious diseases of maize occur in leaves (Munkvold and White 2016), understanding the forces that shape leaf communities is a first step toward identifying ways to harness those communities for global agriculture.

MATERIALS AND METHODS

Field setup and tissue sampling.

The RNA samples for this experiment were taken from a parallel experiment on maize gene expression (Kremling et al. 2018). In brief, 300 members of the extended Goodman Maize Association Panel (Flint-Garcia et al. 2005) were grown on the Musgrave Research Farm in Aurora, New York, in the summer of 2014. Plants were blocked by maturity (five blocks) and planted with delays to try to synchronize flowering as much as possible. Tissue samples were collected on two separate dates (8 and 26 August 2014); the first sampling collected from all plots that had flowered up to that point, and the second on all that had flowered after the first but before the second. A small number of check plots were sampled on both dates.

Samples were collected from the midpoint of the second leaf down from the top (Supplementary Fig. S1). An ∼0.5- to 1-cm-wide strip of leaf tissue from the edge of the leaf to the midrib was cut out with scissors, rolled, and placed into a collection tube. Three plants from each plot were bulked into a single tube, which was then flash-frozen in liquid nitrogen. All samples were collected within a 2-h period. Leaf strips were collected between 11 a.m. and 1 p.m. (“day”), and between 11 p.m. and 1 a.m. (“night”), avoiding the previously collected plants and edge plants when possible. This separate day and night collection was intended primarily to capture circadian gene expression for Kremling et al. (2018); for this experiment they serve as technical replicates.

cDNA library preparation and sequencing.

RNA was extracted using TRIzol (Invitrogen, Carlsbad, CA) with Direct-zol columns (Zymo Research, Irvine, CA). RNA concentration was determined with QuantiFluor RNA dye (Promega) and concentrations were normalized to 100 ng/μl. RNA was then converted to cDNA using ProtoScript II reverse transcription (New England Biolabs) according to manufacturer’s instructions (1 μl template per reaction), except the enzyme amount was cut in half and the incubation time extended to 60 min. The RT primer was 1391R (5′-GACGGGCGGTGWGTRCA-3′) (Lane 1991). cDNA was then PCR-amplified in a 30-μl reaction using 2× Hot Start Taq DNA Polymerase master mix (New England Biolabs), 0.5-μl cDNA reaction as template, and the chloroplast-discriminating primers 799F (5′-TCGTCGGCAGCGTCAGATGTGTATAAGAGACAGACCMGGATTAGATACCCKG-3′) and 1115R (5′-GTCTCGTGGGCTCGGAGATGTGTATAAGAGACAGAGGGTTGCGCTCGTTG-3′) (Shade et al. 2013), which target the V5/V6 region of the gene. (Note that both primers include linker sequenced on the 5′ end to add Illumina adaptors in the next step.) The PCR program was 95°C for 5 min, 15 cycles of 95°C (30 s) – 55°C (30 s) – 72°C (30 s), and a final extension at 72°C for 5 min. The final PCR reaction was diluted 20-fold in a fresh PCR reaction with Nextera barcoded adapters (Illumina); cycle conditions were identical to above except only 8 cycles were performed. Five microliters of each reaction was pooled and purified using a QIAquick PCR Purification Kit (Qiagen). Libraries were submitted to the Genomics Facility at Cornell University for 2 × 200 paired-end sequencing on an Illumina HiSeq 2500. Raw sequence data are available through the NCBI Short Read Archive (SRP132042).

OTU calling.

Paired-end reads were joined using fastq-join through the join_paired_ends.py script from QIIME (Caporaso et al. 2010b) with a minimum overlap of 40 nt (‘–min_overlap 40’) and 15% maximum difference between reads (‘–perc_max_diff 15’). Operational taxonomic units (OTUs) were picked using the pick_open_reference_otus.py script in QIIME with the ‘–suppress_step4’ option. (This option skips the final de novo OTU picking step since preliminary analyses indicated it found a very large number of rare OTUs that were subsequently filtered out.) The default phylogenetic comparisons among OTUs were recalculated using Clustal Omega (Sievers et al. 2011) to align sequences because one clade of Methylobacteria consistently misaligned with the default methods, resulting in extremely long branch lengths. The raw OTU table was then filtered to include only samples with at least 10% more reads than the included blanks (ignoring one blank that had obvious nonsample contamination), and to exclude any chloroplast and mitochondria OTUs. This resulted in a BIOM file with 540 samples and 9,057 OTUs. All bioinformatic scripts and pipelines used in this analysis are available as part of the associated workflow package via Figshare (https://doi.org/10.6084/m9.figshare.5886769.v2).

OTU normalization.

OTUs were normalized to 10,000 reads by converting OTU counts to relative abundances, multiplying by 10,000, and rounding the result. (This is essentially taking the average of an infinite number of rarefactions to 10,000.) Any OTUs with counts of 0 across all samples were then filtered out, as were any samples with fewer than 10,000 counts in the original dataset. The nested pie chart of OTU abundance was created with Krona (Ondov et al. 2011) and modified in Inkscape (https://inkscape.org/en/) to improve clarity.

Higher-level taxonomy.

Counts for higher-level taxonomic clades were constructed from the raw (not normalized) OTU table by summing all read counts within a taxon (species, genus, etc.) based on their uclust-assigned taxonomy. This new OTU table was then normalized to 10,000 reads as above. The read depth for normalization was based on the raw OTU file, not the one with higher-level taxonomy, since the higher-level taxonomy includes each read multiple times (once for OTU, once for species, once for genus, etc.). This resulted in a total of 381 OTUs or higher-level clades across 425 samples.

Diversity statistics.

Alpha and beta diversity were calculated from the normalized OTU table using QIIME. All available alpha diversity metrics were used except for Michaelis-Menten, since this one caused the program to hang. Beta diversity was calculated with weighted UniFrac, unweighted UniFrac, and Bray-Curtis distance metrics, and the first five principal coordinates of each were used for downstream analysis. This resulted in a total of 49 diversity metrics across 425 samples.

Predicted metagenome.

The metagenome content of samples was predicted using PICRUST (Langille et al. 2013) according to its recommended protocol. Briefly, the normalized OTU table was filtered to remove nonreference OTUs and then normalized by 16s copy number. Metagenome content was predicted using both COG and KEGG annotations, which were then categorized by function into higher-level groupings. This resulted in 9,041 predicted metagenome annotations across 425 samples.

Maize genotypes.

Maize genotype data were taken from the Maize Diversity Project (https://www.panzea.org/). Genotypes consisted of the public Maize GBS v2.7 dataset mapped onto maize genome AGPv3 coordinates. Genotypes for the Goodman Association Panel were subsetted out of the full dataset. When more than one genotype sample was available for an inbred line, the one with the highest coverage was chosen. Genotypes were filtered to remove sites with minor allele frequencies <5% or with fewer than 100 samples with called (= not missing) genotypes. These resulted in a genotype file with 378,686 single-nucleotide polymorphisms (SNPs) across 323 taxa (samples). This genotype file is included in the workflow dataset on Figshare (https://doi.org/10.6084/m9.figshare.5886769.v2).

Calculation of best linear unbiased predictors (BLUPs).

The 425 leaf samples consisted of many plots that were sampled twice (day and night) and a few maize inbred lines that were present at multiple plots in the field. To reduce these to a single value per maize line, BLUPs were calculated using mixed linear models in R (R Core Team 2017). Specifically, the following model was used with the R package lme4 (Bates et al. 2015): y ∼ date + collector + maize_line, where y is each individual phenotype (OTU, diversity metric, metagenome annotation, etc.), date is the date the sample was collected on, collector is the names of the people who collected each sample, and maize_line is the name of each maize genotype (B73, Mo17, etc.) All explanatory variables were fit as random effects. Additional cofactors were tested, but all were highly collinear with date and/or collector, and were not significant once these two were included.

Since many phenotypes have highly nonnormal distributions, both the raw value and the log-transformed values were fit, and the model with smaller residual error was kept. Other models, including zero-inflated negative binomial and logistic regression models, were tested but resulted in high rates of the generalized linear model failing to converge.

Narrow-sense heritability analysis.

Narrow-sense heritability was calculated for each of the phenotypes using TASSEL (Bradbury et al. 2007). A kinship matrix was calculated from the maize genotype file and fit in a mixed linear model with flowering time as a covariate. (The heritability of flowering time itself was calculated without covariates.) Narrow-sense heritability was estimated as the amount of BLUP phenotypic variance explained by the kinship matrix relative to the total phenotypic variance. Statistical significance of this estimate was determined by randomly permuting the phenotypes and recalculating heritability to determine a null distribution. For the sake of computation, 100 random permutations were run first to filter out traits with P > 0.02, and then 10,000 permutations run on the remaining traits. Empirical P values were calculated as the fraction of random permutations with heritabilities equal to or greater than the true value.

Genome-wide association (GWAS) analysis.

GWAS for each trait was performed in two steps. First, a mixed linear model (MLM) was fit for each chromosome including both flowering time as a covariate and a kinship matrix constructed from all other chromosomes. (Meaning, for example, that the kinship matrix for chromosome 1 was constructed from markers on chromosomes 2 to 10; this is the leave-one-out method of Yang et al. (2014).) Since missing values in the flowering time data would have removed ∼10% of samples from the analysis, these missing values were interpolated based on established flowering times at https://www.panzea.org/.

After running the initial MLM, empirical P values were estimated by randomly permuting genotype labels 100 times and rerunning the analysis. (Functionally, this is randomizing genotypes among samples while keeping phenotypes and kinship constant; the latter is important to preserve the covariance between phenotypes and kinship.)

GWAS was run on all 255 traits with h² ≤ 0.001, consisting of 247 metagenome traits, five OTUs or higher-level clades, two diversity metrics, and flowering time. Twenty-five metagenome traits were excluded from subsequent analyses both because of their exceptionally high heritabilities (>0.98) and because manual inspection indicated that they were likely artifacts caused by outliers.

Clustering of GWAS hits.

Since several nearby SNPs can be in high linkage, GWAS hits were grouped according to the method of geneSLOPE (Brzyski et al. 2017), picking only the most significant SNP within each linkage window. (Linkage windows were defined by a linkage disequilibrium score of >0.3 as calculated by TASSEL.) Clusters of these independent GWAS hits were then identified by dividing the genome into nonoverlapping windows of 100,000 base pairs and counting the number of unique traits with a GWAS hit in that region. (“Hits” were defined as a marker with an empirical P value below the specified threshold, such as 0.05 or 0.01.) Since many traits were correlated with each other, these counts were corrected based on the correlation among traits in each window as follows:

1. The first trait in a window adds a score of 1, representing 1 hit.

2. Each additional trait adds a score of (1 − r²_max), where r²_max is the highest correlation between it and any other trait already included in this window. For example, if trait X is the fourth trait counted in a window and has r² values of 0.65, 0.42, and 0.87 with the previous three traits, it would add (1 − 0.87) = 0.13 to the score for that window.

Simulated QTL analysis.

To simulate QTL, we first selected flowering time as a trait to base simulations off of, primarily to confirm that no real QTL were identified in the simulation. We fit a mixed linear model using a kinship matrix as the only explanatory variable and extracted the residual phenotypes from this model (that is, the phenotypic variation not explained by the kinship matrix). We then randomly permuted these residuals 10 times. For each permutation, we chose a random 5, 20, or 100 SNP markers to be new QTL, with effect sizes drawn from a geometric distribution and scaled to account for 20, 50, or 80% of residual variance. These simulated residuals were then used to regenerate simulated phenotypes, which were run through our GWAS pipeline. To save computation time, only SNPs that had a raw P value of 0.001 or better in at least one actual GWAS were used in the simulation and GWAS, the logic being that SNPs that never showed significance in any of the 255 real traits may have statistical properties that make them unlikely to be detected regardless of what QTL are simulated at them. A QTL was considered “found” if the GWAS analysis identified it with an empirical P value of 0.01 or better.

Correlations between OTUs and disease resistance.

Breeding values for gray leaf spot (GLS), southern leaf blight (SLB), and northern leaf blight (NLB) were taken from Wisser et al. (2011). Chao1 upper bound and whole-tree phylogenetic distance (PD) measures of alpha diversity were correlated with these breeding values by linear regression with the lm() function in R. Pearson correlations between breeding values for each disease and counts of OTUs present at a relative abundance of >0.01 across samples were calculated in QIIME using the observation_metadata_correlation.py script; P values were determined based on 1,000 bootstraps of the data.

Software.

The following software was used as part of this analysis.

• ● BIOM 2.1.5 (McDonald et al. 2012a).

• ● Clustal Omega 1.2.4 (Sievers et al. 2011).

• ● FASTX-toolkit 0.0.14 (http://hannonlab.cshl.edu/fastx_toolkit).

• ● GNU parallel 20130922 (Tange 2011).

• ● KronaTools 2.7 (Ondov et al. 2011).

• ● QIIME 1.9.1 (Caporaso et al. 2010b), which wraps fastq-join (Aronesty 2011), uclust (Edgar 2010), PyNast (Caporaso et al. 2010a), the Geengenes reference alignment (DeSantis et al. 2006; McDonald et al. 2012b), RDP classifier (Wang et al. 2007), FastTree (Price et al. 2010), UniFrac (Lozupone and Knight 2005), and Emperor (Vázquez-Baeza et al. 2013).

• ● Python 3.4.3 and add-on packages argparse 1.1, biokit 0.1.2, biom-format 2.1.5 (McDonald et al. 2012a), json 2.0.9, matplotlib 1.5.0 (Hunter 2007), networkx 1.11 (Hagberg et al. 2008), numpy 1.10.1 (van der Walt et al. 2011), pandas 0.16.2 (McKinney 2010), re 2.2.1, scipy 0.16.1 (Jones et al. 2016), seaborn 0.8.1, statsmodels 0.6.1. Packages not included in the standard python library were obtained from the python package index (https://pypi.python.org/) using pip.

• ● R 3.3.3 (R Core Team 2017) and packages argparse 1.0.2 (Davis 2018), car 2.1-4 (Fox and Weisberg 2011), colorspace 1.2-6 (Ihaka et al. 2015), corrplot 0.77 (Wei and Simko 2016), gplots 3.0.1 (Warnes et al. 2016), Hotelling 1.0-4 (Curran 2017), KEGGREST 1.12.3 (Tenenbaum 2016), lme4 1.1-12 (Bates et al. 2015), parallel 3.4.0 (R Core Team 2017), and plotrix 3.6-4 (Lemon 2006).

• ● TASSEL 5.2.26 (Bradbury et al. 2007).

Most analyses were run on a desktop workstation with a 3.5 GHz Intel Xenon E5-1620V3 processor (four physical cores with two virtual cores each) and 64 GB of RAM running Linux Mint 17.2. Analyses with large workloads were run on UGA’s Sapelo computing cluster using two 48-core, 256 GB RAM compute nodes. The bioinformatic pipeline for this analysis is available on Github (https://github.com/wallacelab/2014_maize_phyllosphere), and the pipeline plus support files and key intermediate files are available on Figshare (https://doi.org/10.6084/m9.figshare.5886769.v2).

RESULTS

Diversity of the maize leaf microbiome.

We sampled the bacterial communities of approximately 300 diverse maize lines growing in Aurora, New York, in the summer of 2014. Each sample consisted of 0.5 to 1 cm strips of leaf tissue bulked from three plants per plot and flash-frozen in liquid nitrogen (see Methods for details). Metabolically active bacteria were profiled by extracting total RNA from each sample and using RT-PCR to amplify the V5/V6 region of bacterial 16s rRNA using antichloroplast primers (see Methods); a similar focus on RNA instead of DNA was shown to increase power for analyzing the mouse skin microbiomes (Belheouane et al. 2017). Operational taxonomic units (OTUs) were clustered with the QIIME pipeline (Caporaso et al. 2010b) with the default (97%) similarity cutoff. This resulted in a total of 9,057 OTUs across 540 samples, with per-sample read depth ranging from 1,912 reads to 3.17M (mean 84,142; median 32,302). Due to the nature of sampling, these reads likely represent a mix of bacteria living on either the leaf surface (epiphytes) or inside the leaves (endophytes).

Previous results indicated that leaf communities are generally low-diversity (Bodenhausen et al. 2013; Delmotte et al. 2009), and our results confirm this (Fig. 1). The vast majority of 16s reads map to the Alphaproteobacteria (74%), most of them belonging to just two families: the Sphingomonads (55% of total) and the Methylobacteria (17% of total). The remaining OTUs fall into the Betaproteobacteria (4%), Gammaproteobacteria (6%), Actinobacteria (7%), Bacteroidetes (6%), Deinococcales (0.9%), candidate division TM7 (0.9%), and other groups (1.2%).

Download as PowerPoint

We defined the “core” leaf microbiome as those OTUs that were present in at least 80% of samples. Only 16 OTUs fit this category, consisting of 10 Sphingomonads (Alphaproteobacteria; 46.4% of total dataset), two Methylobacteria (Alphaproteobacteria; 10.7%), one Enterobacteriaceae (Gammaproteobacteria; 2.2%), one Comamonadaceae (Betaproteobacteria; 2.6%), and two Microbacteria (Actinobacteria; 2.8%). Together these OTUs constitute 64.7% of the total dataset. Other OTUs present at high amounts (>1% of the total dataset) but not in the core community were two Methylobacteria (Alphaproteobacteria; 2.2%), a Hymenobacter species (Bacteroidetes; 1.9%), an Erwinia species (Gammaproteobacteria 1.1%), a Chryseobacterium species (Bacteroidetes; 1.9%), and an Agrobacterium species (Alphaproteobacteria; 1.6%) (Fig. 1). These abundant OTUs make up another 8.6% of the dataset, with the remaining 26.7% of data coming from many (∼8000) OTUs at <1% overall abundance. Some of these low-level OTUs probably derive from sequencing errors or other artifacts, though many of them are probably true individuals present at very low frequency.

The mean fraction of core OTUs in each sample was 64.7% (median 66.1%), with a range from 20.6 to 88.4% (Fig. 2A). The distribution of major bacterial clades within individual samples followed the same general pattern of the overall dataset; specifically, the majority of sample reads map to Sphingomonads or Methylobacteria, with other clades less common. A few outliers are apparent, however, where normally rare clades make up the majority of reads within a small number of samples. Individual samples show this same pattern (Fig. 2B), although sample-to-sample variability is still quite strong.

Download as PowerPoint

Alpha diversity.

We performed several analyses of the bacterial alpha diversity, or diversity present within individual samples. There are many different ways to measure alpha diversity, so rather than focus on any individual measure we looked at patterns across them. We first calculated alpha diversity using all metrics available in QIIME except for Michaelis-Menten’s fit (because it hung the program without ever terminating); the Kempton-Taylor Q metric was later dropped because it generated only NA results. Many of the available metrics measure highly similar properties of the data, so we pruned them to an informative subset where no pair of metrics had a correlation ≥0.9 or ≤−0.9.) (Supplementary Fig. S2). This left seven metrics: doubletons, Strong’s index, Berger-Parker dominance, Heip’s evenness, Simpson’s evenness, Chao1 upper bound, and whole-tree phylogenetic distance.

We tested several cofactors for potential influence on these alpha diversity metrics: the date and time of sampling, the collection team, the specific plate the RNA was stored in (for batch effects), and the maize subpopulation (tropical, stiff-stalk, non-stiff-stalk, popcorn, sweet corn, or mixed). The cofactors explained between 0 and 10% of variance individually and between 4 and 15% of variance when combined (Supplementary Fig. S3). Hidden factor analysis (analogous to principal components analysis [Lawley and Maxwell 1962]) identified hidden factors that could explain 0 to 5% more variance, depending on the trait (Supplementary Fig. S4). A Type II sum-of-square analysis was performed to identify individual cofactors that underlie these patterns. Only one cofactor was significantly associated with a metric after Bonferroni correction (Supplementary Table S1): the collection team has a Bonferroni-corrected P value of 2.1 × 10⁻⁵ for whole-tree phylogenetic distance. Rarefaction curves of alpha diversity support higher diversity levels for two specific teams, but also for samples from tropical maize lines and plants sampled later in the season (Supplementary Fig. S5). Due to the nature of the sampling scheme these factors are highly confounded in the data, so we cannot say which of them are actually driving this pattern.

Beta diversity.

To look at diversity across samples, we performed beta diversity analysis as implemented in QIIME. Beta diversity was calculated by three metrics: unweighted UniFrac, which calculates distances based on phylogenetic dissimilarity among OTUs (Lozupone and Knight 2005); weighted UniFrac, which does the same but also weights the measures by OTU abundance (Lozupone and Knight 2005); and Bray-Curtis dissimilarity, an ecological measure comparing species counts in two samples (Bray and Curtis 1957).

We compared the principal coordinates of these three datasets to possible confounding factors to identify technical artifacts that could be driving the apparent community structure (Supplementary Fig. S6). Neither Bray-Curtis nor weighted UniFrac show any clear association with batch effects. The unweighted Unifrac measure, however, shows samples clearly dividing into two groups corresponding to the two batches samples were processed in, which are also strongly but not completely correlated with day versus night samples. Unweighted Unifrac is more sensitive to rare OTUs than the other metrics, and the lack of similar clustering for weighted UniFrac and Bray-Curtis dissimilarity implies that the more common OTUs were not affected by batch processing. Multivariate analysis of variance of the first 10 principal components confirms that the RT-PCR batch has the strongest effect on unweighted UniFrac principal coordinates, whereas both Bray-Curtis and weighted Unifrac were most strongly correlated with the collection team (Supplementary Table S2). Although many cofactors show some level of statistically significant association with beta diversity metrics, visual inspection shows little practical clustering for most of them.

Only the weighted UniFrac measure showed significant correlation with bacterial clades (Supplementary Table S3). Correlation values for the most likely clades driving the first five principal coordinates are shown in Table 1, and they are graphed for the first three principal coordinates in Figure 3. In several cases multiple taxonomic levels are highly correlated with a single principal coordinate, usually because higher-level clades consist almost entirely of reads from a single lower-level clade. Not surprisingly, the two most abundant clades—Methylobacteria and Sphingomonads—have a large role in driving community structure. Principal coordinate 1 (PC1) appears to be driven by the family Methylobacteriaceae (r² = 0.845), while PC3 is driven by the Sphingomonadaceae (r² = 0.552). (Although PC3 is slightly more correlated with the Gammaproteobacteria–r² = 0.571–it seems more likely that the much more abundant Sphingomonads are actually driving that principal coordinate). Surprisingly, PC2 is driven by the Bacteroidetes (r² = 0.602), a relatively low-abundance clade that comprises only 6% of the overall dataset. It does, however, contain the two most-abundant noncore OTUs in the dataset (Fig. 1), though one of these (Hymenobacter) also appears to be the major driver of PC4 (r² = 0.699). Finally, PC5 appears to be driven by the Actinobacteria (r² = 0.609), the other major phylum in the data and which contains two core OTUs (Fig. 1).

TABLE 1 Bacterial clades likely driving beta diversity

Principal coordinatea	Clade	r2
1	Methylobacteriaceae (family)	0.845
2	Bacteroidetes (phylum)	0.602
3	Sphingomonadaceae (family)	0.552
4	Hymenobacter (genus)	0.699
5	Actinobacteria (phylum)	0.609

^a All principal coordinates are from weighted UniFrac analysis.

View as image HTML

Download as PowerPoint

Consistency within field plots.

Prior research has shown that individual leaf-dwelling species can have patchy distributions (Lindow and Brandl 2003; Monier and Lindow 2004). It was unclear whether our sampling scheme would capture enough leaf surface to be representative of the entire field plot. To determine this, we compared paired day-night samples within each plot to other plots in the experiment. Within-plot samples are visibly more similar to each other than to neighboring plots (Fig. 4). We further analyzed the distance between samples through various beta-diversity metrics. These show that within-plot sample pairs are significantly (P < = 0.01) more similar to each other than they are to corresponding samples from neighboring plots, to any nonneighboring plot, or to an equivalently sized random subset of nonneighboring plots. The only exception was among samples taken on 26 August and measured with weighted UniFrac distance, in which case the within-plot and neighboring-plot distances were not statistically different (Supplementary Fig. S7). The distances to neighboring-plot samples were significantly smaller than distances to random nonneighbor plots in five of nine comparisons, indicating there may be some spatial structuring within the field samples.

Download as PowerPoint

Heritability of microbiome community characteristics.

One of our chief goals was to estimate the extent to which the genetics of maize shape its leaf microbiome. For this, we estimated the heritability of alpha and beta diversity metrics using the BLUPs from the data. (Some data were log-transformed due to giving a better fit; see Methods.) Heritability in this sense does not refer to vertical transmission but rather to how much of total phenotypic variation is due to plant genetics. In essence, we are treating characteristics of the leaf microbiome as a plant phenotype and estimating how much control the plant exerts over the trait. When testing the overall microbial community instead of individual microbes, this is the same as the community heritability (H²_c) advocated by van Opstal and Bordenstein (2015), although our methods estimate narrow-sense heritability (h², the proportion due to additive genetic effects) instead of broad-sense heritability (H², the proportion due to total genetics). Importantly, narrow-sense heritability can be estimated in unrelated individuals using just their genetic and phenotypic information, even in the absence of full replications or family relationships (Yang et al. 2017a).

We tested all alpha diversity metrics available in QIIME, along with the first five principal coordinates from Bray-Curtis, unweighted UniFrac, and weighted UniFrac beta diversity analyses. Each metric was individually fit as the response variable in a mixed linear model in TASSEL (Bradbury et al. 2007). Narrow-sense heritability (h²) was estimated by including a kinship matrix generated from public genotype data on these maize lines and estimating the proportion of total variance explained by the kinship matrix. Alpha diversity metrics were extremely poorly heritable, with values ranging from 0 (most traits) to 0.05 (log-transformed Strong’s dominance index) (Supplementary Table S4). Beta diversity was significantly more heritable, with maximum heritability of 0.598 for weighted Unifrac PC1, 0.521 for unweighted UniFrac PC2, and 0.19 for Bray-Curtis PC1 (Supplementary Table S5).

To determine the significance of these heritability scores, we performed a permutation analysis where the phenotype data were randomly shuffled relative to the kinship matrix and the heritability recalculated. We initially screened all metrics against 100 random permutations; metrics with an empirical P value of ≤0.02 were then rerun against 10,000 random permutations to more precisely estimate the null distribution. Metrics with a final empirical P value of ≤0.001 were considered significant. (That is, at most 1 random permutation in 1,000 had a heritability as high as or higher than the actual data.) None of the alpha diversity metrics passed the initial filter, while only two beta diversity metrics passed the second: weighted Unifrac PC1 (h² = 0.598, P = 0.0003), and unweighted Unifrac PC2 (h² = 0.521, P = 0.0003).

Heritability of bacterial OTUs and clades.

Using the same methods as for alpha and beta diversity, we performed heritability analysis on 185 individual OTUs and 196 higher taxonomic units (species, genus, etc.). We tested several methods to deal with the nonnormality of many OTU counts, including log transformation and several variations on negative binomial regression. In the end, we chose log-plus-1 transformation of normalized OTU counts because it was the most robust across OTUs. (The various negative binomial variants had a high tendency to fail due to the generalized linear model not converging.)

Most taxa are not significantly heritable, with only two OTUs and three higher-level taxonomic clades showing significant heritability at an empirical P value of ≤0.001 (Fig. 5). All five of these are in the class Alphaproteobacteria; four of the five are in the order Rhizobiales, and three are in the family Methylobacteria. An additional 376 OTUs or higher-level clades did not show significant heritability, indicating that the majority of taxa living on the maize leaf surface are minimally influenced by the genotype of their host.

Download as PowerPoint

Of the 16 core OTUs, only one showed significant heritability (h² = 0.356, P = 0.0008): OTU 575956, an unnamed Methylobacterium. The heritability of other core OTUs ranged from 0 to 0.411, and only two others come close to statistical significance: an unnamed Comamonadaceae at h² = 0.411, P = 0.0013, and Microbacterium chocolatum OTU 1045797 at h² = 0.333 and P = 0.0027 (Supplementary Table S6).

Heritability of inferred metagenome content.

It is increasingly recognized that the functional capacity of a microbiome may be both more important and more consistent than its taxonomic makeup (Human Microbiome Project Consortium 2012). We used PICRUST (Langille et al. 2013) to infer the metagenome content of the maize leaf microbiome using annotations from both the Kyoto Encyclopedia of Genes and Genomes (KEGG Orthology, KO) (Kanehisa and Goto 2000) and Clusters of Orthologous Genes (COG) (Tatusov et al. 2000). We then calculated the heritability of the inferred metabolic annotations as above, including binning individual metabolic annotations into groups according to their respective hierarchies.

The raw number of metagenome annotations (9,041) was much larger than either the diversity metrics (48) or the OTUs and taxonomy (381). A correspondingly larger number of metagenome annotations were identified as significantly heritable (empirical P value of ≤0.001), with 120 KO terms and 127 COG terms (Supplementary Table S7). Twenty-five of these terms (11 KO terms and 14 COG terms) appeared to be artifacts due to highly skewed distributions (Supplementary Fig. S8) and were excluded from further analysis.

The significantly heritable terms spanned a large range of metabolic and cellular functions (Fig. 6, Supplementary Fig. S9). To identify any functions that were enriched in the leaf microbiome, we used the hierarchy of KO or COG terms embedded in the PICRUST-annotated file to search for enriched terms. Parent terms were tested by both Fisher exact test and by annotation enrichment analysis (AEA) (Glass and Girvan 2014), which uses permutation of the hierarchy to correct for highly connected terms being more likely to be found significant by chance. After correcting for false discovery rate (FDR) using the Benjamini-Hochberg method (Benjamini and Hochberg 1995), only four KO terms and zero COG terms were significantly enriched by Fisher’s exact test (Table 2). No terms were significantly enriched by AEA, although the rank-order of P values was similar (Supplementary Table S8).

Download as PowerPoint

TABLE 2 Significantly enriched metagenome terms

Term	Fisher’s exact FDRa	AEA empirical P valueb	Individual h2c	h2 empirical P valued
Methane metabolism	0.011	0.17176	0.201	0.060
Nitrotoluene degradation	0.011	0.21828	0.677	<0.0001
Glyoxylate and dicarboxylate metabolism	0.030	0.12159	0.137	0.070
Secretion system	0.041	0.15167	0.479	0.003

^a Fisher’s exact test for enrichment after Benjamini-Hochberg (false discovery rate [FDA]) correction.

^b Empirical P value from 100,000 random permutations with annotation enrichment analysis (AEA).

^c Estimated narrow-sense heritability of the individual term.

^d Empirical P value of the heritability estimate based on either 100 or 10,000 random permutations.

View as image HTML

GWAS.

After determining overall heritability, our next goal was to identify specific regions of the maize genome that affect the leaf microbiome. For this, we performed GWAS analysis for all traits with an empirical heritability P value of ≤0.001 (n = 230 traits total, including two OTUs, three higher-order taxa, two diversity metrics, 222 metagenome predictions, and flowering time as a positive control). For each trait, a mixed linear model was run separately for each chromosomes using a kinship matrix constructed from all other chromosomes (Lippert et al. 2011; Yang et al. 2014). Flowering time was included as a covariate for all traits except itself, with missing flowering times (∼10% of samples) interpolated based on published data for this panel available at https://www.panzea.org/. This same analysis was then run with 100 random permutations of genotypes (keeping phenotypes and kinship constant) to estimate empirical P values for all associations.

The results for flowering time and two other traits (log COG0181 and log K02769) are shown in Figure 7. Flowering time associations hit several genes known to be important for this trait from previous studies (Buckler et al. 2009; Dong et al. 2012). Flowering time also has many more hits than the microbiome traits, with 10 loci having empirical P values of ≤0.01 compared with only two loci for the best microbiome traits. Most microbiome traits (177 out of 229) do not have any significant hits at P ≤ 0.01 (Supplementary Table S9); 49 traits have 1 hit and 3 have 2 hits. Even if the threshold is relaxed to P ≤ 0.05 (a level likely to introduce many false positives), only 136 of the 229 microbiome traits have at least 1 hit. Simulating QTL with known locations and effect sizes resulted in a similar pattern where only strong QTL were detected (Supplementary Table S10), indicating that most QTL that affect the leaf microbial community are probably of small effect.

Download as PowerPoint

We also investigated the possibility that multiple traits are controlled by the same locus. For this, we separated the genome into 100 kb windows and counted the number of traits with an association in each window. We then adjusted this score based on the correlations among traits (see Methods for details) to generate a correlation-corrected score for each window. We identified two windows with correlation-corrected scores >2 (roughly equivalent to at least two wholly independent traits), on chromosomes 7 and 10, and both result from a large number of traits having significant associations at a single SNP marker (21 traits at 7:165,427,159 and 6 traits at 10:92,038,994). All but one of the traits at these two locations are for metagenome predictions (Table 3); the cluster on chromosome 7 also includes a hit for a beta diversity metric, namely, the first principal coordinate of weighted UniFrac analysis.

TABLE 3 Genome-wide association hit clusters

Chromosome	Number of traits	Adjusted score	Marker name	Physical position	Traits
7	21	2.26	S7_165387694	7:165427159	COG1033: Predicted exporters of the RND superfamily
					COG1139: Uncharacterized conserved protein containing a ferredoxin-like domain
					COG1164: Oligoendopeptidase F
					COG1692: Uncharacterized protein conserved in bacteria
					COG2319: FOG: WD40 repeat
					COG3184: Uncharacterized protein conserved in bacteria
					COG3257: Uncharacterized protein, possibly involved in glyoxylate utilization
					COG3328: Transposase and inactivated derivatives (log-transformed)
					COG3484: Predicted proteasome-type protease
					COG3526: Uncharacterized protein conserved in bacteria
					COG4103: Uncharacterized protein conserved in bacteria
					d-Arginine and d-ornithine metabolism
					K00824: d-alanine transaminase [EC:2.6.1.21]
					K05946: N-acetylglucosaminyldiphosphoundecaprenol [EC:2.4.1.187]
					K07003: None
					K07395: Putative proteasome-type protease
					K07401: Selenoprotein W-related protein
					K08602: Oligoendopeptidase F [EC:3.4.24.-]
					K09769: Hypothetical protein
					K14977: Ureidoglycine aminohydrolase [EC:3.5.3.-]
					Weighted UniFrac PC1
10	6	3.20	S10_92025313	10:92038994	COG0510: Predicted choline kinase involved in LPS biosynthesis (log-transformed)
					COG3131: Periplasmic glucans biosynthesis protein (log-transformed)
					COG3143: Chemotaxis protein (log-transformed)
					COG3772: Phage-related lysozyme (muraminidase) (log-transformed)
					COG4553: Poly-beta-hydroxyalkanoate depolymerase
					K02276: Cytochrome c oxidase subunit III [EC:1.9.3.1]

View as image HTML

Relationship of the leaf microbiota with disease resistance.

Previous work indicates that plant microbial communities can affect disease resistance phenotypes, such as alpha diversity being positively correlated with resistance to SLB in maize (Manching et al. 2014). Although it was not a primary focus of our investigation, we tested if the leaf communities in our data showed correlations with disease phenotypes by comparing our microbiome data to published breeding values for three major maize diseases: SLB, NLB, and GLS (Wisser et al. 2011). Pearson correlations were then calculated between breeding values and OTU counts for each OTU present at a relative abundance of >0.01 across samples.

Neither alpha nor beta diversity is significantly correlated with disease resistance (Supplementary Fig. S10). However, we did identify 10 individual OTUs that were significantly correlated with SLB, 9 with NLB, and 18 with GLS (Supplementary Table S11). The only OTUs positively correlated with disease resistance breeding value were those in the family Methylobacteriaceae, and the ones that could be classified at genus level were all Methylobacterium spp. All other significant correlations were negative.

DISCUSSION

The maize leaf community is low diversity and distinct from the rhizosphere.

Our results confirm that the maize leaf microbiome is a low-diversity community. Over 70% of sample reads map to just two families, the Sphingomonads and the Methylobacteria, with over half (52%) of all reads mapping to a single genus, Sphingomonas. These results are in line with leaf microbiome data from other species (Delmotte et al. 2009) and in stark contrast with data from the maize rhizosphere, where the overall community shows both more total OTUs and greater taxonomic diversity (Peiffer et al. 2013; Yang et al. 2017b). For example, a study that had previously looked at maize rhizosphere samples at this same experimental station showed that the Sphingomonads and Methylobacteria make up <1% of the rhizosphere (Peiffer et al. 2013). This is in line with other studies that show little correlation between soil and leaf microbial communities, especially later in the season (such as when these samples were taken) (Copeland et al. 2015; Kim et al. 2012). Wind dispersal from other plants appears to play a strong role in leaf bacterial community dynamics (Lindow and Andersen 1996; Lymperopoulou et al. 2016), so that the local soil may not even serve as the inoculum for much of the leaf community.

Drivers of leaf community structure.

One of the major goals of this research was to identify the primary drivers of the maize leaf community. The first task was to rule out any confounding from technical errors. Both Bray-Curtis and weighted UniFrac measures of community diversity show significant associations with possible confounding factors, especially the collection team. It is possible that this reflects operational biases for each team. It is also possible, however, that the “collection team” variable is instead just the best representation of some latent variable in the data because it breaks the data into the most number of groups (13), and these groups tend to be confounded with cofactors like flowering time and field location due to the logistics of sample collection. These possibilities are not mutually exclusive, and with the current data we cannot wholly tell which is correct. Regardless of what process it represents, the collection team was included as a cofactor in the generation of BLUPs to factor it out of downstream analyses.

Figure 4 shows that our paired day and night samples show relatively strong consistency within a field plot, indicating that our sampling method was sufficient to get a good representation of the plant leaf community (specifically, that of upper leaves). The visible differences between neighboring plots of different genotypes (Fig. 4, top) would indicate a potentially strong effect of plant genotype. However, the B73 check line, which was planted at multiple locations throughout the field, shows no more plot-to-plot consistency than neighboring plots with different genotypes (Fig. 4, bottom). This indicates that there is probably a strong environmental and/or stochastic component to leaf community composition. Since all plants in this experiment grew in the same field and within a few dozen meters of each other, we expect the environmental variation from plot to plot to be minimal, so much of the nongenetic variation in this data likely comes from stochastic noise. This is in line with previous results showing that leaf communities are shaped by stochastic events early in colonization (Maignien et al. 2014).

No alpha diversity metrics show any significant heritability, indicating that maize genotype does not have a major effect on the diversity of its leaf community. This is surprising, as previous work in a maize recombinant inbred line (RIL) population showed a significant relationship between plant genetics, community diversity, and disease resistance (Balint-Kurti et al. 2010). Specifically, that analysis mapped loci influencing Simpson diversity to six bins in the maize genome. In contrast, the Simpson’s diversity index of our data has a heritability of only 0.008 and an empirical P value of 0.38, so it was never included in the GWAS mapping analysis. (Although this heritability is very low, it is actually among the highest for alpha-diversity metrics.) Many experimental details are different between the two analyses (RILs versus a diversity panel, looking at terminal restriction fragment length polymorphisms versus 16s rRNA, etc.), so we cannot say which is more likely to be correct overall.

Most beta diversity measures were also not significantly heritable, with the exception of UniFrac principal coordinates. Weighted UniFrac PC1 showed significant heritability (h² = 0.598) and is driven primarily by Methylobacteriaceae (Fig. 3), a bacterial clade which is also significantly heritable on its own (Fig. 5). This indicates that plant genetics can significantly impact the degree to which these organisms colonize the leaves, although they were one of very few clades showing significant heritability. Methylobacteria are thought to subsist off of short-chain carbon molecules that leak out of leaves (Lindow and Brandl 2003; Vorholt 2012). The association between Methylobacteria and plant genetics may thus relate to different metabolic pathways within the plants that generate these compounds, or to different cell wall structures that make leaves more or less permeable to them. The only other diversity metric with significant heritability was unweighted UniFrac PC2. Although it is only weakly associated with specific bacterial clades, its strongest correlation is also with the Methylobacteriaceae (r² = 0.341) reinforcing the idea that this clade is one of the few with a robust effect from maize genotype.

Unfortunately, we could identify no significant GWAS hits for Methylobacteria at any clade level with an empirical P value cutoff of P ≤ 0.01. The first principal coordinate of weighted UniFrac analysis did show a hit on chromosome 7 (7:165,427,159, maize genome AGPv3 coordinates) as part of a cluster of metabolic GWAS hits; see the next section for more details on this cluster. If, however, we relax the P value threshold to empirical P ≤ 0.05, we do identify hits for two separate Methylobacteria OTUs: OTU 575956 (unknown species) at 8:8,617,635 with an empirical P value of 0.04, and OTU 591699 (Methylobacterium adhaesivum) at 6:31,420,753 with an empirical P value of 0.05. The hit on chromosome 8 is in GRMZM2G396541, which is annotated as leucine-rich repeat extensin-like protein 5, a protein that in Arabidopsis is implicated in cell wall formation (UniProt:Q9SN46). GRMZM2G396541 is also annotated as sulfated surface glycoprotein 185, which is involved in the extracellular matrix of Volvox (Hallmann 2003). Meanwhile, the hit on chromosome 6 is in an ABC transporter of unknown function (GRMZM2G099619). Although these annotations could logically have a role in affecting the leaf surface microbiome, the high empirical P value for these hits means they should be interpreted with caution.

Host genetics may influence leaf microbiome metabolism more than taxonomy.

It is increasingly thought that the metabolic functions of a microbial community are more important than its taxonomic identity (Coles et al. 2017; Louca et al. 2016). That is, communities which appear very different based on bacterial OTUs could actually appear quite similar when measured by the metabolic functions present. Along these lines, our analyses showed higher and more significant heritability for inferred metabolic functions than for diversity metrics or individual OTUs. Twenty metagenome terms had empirical P values of <0.0001, meaning zero out of 10,000 random permutations showed a heritability that high or higher. Unfortunately, many of these terms are poorly characterized. Of the ones with assigned functions, some of the most heritable include alkyldihydroxyacetonephosphate synthase (K00803, involved in ether lipid metabolism; h² = 0.70), urea carboxylase (K01941, h² = 0.70), nitrotoluene degradation (h² = 0.68), a tRNA amidotransferase (K02433; h² = 0.67), a circadian clock protein (K08481; h² = 0.66), NifD nitrogenase (COG2710, involved in nitrogen fixation; h² = 0.64), a glycan biosynthesis gene (K05946; h² = 0.61), ureidoglycine aminohydrolase (K14977, involved in purine catabolism; h² = 0.59), and isoquinoline 1-oxidoreductase (K07302 and K07303 for the alpha and beta subunits, respectively, and involved in isoquinoline degradation; h² = 0.53 and 0.52). Some of these functions make sense in terms of leaf microbiome biology, such as a protein likely involved in glyoxylate utilization (COG3257; h² = 0.60), and two proteins implicated in DNA repair (COG1059 and COG1743; h² = 0.50 and 0.49, respectively).

A Fisher’s exact test for enrichment of terms among the highly heritable (P ≤ 0.001) terms identified four higher-level terms that are significantly enriched: “methane metabolism,” “nitrotoluene degradation,” “glyoxylate and dicarboxylate metabolism,” and “secretion system.” Both methane metabolism and glyoxylate metabolism involve the metabolism of short (1- or 2-carbon) carbon compounds. These compounds are known to be an important energy source for leaf-associated bacteria, especially Methylobacteria (Lindow and Brandl 2003). Nitrotoluenes are not normally considered part of the plant environment but are involved in the manufacture of agricultural chemicals (Dugal 2005), so the enrichment of this category may reflect management practices. Interestingly, nitrotoluene degradation was the only higher-level term that was both statistically enriched (P = 0.01 after FDR correction) and significantly heritable on its own (h² = 0.677, P < 0.0001), making it among the most robust hits in the dataset. Finally, membrane secretion systems could be involved in many important processes, including intercellular signaling, manipulation of the host, or biofilm production, but this dataset lacks the detail to know which one(s) may be most significant.

The fact that only four higher-level metabolic terms were found significant by Fisher’s exact test and none by annotation enrichment analysis (which was designed specifically to correct for biases in Fisher’s exact test applied to hierarchical metabolism enrichment) indicates that these results should be taken with caution. This is especially true since these results are based off PICRUST-imputed metagenome content instead of actual metagenome data. For now, these results are suggestive of what processes might be significantly controlled in the leaf microbiome, but further work will be needed to confirm them.

Few leaf microbiome traits show robust GWAS associations.

Our inclusion of flowering time as a positive control confirms that our analysis pipeline is capable of identifying robust, biologically relevant hits. However, the small number of hits for leaf microbiome traits (most with 0 or 1 at empirical P ≤ 0.01) indicates that these traits are more difficult to map than flowering time. Several factors could explain this. Even the most heritable microbiome traits have h² values around 0.7, compared with 0.92 for flowering time, so that the microbiome traits have intrinsically lower power for mapping. They may also have a more distributed genetic architecture, or at least fewer large-effect alleles, so that each individual locus explains less variance and is harder to identify. This may explain why nitrotoluene degradation, one of the most robustly associated metagenome traits, has no GWAS hits. When we simulate artificial QTL in this dataset, we find a similar pattern of only the strongest QTL being detected, which supports the idea that most QTL for leaf microbiome traits are undetected because they are of small effect.

Of the hits we do identify, many are close to or within known genes. Unfortunately, the majority of these genes have no known function, making interpretation of their association difficult. In addition to individual hits, however, we also identified two strong hit clusters on chromosomes 7 and 10 that are associated with multiple traits. The first is located at position 7:165,427,159 (maize genome AGPv3 coordinates) and includes 21 traits with a correlation-corrected score of 2.26, roughly equivalent to just over two completely independent traits mapping to the same location. The other is at 10:92,038,994 and includes 6 traits with a corrected score of 3.2.

Linkage disequilibrium (LD) analysis indicates that the hit cluster on chromosome 10 may be part of an extended haplotype, as there is at least one SNP ∼100 kb away in very high LD (Supplementary Fig. S11). This hit is close (∼65 kb) to only a single gene of unknown function, GRMZM2G158141. In addition, the minor allele frequency at this location is very low (only 5 of 102 individuals with nonmissing genotypes), which may increase the odds of this being a spurious association due to stochastic sampling.

In contrast, the hit cluster on chromosome 7 is much more robust, with 25 minor alleles out of 230 individuals. It is also very high resolution, with no other SNPs in significant linkage disequilibrium within 10 megabases or elsewhere in the genome (data not shown). This hit is located just upstream (∼1 kb) of GRMZM2G031545, a housekeeping gene (elongation factor 1-beta) involved in protein synthesis. The Maize Gene Expression Atlas (Stelpflug et al. 2016) (accessed via MaizeGDB.org) indicates that this gene is indeed a housekeeping gene with almost universally high expression across tissues; it is unknown why or how such a gene would affect the leaf microbiome. We examined the relationship between GRMZM2G031545 and several of the traits associated with it using maize expression data from the same samples (Kremling et al. 2018), but the correlations were extremely weak (Supplementary Fig. S12). Most of the other nearby genes have no known function, making it difficult to speculate what biological processes may be involved in this association.

Implications for disease resistance.

We identified several putative correlations between bacterial OTUs in our dataset and published disease resistance values for this panel (Wisser et al. 2011). The only positive relationships with disease resistance were with Methylobacteria. Beneficial interactions between plants and Methylobacteria are well-documented, and these bacteria have been shown to induce disease resistance via production of antimicrobial compounds (Balachandran et al. 2012; Dourado et al. 2015; Ryan et al. 2008), competition with pathogens (Araújo et al. 2002; Berg 2009; Lacava et al. 2008), and induction of induced systemic resistance (Dourado et al. 2015; Madhaiyan et al. 2006; Nigris et al. 2013; Pavlo et al. 2011; Tamosiune et al. 2017). The existence of a significant association between Methylobacteria and disease resistances is intriguing, especially since the data came from different fields in different years. Confirmation of this relationship would require additional experiments, including controlled experiments to determine if greater Methylobacteria counts are a cause or a result of disease resistance.

Conclusions and future directions.

Our primary goals were to determine how much maize genetics affects leaf microbiome characteristics, which characteristics are most influenced by maize genetics, and what maize loci are responsible for these associations. Our results indicate that only a subset of bacteria in the maize leaf microbiome are significantly impacted by the genetics of their host. This is in line with a prior analysis of the maize rhizosphere (Peiffer et al. 2013), where plant genotype had small but statistically significant effect on the rhizosphere community. In our analysis, we find robust signals for the host plant influencing a few key taxa and traits, including the abundance of Methylobacteria and the metabolism of short-chain carbon molecules and nitrotoluene. In line with previous results indicating that microbiome function may be more important and conserved than taxonomy (Human Microbiome Project Consortium 2012; Langille 2018), we find more and stronger associations with inferred metagenome content than with community diversity metrics or individual taxa.

We also identified genetic loci associated with specific microbiome traits, although only a few loci were associated with any given trait. Based on these results, we suspect that this panel is likely underpowered for leaf microbiome traits. Although a panel of this size (n = 246 individuals after filtering) is adequate for traits with several strong QTL, like flowering time (Fig. 7), it appears that larger panels will be needed to map most microbiome-related traits. The lower heritabilities of these traits are partly responsible for this lack of power. Microbiome-related traits may also be controlled by many genes of small effect, making it more difficult to identify any individual locus as having an effect on the phenotype. We would therefore recommend including at least 350 diverse lines in future panels, and up to 400 or 500 lines if resources permit.

Although we tried to be thorough in our analysis, no single experiment can fully describe the entire leaf microbiome system. There are several lines of inquiry that could complement this analysis, including looking at different leaf communities (such as comparing abaxial versus adaxial surfaces versus endophytes), looking at different leaves down the plant or at different developmental stages, examining stalk communities, etc. One could also explore if host genetics can modify the evolution of leaf communities over time, since these communities are known to change over the season (Manching et al. 2018). Finally, it will be important to expand these analyses to nonbacterial members of the microbiome (fungi, oomycetes, etc.) to further understand leaf community assembly and function.

An important take-home from our analysis is that the leaf community is likely amenable to manipulation through farmer management. This would make community manipulation much faster and easier to accomplish than if new plant varieties had to be bred with specific genetic characteristics. A full systems-level understanding of the leaf community will be necessary to understand what management practices may be beneficial to the community, and how the community as a whole could be harnessed to improve global production.

ACKNOWLEDGMENTS

We thank members of the Buckler lab for help collecting samples (especially S.-Y. Chen, M.-K. Su, J. Pardo, N. Lepak, and J. Budka), the Genomic Diversity Facility at Cornell for help with automation, and members of both the Buckler and Wallace labs for commentary and discussions during preparation. We also thank W. A. Walters for help in designing the 16s rRNA amplification.