Genome-wide association study of Fusarium ear rot disease in the U.S.A. maize inbred line collection

Line means and heritability

Genome-wide association mapping of Fusarium ear rot resistance

Background polygenic effects modeled by K accounted for 31% of the variation among entry means in the full inbred association panel analysis and 42% of the entry mean variation in the balanced subset inbred association panel (Table 2). Principal component decomposition of K revealed little association between mean rot scores in the inbred association panel and large-scale population structure (Figure 1). In the topcross analyses, K accounted for 31% of the variation among B47 topcross entry means and 39% of the variation among PHZ51 topcross entry means (Table 2).

	N a	Groups b	Compression c	${\left(\frac{{\hat{\mathbf{\sigma }}}_{\mathbf{G}}^{\mathbf{2}}}{{\hat{\mathbf{\sigma }}}_{\mathbf{G}}^{\mathbf{2}}\mathbf{+}{\hat{\mathbf{\sigma }}}^{\mathbf{2}}}\right)}^{\mathit{d}}$
Full inbred panel	1687	2100	1.18	0.31
Filtered inbred panel	734	2000	1.24	0.42
B47 topcrosses	243	1760	1.41	0.31
PHZ51 topcrosses	313	1770	1.40	0.39

^a Total number of entries included in the analysis.

^b Number of groups determined by optimum compression (note that the complete kinship matrix for 2480 lines was used for all analyses).

^c Compression level is the average number of individuals per group.

^d Polygenic additive background genetic variance divided by total phenotypic variance. This ratio was estimated in GAPIT by fitting the kinship matrix (K) in the mixed linear model without any SNP marker effects.

From the analysis of the full inbred association panel, two SNPs (at bp 64,771,372 on chromosome 5 and at bp 19,532,465 on chromosome 9) were identified as significantly associated with ear rot resistance at a false discovery rate (FDR) of < 0.20 (Table 3; Figure 2). These two SNPs also had the highest RMIP values among SNPs across the 50 data subsamples; the chromosome 9 SNP had an association with ear rot with p-value < 10⁻⁵ in 47 of the 50 data subsamples (Table 3; Figure S1 in Additional file 1; File S6 in Additional file 1).

Chromosome	SNP physical position (bp)	p -value	FDR adjusted P-value	Minor allele frequency	Allele effect (%) a	( R 2) b	Gene containing SNP	SNP effect	RMIP
		Full inbred panel (1689 lines) analysis
5	64,771,372	8.83 × 10−7	0.089	0.07	−0.170	1.3	GRMZM2G060659	mis-sense (A/T)	0.38
9	19,532,465	8.44 × 10−8	0.017	0.15	−0.134	1.5	GRMZM2G035665	mis-sense (V/A)	0.94
		Filtered inbred panel (737 lines tested in three years) analysis
4	7,566,354	7.34 × 10−7	0.074	0.10	−0.230	2.9	GRMZM2G372364	intron variant
4	7,618,125	2.67 × 10−6	0.175	0.10	−0.225	2.6	GRMZM2G012821	mis-sense (N/D)
4	7,618,284	3.96 × 10−6	0.175	0.11	−0.205	2.5	GRMZM2G012821	mis-sense (D/N)
4	9,353,851	6.14 × 10−7	0.074	0.07	−0.254	3.0	GRMZM2G419836	3′ UTR variant
4	124,930,006	4.36 × 10−6	0.175	0.04	−0.271	2.5	GRMZM2G106752	mis-sense (L/S)

^a Allele effects are reported back-transformed to the original 0-100% disease severity scale. Effects are in reference to the minor allele.

^b R ², proportion of total entry mean variance associated with a SNP after accounting for background polygenic variance.

When the analysis was conducted on a filtered data set including only lines with data from all three years, a distinct set of five SNPs, all on chromosome 4, were identified as significantly associated with ear rot resistance (Table 3; Figure 2). No significant SNPs at FDR < 0.20 were identified from either the B47 topcross analysis or the PHZ51 topcross analysis (Figure 3), where the minimum raw P-values among SNP association tests were 1.3 × 10⁻⁵ and 2.3 × 10⁻⁵, respectively.

