Identification of quantitative trait loci for related traits of stalk lodging resistance using genome-wide association studies in maize (Zea mays L.)

Abstract

Background

Stalk lodging is one of the main factors affecting maize (Zea mays L.) yield and limiting mechanized harvesting. Developing maize varieties with high stalk lodging resistance requires exploring the genetic basis of lodging resistance-associated agronomic traits. Stalk strength is an important indicator to evaluate maize lodging and can be evaluated by measuring stalk rind penetrometer resistance (RPR) and stalk buckling strength (SBS). Along with morphological traits of the stalk for the third internodes length (TIL), fourth internode length (FIL), third internode diameter (TID), and the fourth internode diameter (FID) traits are associated with stalk lodging resistance.

Results

In this study, a natural population containing 248 diverse maize inbred lines genotyped with 83,057 single nucleotide polymorphism (SNP) markers was used for genome-wide association study (GWAS) for six stalk lodging resistance-related traits. The heritability of all traits ranged from 0.59 to 0.72 in the association mapping panel. A total of 85 significant SNPs were identified for the association mapping panel using best linear unbiased prediction (BLUP) values of all traits. Additionally, five candidate genes were associated with stalk strength traits, which were either directly or indirectly associated with cell wall components.

Conclusions

These findings contribute to our understanding of the genetic basis of maize stalk lodging and provide valuable theoretical guidance for lodging resistance in maize breeding in the future.

Background

Maize (Zea mays L.) plays an important role in food security, feed provision, and fuel resources. Nevertheless, stalk lodging can lead to 5–20% maize yield loss annually worldwide [1]. Achieving high agricultural yields under different environmental conditions is a major goal of maize breeders. In low-density populations, the yield was improved by selecting taller plants to increase the biomass per plant. In high-density populations, the high yield was obtained by increasing the population density of selected medium height plants through the combination of reasonable panicle height coefficient and lodging resistance. Stable quantitative trait loci (QTLs) are particularly useful in marker-assisted selection [2]. Stalk lodging is a phenomenon whereby plants collapse from the upright state, a complicated and integrated quantitative trait caused by many factors, such as the quality of the stalk itself and the external environmental factors (e.g., climatic and soil conditions, planting density, fertilization and irrigation, pests and diseases) which cause irreversible damage to corn stalks and roots [1, 3]. Maize lodging can be divided into three types: root lodging, stem bending, and stem breaking [4]. Stalk lodging usually occurs at or below the ear node, which consequently influences the regular growth of the ear before harvest and the final yield of maize [5, 6]. Furthermore, grain yield per unit area is highly correlated to the plant’s adaptability to high crop density, but stalk lodging limits planting density and mechanized harvesting [7, 8]. Therefore, improving stalk lodging resistance in maize would benefit future breeding programs and agricultural production.

Stalk lodging resistance is correlated with stalk mechanical strength, hence this variable was used to evaluate lodging resistance in maize [9, 10]. Common methods to quantify the stalk mechanical strength include rind penetration, bending, breaking, and vertical crushing [4, 7, 11]. Most studies have found that the stalk rind penetrometer resistance (RPR) and stalk buckling strength (SBS) are important determinants of crop lodging resistance. Furthermore, RPR did not damage the stalk structure [12,13,14]. Compared with RPR, SBS is more closely correlated to stalk lodging under natural conditions, as stalk lodging happens in case of over-bending [15]. According to previous studies, we found that lodging occurs most frequently at flowering stage or a few weeks after flowering and the third or fourth internode of maize plants is extremely sensitive to stalk lodging in the field [6, 8, 13, 16]. Furthermore, Liu et al. [11] showed that the best period for evaluating stalk strength is the silking phase or stage after silking. The position of the stem lodging mainly occurs between the second and fifth internodes, especially in the third internodes and the fourth internodes above ground (FIAG) were significantly correlated with RPR and SBS [6, 8, 11, 17, 18]. In addition, with the increase of plant density, the length of the base nodes increased significantly, the diameter of the stems decreased significantly, and the content of cellulose, hemicellulose and lignin decreased, resulting in a decrease in the mechanical strength of the stems and an increased risk of lodging [19].

