Genome-wide association study reveals a set of genes associated with resistance to the Mediterranean corn borer (Sesamia nonagrioides L.) in a maize diversity panel.

BACKGROUND: Corn borers are the primary maize pest; their feeding on the pith results in stem damage and yield losses. In this study, we performed a genome-wide association study (GWAS) to identify SNPs associated with resistance to Mediterranean corn borer in a maize diversity panel using a set of more than 240,000 SNPs.
RESULTS: Twenty five SNPs were significantly associated with three resistance traits: 10 were significantly associated with tunnel length, 4 with stem damage, and 11 with kernel resistance. Allelic variation at each significant SNP was associated with from 6 to 9% of the phenotypic variance. A set of genes containing or physically close to these SNPs are proposed as candidate genes for borer resistance, supported by their involvement in plant defense-related mechanisms in previously published evidence. The linkage disequilibrium decayed (r(2) < 0.10) rapidly within short distance, suggesting high resolution of GWAS associations.
CONCLUSIONS: Most of the candidate genes found in this study are part of signaling pathways, others act as regulator of expression under biotic stress condition, and a few genes are encoding enzymes with antibiotic effect against insects such as the cystatin1 gene and the defensin proteins. These findings contribute to the understanding the complex relationship between plant-insect interactions.

Background

Corn borers are the primary maize pest in many environments [1,2]. Corn borers feeding on the pith of the stem results in yield losses because stem damage interferes with assimilate movement to developing kernels. They can also attack the ears, promoting secondary fungal infection, leading to contamination of grain with mycotoxins that may affect human and animal health [3,4].

There are different species of borers that attack maize in different parts of the world. The most economically important species are classified into two families: Crambidae and Noctuidae. Within the Crambidae family, the species with economic importance are: Ostrinia nubilalis Hubner in North America, Europe and North Africa; Ostrinia furnacalis Guenée in Asia; Diatraea saccharalis Fabricius from USA to Argentina; Chilo partellus Swinhoe in Southern USA, Central America and the Caribbean; Diatraea lineolata Walker in Central America, the Caribbean region and South America; and Diatraea grandiosella Dyar in North and Central America. Within the Noctuidae family, the main maize borers are: Sesamia nonagrioides Lefebvre in the Mediterranean region, Busseola fusca Fuller in sub-Saharan Africa, and Sesamia calamistis Hampson in West Africa. This study focuses on the noctuid Mediterranean corn borer (MCB) Sesamia nonagroides Lefebvre, the most important insect pest of maize in the Mediterranean region that includes Southern Europe [2,5,6].

The use of transgenic corn which produces Bacillus thutingiensis (Bt) toxins is a good method for controlling these pests, but transgenic crops are not authorized in several European countries under any agricultural system [7] and are not allowed for organic production [8]. In addition, recent studies have reported a reduction of efficacy of Bt transgenes caused by evolved resistance of some important pests [9-11]. The stacking of several resistant genes has been proposed as one means to delay insect adaptation [12]. Natural sources of resistance to stem borers in maize could reveal promising genes for use in either breeding or transgenic approaches to resistance.

In Spain, there are three MCB generations per year and the second and subsequent generations are able of making significant damage on the stem and, secondarily, on the ear. Following artificial infestation, the level of maize host resistance to stem borer is measured by the tunnel length made by the larvae in the stem as well as by a visual scale for kernel damage. These traits have a complex genetic architecture because resistance depends on the plant-insect relationship, which is influenced by environmental conditions and the developmental stage of the host plant [13,14]. The line mean heritability estimates for tunnel length under corn borer infestation ranged among studies from moderate to relatively high (h2 = 0.50 to 0.78) [15-18], depending on the genetic background.

At present, commercial materials with high levels of native resistance to these insects are not available even though breeding for increasing maize resistance to corn borers has been conducted during the last three decades in different regions around the world. Klenke et al. [19] reduced tunnel length by 4 and 6 cm after four cycles of recurrent selection for resistance to the first and second generations of ECB, respectively. Bosque-Pérez et al. [20] described successful results in the development of materials with resistance to S. calamistis and African sugarcane borer (Eldana saccharina Walker). Sandoya et al. [21] achieved a reduction of 1.8 cm per cycle of recurrent selection for tunnel length by MCB. In addition, a negative relationship between resistance to stem borers and yield has been found [22] when selection for improved yield under infestation was practiced [23]. In summary, classical breeding experiments have demonstrated some successful improvement in corn borer resistance, but natural levels of resistance in elite cultivars remain insufficient to manage the pest. Detection of stem borer resistance QTL could enhance breeding for this trait via marker-assisted breeding or genomic selection.

Genetic effects for resistance against borer attack fit an additive – dominant model, although additive effects appear to be the most important in determining resistance to tunnel length and kernel damage [14,24-28]. Therefore, the study of genetic factors involved in maize resistance to borers can be performed using highly inbred lines.