SNPs identified from either of the two inbred analyses explained relatively small proportions of the observed variance in entry means after accounting for the background polygenic effects (individual SNP R ² values ranged from 1.3% to 3.0%, Table 3), and each SNP also had a small allelic effect (−0.13% to −0.27% back-transformed to the original percentage ear rot scale). All significant associations had negative allelic effects, indicating that the minor allele was associated with lower ear rot (increased diseased resistance) at all loci.

Genes colocalized with associated SNPs

To gauge the resolution of associations, we inspected the local LD structure around the significant associations (Figures 4 and 5). Romay et al. [26] summarized the genome-wide LD characteristics of this panel, noting that LD tends to decay rapidly to below r ² = 0.2 within 1 kb, but that there is substantial variation around this average value among genome regions and germplasm groups. The regions around our associations on Chromosome 4 near 125 Mb and on Chromosome 9 exhibit the typical rapid decay of LD observed in diverse maize. LD was slightly more extensive around the Chromosome 5 association, with a few SNPs about 200 kb away from the significant association having r ² of about 0.5 with the associated SNP. Finally, the region on Chromosome 4 between 7.5 and 9.5 Mb had the most extensive LD, with SNPs separated by almost 2 Mb still having high LD, although much of the region between the ends of this section had much lower LD. Romay et al. [26] observed that Chromosome 4 has particularly high LD. The high LD region reported here is coincident with the interval containing the gametophyte factor 1 (Ga1) locus [29], which is under selection in the popcorn subgroup and may also be more widespread in tropical maize due to selfish gene evolution [30]. These selection effects associated with Ga1 may be involved in maintaining LD in the region.

Genes containing SNPs significantly associated with ear rot resistance were characterized using the filtered predicted gene set from the annotated B73 reference genome [31] (Additional file 1: File S7). All seven SNPs identified across both inbred association panel analyses were within predicted genes on the maize physical map, five of the seven localized to exons (all coding for nonsynonymous mis-sense variations), one to the 3′ untranslated region, and one to an intron (Table 3). The disease associated SNP on chromosome 5 was in a sucrose synthase gene (GRMZM2G060659) located in an LD block extending approximately 0.2 Mbp on chromosome 5 (Figures 4C and 5C). Examination of the lines carrying the minor allele at this locus revealed no relationship between population structure due to kernel type (namely the sweet corn and popcorn groups) and presence of the minor allele. The associated SNP on chromosome 9 was in a DNA replication factor CDT1-like gene (GRMZM2G035665) located at the end of a 0.1 Mbp LD block on chromosome 9 (Figures 4D and 5D). All five SNPs identified in the balanced subset of the inbred association panel analysis were located on chromosome 4 (Figures 4A, B and 5A, B). Four of those SNPs were located in a 1.8 Mbp region between physical positions 7,566,354 bp and 9,353,851 bp, representing a region of high linkage disequilibrium covering a genetic distance of less than 1 cM (Liu et al. 2009) (Figure 4A). The four SNPs in this region were all in high LD relationships with each other (r ² from 0.62 to 0.84; Figure 5A). Two of the SNPs in this region localized to an exon of an F-box domain gene, one localized to a thioredoxin gene, and the last localized to a gene of no known function (GRMZM2G012821, GRMZM2G419836, and GRMZM2G372364, respectively). The fifth SNP identified on chromosome 4 located at position 124,930,006 bp localized to an exon of a loricrin-related gene (GRMZM2G106752).

Heritability and genotypic correlation between experiments

The removal of lines that were not tested in all three years (consisting mostly of 953 unreplicated inbred lines that were present only in the 2010 NCPRIS collection experiment) substantially improved the entry mean-basis heritability. (${\hat{H}}_C=0.21$ in full data set versus ${\hat{H}}_C=0.61$ in filtered data set). This large difference in heritability provided justification for conducting separate GWAS on the complete and filtered inbred association panel data sets. Improved heritability of the mean values from the filtered panel will contribute to increased power of GWAS [32], but this is balanced by the loss of diversity and reduced allele replication in the subset compared to the complete set of inbreds. Analyses on the full versus filtered inbred data sets identified different genomic regions significantly associated with Fusarium ear rot resistance (Table 3). These differing results presumably reflect the tradeoffs between higher heritability and larger sample size that affect GWAS power.