QTL mapping has been widely used in the study of various agronomic traits, including yield-related traits, which is a useful tool for analyzing the genetic structure of complex agronomic traits. In crop, QTL mapping on lodging have been gradually applied in sorghum, wheat, rice, especially in maize. For example, a linkage map with 129 SSRs markers was constructed by Hu et al. [6], and two, three, and two QTLs were detected for the maximum load exerted to breaking (F max), the breaking moment (M max) and the critical stress (σ max), respectively. Li et al. [12] identified seven QTLs associated with RPR in two maize recombinant inbred line (RIL) populations using 3072 single nucleotide polymorphisms (SNP) markers. Zhang et al. [17] identified 44 significant QTLs for SD, SBS, and RPR using the IBM Syn10 DH population in three environments.

The efficiency and accuracy of QTL mapping depend largely on the marker density, the variation range of phenotypes within the population, as well as the population size and type [20]. Genome-wide association study (GWAS) is a powerful tool for analyzing the genetic basis of complex traits. So far, GWAS has been used to analyze many agronomic traits such as plant height, leaf structure and yield-related traits [21,22,23], and other characteristics, i.e. In addition, some genetic studies on crop lodging have also been carried out using GWAS. On the contrary, although there are some GWAS reports on stalk lodging [13, 24], they are still relatively few, and the molecular mechanism of the variation of corn lodging-related traits is still poorly understood. High-throughput SNP markers have been widely used to identify genes controlling quantitative traits [25,26,27,28]. Genotyping by sequencing (GBS) is a relatively inexpensive method to obtain high-density markers for large populations taking the advantage of next-generation sequencing technologies [29,30,31,32].

In this study, an association mapping panel was genotyped by GBS. Based on this, association mapping was used to identify SNPs and excavate potential candidate genes on RPR, SBS, and morphological traits associated with stalk lodging resistance. The objectives of this study were to: (1) identify associated loci for RPR, SBS, and morphological traits of the stalk of maize; (2) ascertain stable SNPs and predict potential candidate genes in these regions; (3) dissect the genetic architecture of stalk lodging resistance-related traits.

Results

Phenotype analysis of the six lodging resistance-related traits

The phenotypes of all lodging resistance-related traits in the association mapping panel are shown in Table 1. The mean values of RPR, SBS, TID, and FID in the low plant density were higher than those in the high plant density. As for TIL and FIL, the mean values in the high plant density were higher than the mean values in the low plant density. For the six traits mentioned above, the skewness and kurtosis were less than 1, indicating that these traits followed a normal distribution. Furthermore, the coefficients of variation (CV) of these traits in the plant densities examined in this study ranged from 5.78–15.78% and 6.49–17.05%, respectively (Table 1).

Table 1 Phenotypic performance for related traits of stalk lodging resistance in the association mapping panel

ANOVA showed that the environment effects, density effects, genotype effects and interactive effects between the genotype and environment were both significant for six traits in the association mapping panel (Table 2). For the association mapping panel, the broad-sense heritability (h²_B) of all traits in low and high plant densities ranged from 0.59 to 0.72 and 0.61 to 0.71, respectively (Table 2), suggesting that variations of stalk strength traits were mainly controlled by genetic factors.

Table 2 Analysis of variance (ANOVA) for related traits of stalk lodging resistance under two plant densities in the association mapping panel

The results of the correlation analysis between the six traits of stalk strength at two densities for the maize inbred lines are shown in Fig. 1. In the correlation analysis, the consistency of all trait correlations between the two densities highly coincided with the results of GWAS. In addition, there was a strongly significant positive correlation between traits between SBS and RPR, SBS and TID as well as SBS and FID.

Fig. 1

Correlation analysis of lodging resistance-related traits under two plant densities in the association mapping panel. A and B stand for low plant density and high plant density, respectively. * Significant at P < 0.05. ** Significant at P < 0.01

GWAS for stalk lodging resistance related-traits

For RPR, a total of 29 significant SNPs were detected and located on chromosomes 1, 2, 3, 4, 5, 6, 7, 8, 9, 10 at all environments, which explained 11.10-16.07% of the phenotypic variation. For SBS, a total of 32 SNPs were detected across all environments, which explained phenotypic variation ranging from 9.29-17.69%. For other lodging resistance traits, the number of SNPs detected for TIL, TID, FIL and FID was 36, 53, 31 and 47, respectively, and accounted for phenotypic variation ranging from 12.31-20.72%, 11.23-18.50%, 13.96-23.59%, and 10.92%-17.44%, respectively (Table S1).