Several studies performed with segregating biparental populations have reported genomic regions containing minor and major quantitative traits loci (QTL) for resistance to stem and leaf attack by European corn borer (ECB; [15-17,29-31]. Fewer studies have reported QTLs for resistance to other borer species such as the Mediterranean corn borer (MCB) [18,32] or some tropical borers [33-35]. In maize, QTLs for resistance to tunneling by corn borers have been detected on all chromosomes, with the most commonly detected regions occurring on chromosomes 1, 2, 3, 5, and 9 [1]. Identification of causal genes underlying the QTL for resistance in these genomic regions could help breeders to transfer the allelic variants that confer resistance from the lines that carry them to elite breeding lines that lack these resistance alleles. The identification of the causal genetic variants or markers in high LD with causal variants in diverse materials would minimize the risk of dragging other genes with negative effect on the agronomic value during the transfer process.

Although conventional QTL mapping based on linkage maps of biparental populations has been an efficient approach to detect regions related to the resistance to corn borers, higher resolution is needed to detect the genes involved in the defense mechanism of the plant. Genome wide association study (GWAS) based on linked disequilibrium (LD) in diverse genetic samples is a relatively new approach which offers higher resolution mapping that under optimal conditions can pinpoint causal genes underlying quantitative trait variation. Exploiting advances in genotyping and sequencing technology, this approach has been successful in detecting genes associated with diseases in humans [36-38], animals [39-42] and different quantitative traits in plants [43-46]. In contrast to the conventional QTL mapping approach based on linkage in a biparental population, GWAS is based on LD among extant lines from different populations, such that a large number of markers covering the whole genome are required [45,47]. In diverse maize samples, LD is low, therefore, many more markers are needed than in autogamous species (with higher LD) to adequately explore the genetic architecture of complex traits [48]. The low LD offers the benefit of better resolution to delineate potential causal genes within small LD blocks.

Many single-nucleotide polymorphism (SNP) markers have been identified and scored on a maize diversity panel (composed of 302 inbred line) that represents the diversity available in public breeding sector around the world [49,50]. The population has been successfully used by the maize community to perform GWAS in economically important quantitative traits such as kernel composition [51], hypersensitive response [52] and Fusarium ear rot resistance [53].

To date, no GWAS for insect resistance in maize has been reported, although a few GWA studies that deal with plant defense mechanisms against insect attack have been reported in other plant species [54,55]. In this study, we performed GWAS to identify SNPs associated with resistance to MCB. GWAS was done in a subset (270 inbreds) of the maize diversity panel using the maize 50 k SNP genotyping array [56] plus a more recently developed set of 425 k SNPs found through genotyping by sequencing [57,58].

Results

Means, analysis of variance and heritabilities

Differences among inbred lines were highly significant (P < 0.01) for all resistance (TL, SD, and KR) and agronomic traits (PH, DTA, and DTS); while significant (P < 0.05) genotype by environment (G × E) interactions were observed for all traits as shown in Additional file 1: Table S1 and S2. The inbred means for TL ranged from 5.2 to 49.2 cm with an overall mean of 20.9 cm, and for SD from 4.1 to 22.6% with an overall mean of 11.9%, respectively. The inbred mean for KR ranged from 5.4 to 9 with an overall mean of 7.8 in the subjective scale. Higher values for TL were observed in 2012 (overall mean = 43.3 cm) compared with 2010 (overall mean = 17.5 cm) and 2011 (overall mean = 14.7 cm, Additional file 1: Table S3 and S4). Intermediate values for heritability on a line mean-basis were estimated for TL (h2 = 0.60) and KR (h2 = 0.52) across the three years, whereas a low heritability value was obtained for SD (h2 = 0.25, Additional file 1: Table S1).

Correlation analysis

A significant and high (r > 0.50) phenotypic correlation was observed between TL and SD. Genetic correlation coefficients were also significant and high between TL and SD and between TL and PH (Table 1). KR showed significant phenotypic correlations with other resistance traits but the correlation coefficients were not higher than 0.5. Genetic correlations between KR and the agronomic traits (DTA, DTS and PH) were significant and high (Table 1).

Table 1

Genotypic (above diagonal) and phenotypic (below diagonal) correlation coefficient estimates for each pair of traits

	TL	SD	KR	DTA	DTS	PH
TL		0.97*	NS	NS	NS	0.51*
SD	0.74*		NS	NS	NS	NS
KR	−0.31*	−0.31*		0.75*	0.71*	0.51*
DTA	NS	NS	0.30*		0.99*	0.66*
DTS	NS	NS	0.28*	0.96*		0.65*
PH	0.24*	−0.06	0.24*	0.51*	0.48*

TL, tunnel length; SD, stem damage; KR, kernel resistance; PH, plant height; DTA, days to anthesis; DTS, days to silking.

1 Genotypic correlation coefficients were considered significant, *, when exceeded twice its standard error; NS, not significant.

2 Phenotypic correlation coefficients were considered significant, *, at 0.05 probability level according to Steel and Torrie [59]; NS, not significant.

Association analysis of maize resistance to MCB and agronomic traits

The compressed mixed linear model computed for each trait in Tassel reduced the pairwise kinship matrix by clustering the 267 lines into 32 groups for TL, 201 groups for SD, 166 groups for KR, 192 for PH, and 240 for DTA, and DTS (Table 2). The proportion of the total phenotypic variation explained by background genetic effects was 39, 13, and 48% for TL, SD, and KR, respectively. By comparison, the background polygenic effects modeled by the K matrix accounted for 60, 94 and 90% of total variation for PH, DTA, and DTS, respectively.

Table 2

Summary of the compressed mixed linear model analysis for three traits related to resistance to MCB attack and three agronomic traits in an inbred association panel evaluated in three years

Trait	n a	s b	Compression Level ( c ) c	σ^g2 d	σ^^2 e	σ^g2σ^g2+σ^^2 f
TL	267	32	8.34	32.80	50.61	0.39
SD	267	201	1.33	1.45	9.51	0.13
KR	265	166	1.60	0.15	0.16	0.48
PH	266	192	1.39	384.70	256.56	0.60
DTA	266	240	1.11	45.03	5.03	0.94
DTS	266	240	1.11	51.19	5.87	0.90

TL, tunnel length; SD, stem damage; KR, borer kernel resistance; PH, plant height; DTA, days to anthesis; DTS, days to silking.

a Total number of inbred lines included in the analysis.

b Number of groups obtained using a clustering approach based on K matrix with the optimum compression level option.

c Compression level is the average number of inbred lines per group estimated as n/s.

d Additive background genetic variance component estimated in Tassel by fitting the K matrix in the MLM without any SNP marker effects.

e Residual genotypic variance component estimated in Tassel.

f Proportion of phenotypic variance explained by the K matrix, estimated as background genetic variance divided by total phenotypic variance.

Ten SNPs were identified as significantly associated with length of tunnels made by MCB (Figure 1, Table 3). Based on the additive effects, the major allele reduce TL for all significant SNPs except for the SNP located on chromosome 10 (Table 3). Four SNPs were significantly associated with SD made by MCB (Figure 1). The minor allele for those SNPs increases the SD from 1.1 to 2%. The total variance explained (R) by each SNP associated with TL and SD ranged from 7 to 9%. Eleven SNPs were significantly associated with KR (Figure 1) and the total variance (R) explained by each of them ranged from 6 to 8%. The minor allele for those SNP reduced from 0.15 to 0.40 points the ratio that accounts for kernel resistance.

Figure 1

GWAS results for the three resistance traits to MCB attack in a maize association panel. Each graph (TL, SD, and KR) represent the P-values of the 246,477 SNPs tested for each resistance trait. Each row indicates the SNP significantly associated (RMIP ≥ 0.30) to each resistance trait analyzed.

Table 3

SNP identification (SNP ID), additive effect and allelic variants for the SNP, proportion of total variance explained by the SNPs significantly associated with resistance traits (TL, SD, and KR), and significance values for the association between the SNP and the phenotype ( -value and RMIP)

Trait a	SNP ID b	Alleles c	(N o ) d	Additive effect e	P -value	( R 2 ) f	RMIP g
TL	S2_168004182	G/T	225/36	3.09	9.51 × 10−6	0.08	0.51
TL	S3_7081859	T/A	207/55	2.52	1.38 × 10−5	0.07	0.50
TL	S4_190444179	G/A	218/38	3.18	3.26 × 10−6	0.09	0.37
TL	ss4_190679094h	G/A	226/41	3.09	3.37 × 10−6	0.08	0.45
TL	S7_154739818	G/C	248/14	4.25	5.28 × 10−5	0.07	0.37
TL	S7_154741622	C/T	217/47	2.84	9.84 × 10−6	0.08	0.56
TL	S7_155702328	T/A	171/67	2.43	1.55 × 10−5	0.08	0.38
TL	S10_133337924	G/T	133/99	−1.30	7.40 × 10−5	0.07	0.30
TL	S10_133337925	G/C	133/99	−1.30	7.40 × 10−5	0.07	0.30
TL	S10_133337950	T/C	133/99	−1.30	7.40 × 10−5	0.07	0.30
SD	S1_208315891	C/A	126/106	1.09	7.86 × 10−6	0.09	0.63
SD	S1_293163491	T/A	207/37	1.37	2.39 × 10−5	0.08	0.38
SD	S2_59729532	A/T	235/17	1.90	2.93 × 10−5	0.07	0.32
SD	ss5_176870721h	G/A	228/38	1.37	1.12 × 10−5	0.08	0.41
KR	S3_187742562	C/A	244/12	−0.38	1.30 × 10−5	0.08	0.48
KR	S3_204458505	T/A	186/65	−0.19	4.14 × 10−5	0.07	0.40
KR	S3_204586960	C/G	251/13	−0.37	1.25 × 10−5	0.08	0.44
KR	S3_222733400	C/T	142/110	−0.15	4.90 × 10−5	0.07	0.34
KR	S4_227101950	C/T	242/17	−0.30	1.26 × 10−5	0.08	0.43
KR	S4_227101985	A/T	242/17	−0.30	1.26 × 10−5	0.08	0.43
KR	S5_93580059	C/T	230/28	−0.24	2.27 × 10−5	0.07	0.31
KR	S6_88149024	A/G	187/55	−0.17	6.83 × 10−5	0.07	0.48
KR	S6_88149036	G/C	187/55	−0.17	6.83 × 10−5	0.07	0.48
KR	S7_15072370	G/A	230/24	−0.26	3.65 × 10−5	0.07	0.46
KR	S7_19347596	A/G	250/11	−0.40	7.47 × 10−5	0.06	0.36

a TL, tunnel length in cm; SD, stem damage in percentage; and KR, kernel resistance on a subjective visual scale of 1 to 9 in which 1 indicates completely damaged and 9 indicates no damage.

b The number before the underscore (_) indicates the chromosome number and the number after the underscore (_) indicates the physical position in bp within the chromosome.

c The letter before the diagonal (/) is the nucleotide more frequent ; and the letter after the diagonal the nucleotide less frequent.

d N° = number of inbred lines homozygous for a determined allelic variant. The number before the diagonal (/) represents the number of individuals with the mayor allele; and the number after the diagonal represents the number of individuals with the minor allele.

e The additive effect was calculated as half the difference between the mean of the homozygous for the minor and the mean of the homozygous for the major allele.

f R , proportion of the phenotypic variance explained by the SNP.

g RMIP, resample model inclusion probability.

h Based on SNPs from Illumina chip, the remaining locations without a superscript are based on SNPs obtained by GBS.

Nineteen SNPs were significantly associated with PH and 48 and 43 SNPs were significantly associated with DTA and DTS, respectively (Additional file 1: Table S5, S6, and S7). But none of those SNPs coincide or are close to those detected for resistance traits.

Candidate genes selection

The filtered predicted gene set from the annotated B73 reference maize genome [60] was used to characterize the gene containing or nearby the SNP declared significant. Seven candidate genes containing or adjacent to the SNPs significantly associated with TL, four candidate genes containing or adjacent to the significant SNPs associated to SD, and ten candidate genes containing or adjacent to the significant SNPs associated to KR were proposed (Table 4).

Table 4

Candidate genes for each significantly SNP associated with TL, SD, and KR and its respective encoding product

Chr a	Trait	SNP physical position (bp)	Gene ID	Encoding
2	TL	168,004,182	GRMZM2G504910	Tetratricopeptide repeat (TPR)
3	TL	7,081,859	GRMZM2G104081b	hex1 (hexokinase1)
4	TL	190,444,179	GRMZM2G013128b	Double Clp-N motif-containing P-loop nucleoside triphosphate hydrolases superfamily protein
4	TL	190,679,094	GRMZM2G033820	Phospholipase A2
7	TL	154,739,818	GRMZM2G077008bc	Histidine kinase, hybrid-type, ethylene sensor
		154,741,622
7	TL	155,702,328	GRMZM2G861541	Expressed protein
10	TL	133,337,924	GRMZM2G057084	Calcium/calmodulin-dependent protein kinase
		133,337,925
		133,337,950
1	SD	208,315,891	GRMZM2G101422	Expressed protein
1	SD	293,163,491	GRMZM2G060702b	Actin depolymerizing factor 4
2	SD	59,729,532	GRMZM2G389097	Leucine-rich receptor-like protein kinase family
5	SD	176,870,721	GRMZM2G325683c	Expressed protein
3	KR	187,742,562	GRMZM2G438551	cystatin1
3	KR	204,458,505	GRMZM2G111666	basic Helix-Loop-Helix (bHLH) transcription f.
3	KR	204,586,960	GRMZM2G091494c	Starch branching enzyme interacting protein-1
3	KR	222,733,400	GRMZM2G055578c	Glycine-rich protein
			GRMZM2G055629	Plant thionin family protein precursor/Defensin
4	KR	227,101,950	GRMZM2G116314	Ubiquitin thiolesterase
		227,101,985
5	KR	93,580,059	GRMZM2G037308c	Phytosulfokine receptor
6	KR	88,149,024	GRMZM5G876960	Polyamine oxidase (propa-1,3-diamine-forming)
		88,149,036
7	KR	15,072,370	GRMZM2G316256	Catalase//L-ascorbate peroxidase
7	KR	19,347,596	GRMZM2G042627	Kinase associated protein phosphatase

a Chr, chromosome; TL, tunnel length; SD, stem damage; and KR, kernel resistance.

b Gene containing the significant SNP within an exonic region.

c Gene containing the significant SNP within an intronic region.

In general, the LD (r) between significantly associated SNPs and SNPs around them decays (r ≤ 0.10) rapidly; but, in some cases, the LD spans as much as 0.5 Mb at an r value greater than 0.20 (Figure 2). For each trait, LD estimates between significant SNPs were always below 0.1 except for SNPs significant for TL located on chromosome 7 as shown in Additional file 1: Table S8.

Figure 2

Linkage disequilibrium heat chart showing LD measure ( ) between the SNP significantly associated with traits related to resistance to MCB attack and the closest 60 SNPs. Each bar represent a region (ranging from ~0.15 to ~1 Mbp) containing each significant SNP associated (black square) to resistance traits. The LD (value of r ) with the 30 upstream SNPs were shown at the right side of the black square, and the LD with the 30 downstream SNPs were shown at the left side of the black square; on each bar, the extreme distances (in kbp) covered by the upstream and downstream SNPs are indicated.

Discussion

Means, heritabilities, and correlations

Substantial differences among yearly means for TL made by MCB were observed across environments, with the highest values of TL observed in 2012. Artificial infestation is used to ensure the contact of the insect with the plant, but the severity of the attack is conditioned by environmental factors that in turn influence natural infestation. The prevalence of the pest can vary greatly between and within locations [5,6]. Data on the monthly samplings at different locations in Pontevedra (unpublished data) indicate a higher rate of natural infestation in 2012 compared with 2010 and 2011 in a plot adjacent to the GWAS trial. In addition, the experiment was harvested earlier in 2011 compared with 2010 and 2012, limiting the time that larvae could damage the plants. Therefore, differences on infestation levels and harvest time could be major causes of the observed differences in TL among years.

Significant genotype × environment interaction was observed for all traits, although no G × E interaction for TL was found in previous studies under infestation with MCB [18,61]. However, G × E interaction for resistance traits was not significant, except for KR, when we made the analysis discarding data from 2012. Therefore, the high values for TL in 2012 could be obscuring the genotype effect. In previous studies, typically no G × E interaction for KR was found [61-63], with some exceptions [18].

The heritabilities for resistance traits ranged from low to moderate while PH and flowering time were higher as expected. The heritability for TL estimated herein is within the range (h2 = 0.40-0.77) of those obtained in other works with biparental crosses under infestation with MCB or ECB [16,18,64].

A significant genetic relationship was observed between TL and PH (r = 0.51) in a panel of diverse origin. Nevertheless, we did not find a single SNP or linked SNPs associated with both traits, suggesting that it is primarily due to the polygenic background. This relationship has also been important in some previous experiments with MCB and ECB [16,32], but negligible in others [17,30].

A high genetic correlation coefficients between KR and flowering time (r = 0.75 for DTA and 0.71 for DTS) and between KR and PH were found, which suggest that late and taller genotypes will be the healthiest; but these results have to be taken with caution because infestation was made simultaneously for all genotypes placing the MCB eggs close to the ground and the higher distance between the eggs and the ear of the taller and later plants could have impeded the larval arrival to the ear. In addition, when stem tissue is more lignified at time of infestation (as would be the case for earlier flowering genotypes), the preference of MCB larvae for the stem tissue compared to ear tissue could be less evident.

Association analysis

There was a minimal variation in model fit of the compressed MLM among different traits because similar compression level values were observed for all traits (c = 1.1 to 1.6), except for TL, which had the highest value for compression level (c = 8.3). However, Zhang et al. [65] shown that this method controls the false positive rate well when the compression levels ranged from 1.5 to 10.

SNPs located on chromosomes 3 and 7 for TL co-localized with previously reported QTLs for TL by corn borers in genome bins 3.02 and 7.03 [30]. The proportion of the phenotypic variance explained (R2 = 7 - 9%) by each SNP was comparable to that explained by QTLs for TL made by ECB and MCB (R2 = 3.5 – 15.7%) in biparental crosses [16-18,64]. No QTLs for SD and KR have been previously reported in biparental crosses in the same regions where we located the significant SNPs, except one QTL for KR made by MCB at the bin 5.04 [66]. Therefore, association mapping uncovers additional genomic regions involved in maize resistance to corn borers that were not detected using biparental populations.

Candidate genes

We used the maize B73 genome v2 (RefGen_v2) available from the Maize GDB [67] (http://www.maizegdb.org/) to identify genes that either include or are close to the significantly associated SNPs. A region of approximately 0.2 Mb around the SNP was checked for annotated genes putatively involved in plant response to wounding and/or damage by microbes (insects or pathogens) based on bibliographic records.

Genes associated with TL made by MCB

The candidate gene adjacent to the significant SNP associated to TL on chromosome 2 encodes a Tetratricopeptide repeat (TPR) protein containing. The TPR is one of many repeat motifs that form structural domain mediating protein-protein interactions in several cellular process including translocation and degradations of proteins [68,69]. A recent study has proposed that the presence of those protein-protein interaction motifs could be acting as a modulator of the gene function and protein expression during the stress response caused by invading pathogens [70]. The hex1 gene containing the significant SNP on chromosome 3 encodes hexokinase, a sugar sensor with numerous physiological functions within the cell including response to oxidative-stress and pathogen resistance [71-73]. A gene located on chromosome 4 that encodes a Double Clp-N motif-containing P-loop nucleoside triphosphate hydrolase superfamily protein contains a significant SNP. This gene shows a weak similarity to the AtHSP101 gen in Arabidopsis, that codifies for a heat shock protein required for acclimation to high temperatures, and probably could be involved in response to other stresses [74]. Furthermore the ofp44 (OVATE-transcription factor 22) gene is relatively close to this SNP. It is known that some transcription factors from the Ovate family protein interact with other transcription factor families such as NAC domain protein1, MYB transcription factors, and KNOX homeodomain protein to regulate the synthesis of the three major components of secondary cell wall (lignin, cellulose and hemicellulose) in A. thaliana [75-78] and plants with a fortified cell wall are more resistance to corn borer attack [79].

The candidate gene for the other significant SNP on chromosome 4 encodes a Phospholipase A2 (PLA2) which plays a very important role in signal transduction in plants since it is the precursor of oxylipins and jasmonated acid, two hormones which regulates defense genes against herbivores [80-82].

A 1 Mbp region on chromosome 7 contains three SNPs that were significantly associated with TL and were in significant LD with at least one of the other associated SNPs in the region (r2 > 0.2; Additional file 1: Table S8). The SNPs at 154,739,818 and 154,741,622 bp are both located within a gene that putatively encodes for a Histidine kinase, hybrid-type, ethylene sensor; five other genes encoding Serine/Threonine kinase receptor and receptor-like ser/thre kinases family proteins (RLK) were also physically nearby. It is well known that both kinases and RLK proteins play a central role in signaling during pathogen recognition and the subsequent activation of plant defense mechanisms [83-85]. They are also involved in wound-mediated defense response [86], and maintenance of plant cell wall integrity [87]. Polymorphisms at the proposed kinase genes could also be responsible for the QTLs at bin 7.03 detected for TL in a biparental population [30]. The candidate gene encoding the maize Calcium/calmodulin-dependent protein kinase (CDPK) close to the three significant SNPs on chromosome 10, could be playing an important role in the activation of defense against the attack of MCB since it has been known that the CDPK is induced by mechanical wounding by herbivore attack inducing accumulation of jasmonic acid in maize and other species [88-91].

Genes associated with SD made by MCB

A gene encoding an Actin polymerizing factor 4 (APF4) contains the second SNP at chromosome 1 significantly associated to SD. One of the functions of APF4 protein is remodeling the actin of cytoskeleton under different stimulus, including wounding and pathogen attacks [92]. Some studies in Arabidopsis indicated that the APF4 mediated defense signal and it is also relates with actin dynamic of cytoskeleton during the innate immune signaling [93-95]. The candidate gene for the SNP significantly associated to SD on chromosome 2 encodes a Leucine-rich receptor-like protein kinase (LRR-RLK). LRR-RLK protein family plays an important role in cell-cell signaling and other signals involving peptide in ligands. They are involved in systemic activation of protease inhibitors in response to wounding by insect feeding [96,97]. On the other hand, the significant SNP located on chromosome 5 is within a gene encoding a protein with unknown function, and it is interesting that the SNP is close to the gene nactf30 which encodes a NAC domain protein transcription factor. As already mentioned, this gene and other transcription factor family members regulate the synthesis of secondary cell wall [98,99]. In addition, other authors have associated the NAC domain proteins with response to stresses made by herbivore attack [100].

Genes associated with KR made by MCB

The SNP on chromosome 3 significantly associated to KR was located nearby the cystatin1 gene, whose product is the corn kernel cysteine proteinase inhibitor//cysteine proteinase inhibitor I (psei1), an anti-metabolic protein synthetized and stored in the maize kernel [101]. The expression of some proteinase inhibitor genes are induced in response to mechanical wounding and insect damage [102], and it has been shown that cysteine proteinase inhibitor interferes with the digestive process of insects [103,104]. LD is low in this region, but it is interesting that the SNP significantly associated with KR was in LD with a SNP positioned in the exonic region of the gene (Figure 3). Therefore, if the association found between the SNP at position 187,742,562 and KR is due to the linkage between that SNP and a certain undetected polymorphism in the cysteine1 gene, the effect of this polymorphism would be expected to be especially large because the linkage between the significant SNP and SNPs in the cystatin1 gene is low, although significant.

Figure 3

A region of approximately 280 kpb in chromosome 3 where a SNP significantly associated with borer kernel resistance (KR) has been found. The red points represent the P-values of the SNPs in the region mentioned above. The orange strip indicates the location of the SNP significantly associated with KR and the dark blue strips indicate genomic locations in LD (r  > 0.2) with the SNP significantly associated with KR. The solid black line and the green rectangle indicate the intronic and exonic region of the cystatin1 gene, respectively. The SNP significantly associated with KR herein is in LD with the exonic region of the cystatin1 gene, which encodes a maize cysteine proteinase inhibitor involved in plant resistance to insects.

Another significant SNP associated to KR was located close to a gene encoding a basic Helix-Loop-Helix (bHLH) transcription factor family member. Recent studies in other plant species have demonstrated that this family gene protein has an important role in regulating jasmonic acid response [105-109]. The third significant SNP associated with KR located on chromosome 3 was within a gene encoding a starch branching enzyme interacting protein, whose deficiency leads to decreased digestibility of maize kernel [110]. No previous evidence about its involvement in resistance to insect was found, however. A gene that encodes a Germin 1–2 protein is 106 kb downstream from this associated SNP; these proteins are implicated in the response to abiotic and biotic stresses including the response to mechanical wounding and insect damage [111-114]. Another significant SNP located on chromosome 3 at position 222,733,400 is within a gene that encodes a Glycine rich protein (GRP), a structural protein which could be playing an important role in the fortification of plant cell wall [115-117], also this protein is wound-inducible [118,119]. A set of plant Thionin family protein precursor genes were found nearby the SNP on chromosome 3 at position 222,733,400. This finding is particularly interesting since these proteins belong to the Defensin family protein which has antibacterial, antifungal, and insecticide activities [120-124].

The two significant SNPs associated to KR on chromosome 4 are adjacent to a gene that encodes for an Ubiquitin thiolesterase. These proteins form complex systems of selective protein degradation [125,126] and mediate the biosynthesis of plant hormone signaling such as salicylic, jasmonic, and abscisic acid and auxins [127,128].

The candidate gene containing the significant SNP associated to KR on chromosome 5 encodes a Phytosulfokine receptor (PSK), a recognition hormone whose level of expression has been increased by pathogens elicitors [129,130].

The candidate gene, adjacent to the significant SNP located on chromosome 6, encodes a Polyamine oxidase (propa-1,3-aimine-forming) (PAO). PAO plays an important role in stress tolerance by generating H2O2 which is a key component in signal transduction pathways leading to stress responses. In Zea mays, the function of PAO in wound-healing is likely due to increased lignin and suberine deposition as consequence of H2O2 release [131]. In addition, it has been described to play a role in cell wall stiffening and mediating abiotic and biotic stresses [132]. The candidate adjacent gene to the significant SNP associated to KR on chromosome 7 encodes a Catalase//L-ascorbate peroxidase, they are two major hydrogen peroxide-detoxifying enzymes whose activity is very important in the reduction of the oxidative stress caused by H2O2 [133,134], and the importance of these detoxifying enzymes in resistance to insect attack has been recently reported [135].

A candidate gene for the significant SNP on chromosome 7 at position 19,347,596 is adjacent to the SNP that encodes a Kinase associated protein phosphatase. This gene family has been proposed to regulate the response to different type of stresses including pathogens and herbivore attack [136]. It is interesting to note that another close gene to this SNP is the ipt4 (isopentenyl transferase4) gene, which is involved in the regulation of cytokinin biosynthesis pathway [137]. The expression of an itp gen (fused with a wound-inducible promoter) in transgenic plants of Nicotiana plumbaginifolia decreased leaf consumption by Manduca sexta (lepidopterous) and reduced survival of Myzus persicae (aphid) [138]. Although there are several reports about the importance of cytokinins in the modulation of plant defense against pathogens and insect attack, their role is not clear [139].

Unlike markers linked to QTLs for resistance to corn borers, the SNPs associated with resistance in the present study could be incorporated in whole-genome predictor models in order to improve genomic selection [140]. Marker–assisted selection has proved an useful tool for improving resistance to the European corn borer [141] of segregating bi-parental populations related to the mapping populations used for QTL detection, but could be inappropriate in non-structured populations. As additive effects are the most important genetic effects for resistance traits, crossing inbreds with improved resistance will render hybrids more resistant to attack by MCB larvae.

Conclusion

We conducted a genome wide association study for resistance traits to MCB with more than 245 kb SNP distributed through the whole genome. We found a set of significant SNPs associated to the three resistance traits to MCB attack. Of which 10 SNPs were significant associated to TL, 4 SNPs were associated to SD, and 11 SNPs were associated to KR. In general, each of these SNPs explain a considerable proportion of the phenotypic variance (R = 6-9%). No co-localized SNPs were found for resistance and agronomic traits that could underlie the genetic correlations found between these traits.

Twenty one candidate genes were proposed for the three resistant traits, they were either containing or adjacent (within a window of ±130 kbp) to each significantly associated SNP.

Most of the candidate genes proposed herein are part of the signaling pathway, others act as regulator of expression under biotic stress condition, and a few genes are encoding enzymes with antibiotic effect against insects such as the cystatin1 gen and the defensin proteins.

The identification of these polymorphisms associated to resistance traits to MCB attack can be useful to understand the molecular mechanisms that affect resistance and susceptibility of host plants to insect attack, in order to contribute to advance in the understanding of plant-insect interactions. Nevertheless further studies are necessary to validate the candidate genes identified herein.

Methods

Plant material and phenotypic data

The maize diversity panel (composed of 302 inbred lines) represents much of the diversity available in public breeding sector around the world. A subset of the maize diversity panel (henceforth we will refer to this population as “association panel”) was evaluated for resistance to MCB attack at Pontevedra (42°24’ N, 8°38’ W, and 20 m above sea level), Spain, through three years (2010, 2011, and 2012). A subset of 270 inbred lines was assayed in an 18 × 15 α-lattice design with two replications in 2010 and 2011. In the third year a subset of 255 inbred lines was assayed (because we did not have enough seed for the remaining15 lines) in a 17 × 15 α-lattice design with two replications.

The trials were hand-planted and each experimental plot consisted of one row spaced 0.8 m apart from the other row with 29 two-kernel hills spaced 0.18 m apart. Plots were overplanted and thinned, obtaining a final density of ~70,000 plant ha−1. The evaluations were performed under artificial infestation with eggs of MCB. The eggs for inoculation were obtained at the Misión Biológica de Galicia by rearing the insect as described by Eizaguirre and Albajes [142]. Five plants of each plot were infested with ~ 40 MCB eggs placed between the stem and the sheath of a basal leaf.

Data collected were: tunnel length (TL), the mean length of stem tunnels made by borers on the five infested plants, which were longitudinally split at the time of harvest; stem damage (SD) as the percentage of the stem damaged by MCB larvae on the five infested plants; kernel resistance to borer attack (KR) recorded at harvest as the damage on the main ear of the five infested plants according to a subjective visual resistance scale of 1 to 9 in which 1 indicates completely damaged and 9 indicates no damage; days to anthesis (DTA) and days to silking (DTS) as the days from planting to the date on which 50% of plants were shedding pollen or showing silks, respectively; and plant height (PH) on five representative plants as the distance from the ground to the top of the plant.

Genotypic data

We used a set of unique SNP markers derived from a Illumina maize 50 k array [56] and a genotyping-by-sequencing (GBS) strategy [48]. The two data sets were combined and filtered to exclude SNPs with more than 20% missing genotype data and minor allele frequency (MAF) less than 5%. Heterozygous genotypes were considered as missing data in the analysis. After filtering, a total of 246,477 SNPs (Additional file 2) distributed across the maize genome were used in this study.

A genetic kinship matrix (K) previously published by Olukolu et al. [52] was used for GWAS. The kinship matrix was estimated using a subset of 5000 SNPs without any missing genotypes and distributed approximately uniformly across the entire genome (at least 60 kbp between any two markers).

Statistical analyses

Best linear unbiased estimator (BLUE)

Each trial was analyzed separately with the SAS mixed model procedure (PROC MIXED) in SAS software version 9.3 [143] considering inbred lines as a fixed effect and replications and block within replication as random effects. Then, trials were combined using a mixed linear model across years and considering inbred lines as the only fixed effect. As large predicted values for stem damage and tunnel length were associated with larger residuals, a natural logarithmic transformation of TL and SD scores was used for obtaining BLUEs. The logarithmic transformation eliminated the relationship between residual variance and predicted values. Line BLUEs were back-transformed and then used to perform GWAS.

Heritabilities

Heritabilities (ĥ2) for each year were estimated for each trait on a family-mean basis as described previously by Holland et al. [144]. The model for these analyses was similar to the model mentioned above with the exception that inbred lines were treated as random effects. The genetic (r) and phenotypic (r) correlations between traits were computed using REML estimation in SAS mixed model procedure following Holland [145].

Genome-wide association analysis based on mixed linear model (MLM) was performed in Tassel 4.1.26 [146]. The MLM used by Tassel waswhere y is the vector of phenotypes (BLUEs), β is a vector of fixed effects, including the SNP marker tested, u is a vector of random additive effects (inbred lines), X and Z represents matrices, and e is a vector of random residuals. The variance of random line effects was modeled as , where K is the n × n matrix of pairwise kinship coefficient and is the estimated additive genetic variance [147].

Restricted maximum likelihood estimates of variance components were obtained by using the optimum compression level (compressed MLM) and population parameters previously determined options (P3D) in Tassel [65]. The optimum compression level option reduces the computation demand by clustering the (n) total individuals into (s) groups based on their realized genomic relationships, allowing the original K matrix to be replaced by a smaller relationship matrix. The P3D option uses iteration to estimate population parameters such as genetic and residual variance only once in a model with no fixed marker effects, then uses those estimates without iteration in subsequent association tests for each marker. The combination of these two methods reduces computational time and improves model fit [65].

Threshold for GWAS

To identify SNPs with the most robust associations with traits, a subsampling or subagging procedure was employed in GWAS analysis [148,149]. Each of 100 subsampled datasets generated using the R software [150] comprised a random sample of 80% of inbred lines from the diversity population. Only SNP markers determined as significant at P < 1 × 10−4 and subsequently detected in ≥ 30 subsamples, i.e. resample model inclusion probability (RMIP) threshold of 0.30, were considered as significantly associated to the trait under study.

Linkage disequilibrium and candidate gene selection

We examined the linkage disequilibrium (LD) measure (r2) with each SNP significantly associated with resistance traits in a region of 60 SNPs (ranging from 0.15 to 1 Mbp). For each trait, the linkage disequilibrium between significant SNPs was also calculated. The genes containing or adjacent to SNPs significantly associated with traits were identified and characterized by the use of MaizeGDB genome browser [67]. We examined a region around each significant SNP in order to identify candidate genes of interest. For most SNPs, more than two genes that could be involved in plant defense mechanism (based on previously published evidence) were preselected; the closest gene to the SNPs significantly associated with resistance was selected as candidate gene and the remaining genes were presented as other interesting genes that could be putatively involved in plant defense against insect herbivores.