Although the heritability estimate for ear rot resistance averaged across testers in the topcross experiment $\left({\hat{H}}_C=0.63\right)$ was comparable to that of the filtered inbred data set, no SNPs were identified as being significantly associated with ear rot resistance in either the B47 or PHZ51 topcross data sets. Estimates of genetic variance in the heritability calculations revealed reduced genetic variance in the topcross experiment compared to the inbred experiments (Table 1). Smaller genotypic sample size of the topcross experiment also contributes to reduced power of detection of SNP associations. In addition, genotypic correlations between inbred per se resistance and hybrid performance in the two sets of topcrosses were moderately low (r _g ≤ 0.42).

Association mapping

Two SNPs significantly associated with ear rot resistance, located on chromosomes 5 and 9, respectively, were identified in the full inbred association panel analysis, and five additional SNPs (representing two different LD blocks) were identified on chromosome 4 in the filtered inbred panel analysis (Table 3). Although all SNPs localized to genic regions, no obvious relationship exists between the predicted functions of these genes and Fusarium ear rot resistance; however the currently limited understanding of pathways contributing to resistance restricts our ability to predict what genes might be involved in resistance to this complex disease.

These SNP associations are different than those previously reported by Zila et al. [22] based on analysis a subset of 267 lines with a smaller and largely distinct set of SNPs. The closest pair of associations between the two studies were the SNPs on chromosome 5, which localized to the same genomic bin; however, they are 34 Mbp distant from each other physically, and 14.4 cM apart genetically [33]. The differences between the results presented here and those reported by Zila et al. [22] may be due to sample size and sampling of alleles and also due to differences in the SNPs tested for association. None of the three SNPs reported as associated with ear rot resistance by Zila et al. [22], located on chromosomes 1 (63,540,590 bp), 5 (30,997,717 bp), and 9 (151,295,233 bp), were present in the filtered GBS Romay et al. [26] marker set, and thus we had no potential to detect them in this study. The nearest neighboring filtered GBS SNP to each of the three SNPs reported by Zila et al. [22] were located 82 bp (raw p = 0.44), 2902 bp (raw p = 0.74), and 299 bp away (raw p = 0.11), respectively. However, the chromosome 9 SNP from the Zila et al. [22] study located was present in the original unfiltered Romay et al. [26] marker set, but a follow-up analysis of this single marker in GAPIT using the full inbred panel found it insignificant (raw P = 0.78). Finally, only the three SNPs in the LD block from 7.5 – 9.2 Mb on Chromosome 4 in this analysis colocalized with any QTL intervals identified in two biparental families by Robertson-Hoyt et al. [15]. QTL positions for Fusarium ear rot are not consistent among biparental families [15],[16], but this one QTL region on Chromosome 4 was exceptional in being identified by linkage in two families by Robertson-Hoyt et al. [15] and by association in this study.

The variability in SNP association results among different germplasm samples may be due in part to the relatively small effect sizes of the potentially many underlying causal variations, coupled with low frequency of many variants and rapid decay of LD in diverse maize germplasm. This could result in a situation where even SNPs physically close to a causal variant are not likely to be associated with enough phenotypic effect to permit their robust and reliable detection through association analysis in diverse populations. The high frequency of detection of the chromosome 9 SNP (in nearly all random subsamples of 80% of the full data set; Figure S1 in Additional file 1; File S6 in Additional file 1), and the consistency of its effect even in the filtered subsample (where although it did not pass the FDR threshold of 0.2, its raw p-value was 2.15 × 10⁻⁵, Additional file 1: Table S1), suggest that its association in very diverse maize is reliable.