In total, 33 SNPs detected of different traits under same environment and density and explained phenotypic variation ranging from 11.23% to 20.70% (Table 3). Moreover, 2 significant SNPs for TIL were commonly detected across different environments, among which, Chr1_289271328 were identified in 2015BD, 2016BD and 2016SJZ at under high density and Chr2_54407952 were identified in 2016SJZ under low density and high density, with explanation of phenotypic variation range from is 14.97% to 18.14%. Moreover, one SNP, Chr2_233691764, was collocated for SBS, TID and FID on chromosomes 2 (Table 3).

Table 3 Important SNPs detected of different traits under same environment and density

To minimize the effect of environmental variation, the BLUP values were used to examine associations. In total, we identified the number of SNP for each trait by BLUP data, 6 for RPR, 3 for SBS, 10 for TIL, 8 for TID, 8 for FIL, 7 for FID at low plant density and 5 for RPR, 9 for SBS, 7 for TIL, 5 for TID, 7 for FIL, 6 for FID at high plant density (Fig. 2 and Table S2). The percentage of phenotypic variation explained by the identified SNPs (R²) for six traits ranged from 13.30 to 21.13% and from 10.10 to 21.01% at low and high plant densities, respectively (Table S2). The Manhattan plots and Quantile–quantile (Q-Q) plots between the six related traits of stalk strength at two densities are shown in Figs. 3 and 4. In addition, 14 important SNPs was detected of different traits at same density by BLUP value, which were located on chromosomes 2, 3, 4, 5, 8, 9 and 10 (Table 4).

Fig. 2

Stable SNPs were repeatedly detected in the two planting densities and the BLUP model, which were associated with six stalk lodging resistance-related traits. The significance threshold is –log10 (P-value) = 4.0. LD represent low plant density, HD represent high plant density, respectively. Purple represents third internodes length, Red represents fourth internode length, Blue represents third internode diameter, Orange represents fourth internode diameter, Yellow represents rind penetrometer resistance and Green represents stalk buckling strength, respectively

Fig. 3

Manhattan plots and QQ plots for the six traits at the low plant density. A Rind penetrometer strength. B Stalk bending strength. C Third internode length. D Third internode diameter. E Fourth internode length. F Fourth internode diameter

Fig. 4

Manhattan plots and QQ plots for the six traits at the high plant density. A Rind penetrometer strength. B Stalk bending strength. C Third internode length. D Third internode diameter. E Fourth internode length. F Fourth internode diameter

Table 4 Important SNPs detected of different traits by BLUP value

Candidate genes associated with significant SNPs

The physical locations of the SNPs were recorded using the B73 RefGen_v2 (www.maizesequence.org) based on the LD decay distance. A total of 346 candidate genes with gene descriptions were found (Table S3). The number of candidate genes involved in the six stalk lodging resistance related-traits of RPR, SBS, TIL, TID, FIL, and FID were 55, 78, 117, 37, 51, and eight, respectively. From the GO analysis results of the candidate genes in biological processes are mainly concentrated in the metabolic and cellular process, those influencing cellular component are mainly found in the intracellular and cellular anatomical entity, and those influencing molecular functions are mainly found in catalytic activity and binding (Fig. 5). As for the KEGG analysis of the candidate genes, a total of 13 pathways were identified (Fig. 6). These pathways included the carbon metabolism, ubiquitin mediated proteolysis, starch and sucrose metabolism, beta-alanine metabolism, pyrimidine metabolism, etc., which could be related to the stalk lodging. Among them, the pathway with the largest number of genes is the metabolic pathways, which have 36 candidate genes. Furthermore, we identified seven candidate genes to be associated with stalk lodging resistance (Table 5). Annotation information suggested that these candidate genes may control multiple traits during maize growth and development.

Fig. 5

GO-second class of candidate gene

Fig. 6

Analysis of KEGG pathway based on candidate genes (The figure was created by R version 3.6.1 based on KEGG pathway database www. kegg. jp/ kegg/ kegg1. html)

Table 5 Putative candidate gene of stalk lodging resistance-related traits

Discussion

Phenotypic variation, heritability, and correlations of traits

In general, obtaining an accurate measurement of phenotypic traits is essential to obtain reliable association results. The six traits investigated in this study exhibited large phenotypic variations with a normal distribution. A previous study showed that relatively high heritability will determine the power of QTL detection [33]. Our genetic analysis shows that the heritability of RPR and SBS ranged from 0.61 to 0.80. It was much higher than the range of 0.08–0.34 in a nested association population of maize [1]. The relatively high heritability in this study shows the predominant role of genetic factors for these traits.

There were significant correlations between each pair of stalk lodging resistance-related traits in this study, for instance: between RPR and SBS, which is consistent with previous results [13, 17]. Our study showed that the stalk strength traits decreased gradually with increasing density, which was consistent with previous findings [11, 34]. In the association mapping panel, a significant correlation was detected between SBS, TID, and FID. By contrast, the correlation between SBS, TIL and FIL was significantly negative, indicating that stalk strength traits are negatively associated with internode length and width at the population level. The above results suggest that some genetic factors were shared among these stalk lodging resistance-related traits.

Mapping analysis

Compared with traditional QTL mapping, GWAS covers a wide range of genetic diversity and more allelic polymorphisms, which could exploit the short linkage disequilibrium distance and help to pinpoint the functional genes of target traits using high-density molecular markers.

Hu et al. [8] detected ten QTLs for RPR and three QTLs for Internode diameter (InD) by applying the RIL population. In this study, we used GWAS to identify some RPR-related SNPs, among which Chr7_163048364 (bin7.04) and Chr8_88680106 (bin8.03) were located in the chromosomal region with Hu et al. [8]. In addition, Chr4_203233149 (bin4.08) and Chr8_67356036 (bin8.03) for TID and FID identified by the GWAS analysis locates exactly in the interval of the InD QTLs detected by Hu et al. [8]. Liu et al. [11] identified pleiotropic QTL, pQTL6-2, was association with RPR, whose confidence interval encompassed 16 QTLs, its genomic region is coincided with the physical position Chr6_158343036 (158 Mb) in this study. In addition, the SNP Chr1_272576164 (272 Mb) was detected association with SBS in this study also have same physical position with Liu et al. study. The remaining SNPs in this study were first reported to be associated with lodging resistance-related traits in maize.

Co-localization of SNPs for stalk lodging resistance traits

The SNP repeatedly detected in multiple environments is generally considered a stable SNP. Stably expressed SNPs detected in this study, five co-localized SNPs (Chr4_66017316, Chr4_16211307, Chr4_203233149, Chr4_236385528 and Chr8_130686461) were simultaneously identified under two plant densities. These stable SNPs were insensitive to the external environment and were hence considered to be important loci for the improvement of stalk lodging traits, as such, they can provide references for further gene cloning. Meanwhile, some specific SNPs were detected at high or low plant densities, respectively, which may be environmentally-specific loci requiring further genetic mapping.

From the comparison, we found some co-located locus in different densities in the same environment, but extremely few stable sites in different environments. The reason we detected less consistent loci in different environments may be because stalk strength trait itself is a relatively complex quantitative trait and is greatly affected by the environment. In addition, we found that the heritability of these traits is relatively low. This reason was further confirmed. From the results of the phenotypic correlation analysis, the correlation coefficient of both TID and FID was as high as 0.97 at both densities. Similarly, we located three SNPs (Chr4_16211307, Chr4_203233149, Chr8_130686461) associated with both TID and FID at both densities, this confirms the views of previous, phenotypic correlations between quantitative traits may derive from the correlation between QTL controlling them [35]. However, there were a large number of SNPs that did not co-located, indicating that lodging-related traits in maize seem to be controlled not only by several major QTLs but also by multiple micro-effect QTLs in specific locations or environments [36].

Candidate genes analysis