The five SNPs on chromosome 4 that were detected in the filtered but not the complete inbred panel had substantially higher allele effect estimates in the filtered panel, but similar allele frequencies across panels (Table S1 in Additional file 1). The difference in these results may be due to a reduction in the influence of many line means with only a single environment observation associated with a lower heritability in the full inbred line panel, and possibly greater precision of the resulting allele effect estimates. In contrast, the two SNPs detected in the full panel had consistent allele frequencies and effect estimates across the two analyses, but simply did not have sufficient statistical significance to stand out among the hundreds of thousands of tests performed (Table S1 in Additional file 1).

Although this study used four times the number of SNP markers (200 k versus 47 k) and an association panel almost six times as large as those used by Zila et al. [22], the number of genic regions identified as significantly associated with ear rot was about the same for the two studies (four and three, respectively). Furthermore, the proportion of phenotypic variance among entry means explained on average by the K matrix across the two inbred analyses and two topcross analyses was similar to results reported by Zila et al. [22]. These results suggest that the genetic architecture of resistance to Fusarium ear rot is highly polygenic, with substantial genetic variability generated by a large number of effective variants, each with individually small effects. Even with increased marker coverage and a larger association panel, the results of this study highlight the limitations of GWAS to precisely identify allele variants with small effects on complex traits.

Marker coverage in this study is still insufficient to provide SNPs in high LD with all segregating sequence variants; Romay et al. [26] suggested that more than 700,000 SNPs would be required to tag almost all variant regions in diverse maize. Examination of the annotated genes around the significant associations reveals a number of genes nearby that contain no SNPs in our data set (Figure 5; Additional file 1: File S7), suggesting that we are likely to miss some true associations. Thus, it is possible that a further increase in marker density might reveal more SNP associations and possibly some genetic variants with larger effects. However, if the genetic architecture really is highly polygenic, then the benefit of increasing marker density on increasing the likelihood of tagging additional causal variations by LD association is likely to offset by the increasingly stringent significance thresholds imposed by the larger number of association tests conducted. The additional benefit of adding markers is also somewhat limited if most of the markers have low minor allele frequency (MAF), as is the case for the GBS markers used here [26]. The SNP associations detected in this study had minor allele frequencies ranging from 0.04 to 0.15 (missing phenotypic observations caused some markers to have MAF < 0.05 in the GWAS), compared to minor allele frequencies below 0.05 for more than half of the complete GBS marker set. Besides having low power of detection just due to reduced allele replication, rare alleles tend to be highly associated with population structure since they are usually limited to a single subpopulation, thereby further reducing their potential for trait association following correction for population structure. In this study, we removed SNPs with MAF < 0.05 to ensure reliable associations based on sufficient replication across lines. If rare alleles are a major component of the genetic architecture, however, we may have missed many important associations by dropping SNPs with low allele frequencies that would represent the best possible associations with rare functional alleles. Further studies would be required to better understand the compromises between improving reliability of results by removing rare SNPs versus potentially missing important but rare functional variants.

No significant SNPs were identified in either topcross analysis, and examination of the empirical distribution of P-values from the four analyses revealed a tendency towards higher P-values in the two topcross analyses compared to the two inbred panel analyses (Figures S2, S3, S4, and S5 in Additional file 1). Heterosis plays a significant part in Fusarium ear rot resistance, reducing both genetic variance and the mean level of disease in F₁ hybrids compared to inbred parents [34], which can reduce the ability to discriminate levels of disease resistance in topcross hybrids. Further, within a set of hybrids created from crosses to a common tester, each topcross hybrid has an equal contribution of half of its alleles at all loci from the common tester, which also reduces genetic variation among the hybrids. The reduction of genetic variance, along with the smaller sample sizes, reduced the power of detection in hybrids relative to inbreds.

Candidate genes for Fusarium ear rot resistance

Genetic and biochemical pathways leading to resistance to Fusarium ear rot are entirely unknown. Therefore, GWAS provides a forward genetics approach to screen efficiently many thousands of genes for association with the phenotype without requiring assumptions about what gene functions might be involved in resistance. The SNP associations reported here may help suggest and prioritize candidate genes for resistance to Fusarium ear rot, although we emphasize that associations between genetic variants and phenotypes do not imply that either the SNP is a functional variant or even that the gene containing the SNP is causally involved in resistance. Independent studies, particularly focusing on the biology of the gene functions in relation to infection of maize seeds or other plant tissues, will be required to determine if any of the genes identified here have a role in Fusarium ear rot resistance. Conversely, we expect that GWAS was unable to identify some true functional variants because of the combined effects of small effect size, allele frequency, limited LD in maize, and insufficient SNP density.