We identified 346 candidate genes in total located around common loci for stalk lodging resistance-related traits, which are involved in a variety of biochemical metabolic pathways. Based on the information of the gene model on MaizeGDB (Table S3), seven potential candidate genes related to RPR, SBS, TIL, FIL and FID were obtained (Table 5). Notably, some candidate genes correlated to stalk lodging-related traits were related to cellulose and lignin biosynthesis, essential for the cell wall development in the plant stem. For instance, beta-amylase (AMY), beta-glucosidase (GLU), UDP-glycosyltransferase (UGT), and protein kinase played an essential role in the synthesis of cell wall components [37]. Indeed, modify the expression of a transcription factors by changing the mRNA abundance of downstream target genes to change the biosynthesis of lignin and he lodging resistance of stalk can be increased [38]. Interestingly, seven candidate genes were found to be related to cell wall components in this study (Table 5). GRMZM2G074792, which is located in Chr6_158343036 of RPR, encodes xyloglucan glycosyltransferase and related to plant cell wall cellulose synthesis, which is the major source of cellulose-harbours enzyme [39]. GRMZM2G300412, encoded for UDP-glucuronic acid decarboxylase, was located in Chr1_272576164 of SBS, involving in metabolic pathways and amino sugar and nucleotide sugar metabolism. GRMZM2G072526 was located in Chr7_160255239 and Chr7_160255241, controlling SBS, whose encoded glucan endo-1,3-beta-glucosidase is mainly involved in carbohydrate metabolism, it is associated with cell wall synthesis, which may be related to maize lodging. Previous studies demonstrated that UDP-glucuronic acid decarboxylase was a key enzyme in the synthesis of UDP-xylose for the formation of xylans during cell wall biosynthesis [40]. GRMZM2G111344, was located in Chr5_15958677 of TIL, encoding for UDP-glycosyltransferase (UGT), involved in flavonoid biosynthesis and biosynthesis of secondary metabolites. According to previous studies, UGT was the key precursors of cell wall carbohydrates [37]. These descriptions indicate that regulation of the expression of these genes may affect cell wall formation. The candidate genes GRMZM2G007899 and GRMZM2G311059, were located in Chr10_139852648 of TIL, showed high expression of MYB transcription factor had increased ectopic lignin and the xylem vessels were regular and open, are related transcriptional activators of the lignin biosynthetic pathway during secondary cell wall formation in Arabidopsis [41, 42]. In rice, GRMZM2G021051 was located in Chr2_233691559 of FIL, whose the homologous with shortened basal internodes, is a new rice lodging-resistance gene and encodes a gibberellin (GA) 2-oxidase and can control the elongation of internodes at the base of the stem by regulating the activity of the GA [43]. GRMZM2G408462, which is located in Chr3_212705423 of FID, encoded for WRKY transcription factor, whose directly regulate expression of the major monolignol biosynthetic genes and genetic modification of genes involved in lignin biosynthesis [44, 45]. Although the role of these genes in maize requires further investigation, they should be used as target sites for the development of maize lines resistant to lodging.

Conclusion

In this study, we identified 6, 3, 10, 8, 8, 7 SNPs associated with RPR, SBS, TIL, TID, FIL, FID at low plant density and 5, 9, 7, 5, 7, 6 SNPs associated with RPR, SBS, TIL, TID, FIL, FID at high plant density, respectively, via GWAS. Most markers were located within or close to QTLs identified in previous studies. We were particularly interested in the seven potential candidate genes that were predicted based on functional annotations, but further investigation is needed for verification of this hypothesis. These findings shed light on the genetic basis of six stalk lodging resistance related-traits, and candidate genes could be used for further positional cloning.

Materials and methods

Plants materials and field experiments

A total of 248 diverse maize inbred lines were used to form an association mapping panel. All lines were grown according to the split-plot set two densities, two replicates for each density, and a low density of 75,000 plants/ha and a high density of 105,000 plants/ha. The work was performed at the Experimental Station of Hebei Agricultural University in Baoding and Shijiazhuang in 2015 and 2016. For each replicate, each line was grown in a 3-m long single-row plot, with a 0.6-m inter-row spacing. All of the plant materials used in our study were derived from the China Agricultural University and National Maize Improvement Center of China.

Phenotype evaluation

Based on previous studies on stalk lodging resistance in maize, we decided to measure morphological traits and stalk strength during one week after grain filling [11]. Five representative plants of each line from each replicate were selected for evaluation and the mean values for each line were computed for each trait. The TIL, FIL, TID, and FID were measured using electronic micrometers. At the same time, morphological characters were measured with the same material, RPR and SBS were measured in the middle of the flat side of the third and fourth internode of the stalk using a stalk strength appliance YYD-1 (Zhejiang TopuYunnong Science and Technology Co., Ltd, Zhejiang, China). At the base of the stem, the middle part of the third and fourth internodes is inserted at a constant speed and perpendicular to the direction of the stem, and the maximum penetration of the stem epidermis is read. Similarly, the bending strength of the stalk is also pressed at the center of the stalk at a uniform speed, and the force should not be too strong and record the value. The range of measurement was between 5 and 500 N, with a resolution of 0.1 N; reported units of RPR and SBS are in N/mm² and N, respectively.

Statistical analysis of phenotypic data

The mean value of each inbred line for each trait was used for descriptive statistical analysis. Analysis of variance (ANOVA) was carried out with SPSS19.0 for related traits of stalk lodging resistance under two plant densities in the association mapping panel. Broad-sense heritability (h²_B) was calculated according to Knapp et al. [46].

$${h}_{B}^{2}={\sigma }_{g}^{2}/({\sigma }_{g}^{2}+{\sigma }_{ge}^{2}/e+{\sigma }_{\varepsilon }^{2}/re)$$

where ${\sigma }_{g}^{2}$ is the genetic variance, ${\sigma }_{ge}^{2}$ is the interactive effect of genotype × environment, ${\sigma }_{\varepsilon }^{2}$ is the error variance, e is the number of environments, and r is the number of replications in a given environment.

The best linear unbiased prediction (BLUP) of the phenotypic values of each line was calculated across all environments using the R package “lme4” [47]. The BLUP value of each line was used for the GWAS analysis. The correlation analysis was performed using the “Performance Analytics” package in R.

Genotyping

The GBS method was used to genotype the 248 inbred lines of the association panel [29]. First, the genomic DNA was extracted from leaves of maize under normal growth conditions using the cetyltrimethylammonium bromide (CTAB) method [48]. The DNA concentration and integrity were measured with NanoDrop 2000 instrument (Thermo Fisher Scientific, Waltham, MA, USA) and agarose gel electrophoresis, respectively. The extracted DNA of each line was digested using the restriction enzyme ApeKI and ligated with bar code. The DNA samples of certain numbers were mixed, purified, amplified, purified again, and chosen according to fragment length. Those fragments were evaluated using the length test, Paired-End-Tag by Illumina Hiseq2000. Then selected sequences were aligned to the B73 reference genome (the second version) using the BWA software, followed by SNP calling using Samtools [49]. SNPs with a missing rate < 0.2 and minor allele frequency (MAF) > 0.05 were selected. Finally, a total of 83,057 SNPs were used for the GWAS analysis. The PLINK 1.90 beta software was used to estimate LD between pairs of SNPs within 200 kb in the genomic region based on the Hill and Weir method [50, 51]. The LD decay distance for this association mapping panel was 120 kb (r² = 0.1) based on previous study [52, 53]. The population structure (Q) was estimated using the software Admixture 1.3, while kinship (K) was estimated using Analysis-Kinship in Tassel 5.0.

Genome-wide association studies

GWAS data was analyzed with the mixed linear model (MLM) using the “GAPIT” package in R. The SNP markers of six stalk lodging resistance related-traits in the association mapping panel together with the Q and K matrix were used as covariates to decrease spurious association and detect marker loci combining with target traits. The GWAS analysis is performed with a Bonferroni correction, however this was found to be too strict for less significant trait associations. Therefore, we reduced the significance threshold to–log10 (P) ≥ 4 for all traits [28].

Prediction of candidate genes

The candidate gene analysis was based on the maize inbred line B73 reference genome version v2 (centering on the marker site and extending 120 kb upstream and downstream) and searching for the information and functions of the candidate genes on the MaizeGDB genome browser (http://www.maizegdb.org/) and NCBI website (https://www.ncbi.nlm.nih.gov/). Gene ontology (GO) enrichment analysis was performed using the Gene ontology website (http://www.geneontology.org/). Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment analysis was performed using the KOBAS version 3.0 (http://kobas.cbi.pku.edu.cn/kobas3/?t=1) [54].