The genes containing significant SNP associations in this study include a thioredoxin gene, an F-box gene, a loricrin gene, a sucrose synthase, a CTD1 gene, and a gene of unknown function. The common theme among the likely functions of these genes is that they are very generally important for a variety of cellular processes. The thioredoxin protein family is involved in redox signaling for nearly every plant cellular process [35]; F-box genes are one of the most abundant gene superfamilies in plants and their protein products are involved in uniquitination and degradation [36]; loricrin is likely to be involved in cell membrane function, sucrose synthase is a key enzyme in plant metabolism, and CTD1 is involved in DNA replication. Because of the generality and importance of these gene classes, variation in their function is expected to affect a variety of cellular mechanisms, complicating their possible functional relationship to Fusarium ear rot resistance. Thus, our ignorance of the pathways to resistance to Fusarium suggests that the gene containing a SNP association but no known function should have similar priority for further research as the other candidate genes.

In addition to the genes containing the associated SNPs, there are some cases where LD appears to be sufficiently extensive as to suggest other genes in the region may be important. Around the associations reported between 7.5 and 9.4 Mb on Chromosome 4, for example, it is clear that SNPs in a number of genes across this nearly 2 Mb region are in high LD and will share the association signal with the functional variant in this region. There is a nearby cluster of defense-related and wound-induced proteins around 9.6 to 9.7 Mb that might be considered as putative candidates for further research. Two of those genes had no SNPs in our data set, so we cannot test their associations directly with these data. A few other genes very close to some of significant association also lacked SNPs for testing (Figure 5; Additonal file 1: File S7), and these could not be ruled out as potential candidates. Outside of the region on the short of Chromosome 4, however, the LD decay appears so rapid that it seems unlikely that the SNP associations are more than a few kb from a functional variant. Finally, we also note that there are some larger intergenic regions that lack SNPs (Figure 5), and some sequence variation in these regions may impact gene regulation important to ear rot resistance, but we are likely to miss many such variants in our GWAS scan.

Germplasm and experimental design

In 2010, the NCRPIS collection of inbred lines [26] was evaluated for disease resistance at the Central Crops Research Station in Clayton, NC. The 2010 field experiment consisted of 2572 inbred line entries and was arranged in an augmented single replicate design. Experimental entries were divided into 18 sets of differing sizes based on maturity and field assignment, and sets were then randomly subdivided into incomplete blocks (where the maximum block size across sets was 23 plots). Each block within each set was augmented with a B73 check plot in a random position, and five other checks of varying maturities (IL14H, Ki11, P39, SA24, and Tx303) were included once per set in a random position.

In 2011 and 2012, a novel association mapping panel consisting of 771 diverse inbred line entries was evaluated for disease resistance in Clayton, NC. Based on phenotypic information from the 2010 field experiment, a subset of 681 inbred lines from the NCRPIS collection representing a range of both pedigrees and disease severity scores was chosen for the panel. An additional 90 lines, mostly modern public lines available from North Carolina State University as well as a few lines developed by private industry with recently expired Plant Variety Protection Act (exPVPA) coverage that had become available through the NCRPIS in the spring of 2011 were included. The complete panel of 771 entries was divided into eight sets based on maturity and replicated across the two years using an augmented design. Within years, sets were randomized within the field, and each set was blocked using an α-lattice design [38]. Similar to the NCRPIS evaluation, each block was augmented by a randomly assigned B73 check plot, and five other checks representing a range of maturities and disease reactions (GE440, NC358, 794, B47, and Tx303) were included once per set.

Topcross F₁ hybrids representing a subset of inbred lines from the 2011–2012 association panel were also evaluated in Clayton, NC in 2011 and 2012. Due to seed availability, topcross seed was limited to a sample of 405 inbred lines from the total 771 entries of the association panel. F₁ hybrid seed was generated by crossing inbred lines to either the stiff stalk exPVPA inbred tester PHB47 or the non-stiff stalk exPVPA inbred tester PHZ51 (or both). Overall, 92 lines were crossed only to B47, 162 lines were crossed only to PHZ51, and 151 lines were crossed to both testers, resulting in a total of 556 F₁ hybrid entries in the topcross panel. In the 2011 and 2012 field experiments, topcross entries were classified by tester and maturity (early or late, for a total of four tester × maturity combinations), and each tester × maturity combination was randomly subdivided into three groups. One random group of each tester × maturity combination was assigned to a set, for total of three sets (with four groups per set). Similar to the inbred association panel, sets were randomized within the field in each year, groups were randomized within set, and each group was then subdivided into incomplete blocks, but the topcross hybrids were grown in different field blocks than the inbreds. Each block was augmented with a B73 × PHZ51 topcross check plot in a random position, and two other hybrids that exhibited relatively good resistance to Fusarium ear rot in previous experiments (Pioneer 31G66 and NC478 × GE440) were included once per group. Lastly, one additional check plot of P39 × PHZ51 or CML52 × PHZ51 was included once per group depending on maturity (early or late, respectively).

Inoculation and phenotyping methods

The 2010 NCPRIS collection experiment and the 2011/2012 inbred association panel experiments were inoculated with local toxigenic Fusarium verticillioides isolates using the toothpick method [12],[22]. Approximately one week after flowering, a toothpick containing dried F. verticillioides conidia was inserted near the base of the primary ear of five plants in each plot. At maturity, inoculated ears were harvested and visually scored for Fusarium ear rot symptoms. Scores were assigned to each ear in increments of 5% from 0% to 100% diseased based on the percentage of the ear displaying disease symptoms [19].

Topcross hybrid experiments in 2011 and 2012 were inoculated with a suspension of F. verticillioides conidia using the method described by Robertson et al. [19]. Approximately one week after flowering, 5 mL of a liquid suspension containing 2 × 10⁶ conidia mL⁻¹ was injected into the silk channel of the primary ear of five plants in each plot. One week following the first inoculation, 5 mL of the conidia suspension was injected near the base of the primary ear of the same plants inoculated in the first week. At maturity, inoculated ears were harvested and visually scored using the same protocol as the inbred disease experiments. Raw data from both the inbred and topcross experiments are provided in supplemental datasets File S1 and File S2 in Additional file 1, respectively.

Genotypic data

The genotypic data used in this study consisted of 200,978 SNPs filtered from the GBS markers developed by Romay et al. [26]. The original set of markers consisted of 681,257 SNPs generated by the approach described by Elshire et al. [27] and Glaubitz et al. [28] with missing data imputed using the haplotype-based imputation method described by Romay et al. [26]. SNP data are available at http://panzea.org/db/gateway?file_id=Romay_etal_2013_imputed_geno_data. In addition, the Romay et al. [26] marker set was augmented with GBS data for the ninety inbred lines in the 2011/2012 association panel that were not present in the NCPRIS collection in 2010. GBS data for the aforementioned lines were obtained through the Institute for Genomic Diversity at Cornell Unversity, Ithaca, NY (http://www.igd.cornell.edu). Even after haplotype-based imputation, some missing genotypes exist because the imputation method of Romay et al. [26] does not impute missing data when the observed scores within a test haplotype window do not sufficiently match the reference haplotype set. Therefore, the augmented SNP marker set was then filtered to include only those markers that had less than 20% missing data (after haplotype-based imputation) and a minor allele frequency (MAF) greater than 5%. Duplicate samples present in the Romay et al. [26] data set were also removed from the augmented data set; after this filtering step, genotypic data were available for a total of 2480 inbred lines from across all years combined. The final genotypic data set used in the GWAS analyses is provided in supplemental dataset File S3 in Additional file 1.

Statistical analyses

Availability of supporting data

The data sets supporting the results of this article are available at the Panzea.org repository: http://www.panzea.org/db/gateway?file_id=Zila_etal_2014_data_and_supp.