Identification of Candidate Genes for Drought Tolerance at Maize Seedlings Using Genome-Wide Association.

BACKGROUND: Drought stress is a serious threat that limit maize growth and production.
OBJECTIVES: The assessment tolerance level of maize by measuring changes in the main biochemical and physiological indicators under drought stress.
MATERIAL AND METHODS: We performed a genome-wide association analysis of biochemical and physiological indicators using an elite association panel.
RESULTS: The results revealed that eight significant SNPs (p<0.05/N) located in eight genes that are distributed on different chromosomes were associated with drought resistance indices under drought stress. Among these genes, four genes were linked via the associated SNPs with drought-resistance indices of the malondialdehyde activity (MDA), three genes were linked with drought resistance indexes of the superoxide dismutase activity (SOD), and one gene was linked with drought resistance indexes of relative conductivity (REC). The candidate genes functioned as transcription factors, enzymes, and transporters, which included trehalase, the AP2/EREB160 transcription factor, and glutathione S-transferase and also encoded a gene of unknown function. These genes may be directly or indirectly involved in drought resistance. The expression levels of ZmEREB160 responded to ABA and drought stress.
CONCLUSIONS: These results provided good information to understand the genetic basis of variation in drought resistance indices of biochemical and physiological indicators during drought stress. Copyright:

1. Background

Maize irrigation requires a large amount of water at all stages, and water is one of the most important factors for limiting maize growth during all developmental stages from seed germination to reproduction. However, water deficit is a serious problem worldwide. It is necessary to understand the genetic mechanisms of drought tolerance at the seed germination and seedling stage in maize ( 1 - 2 ). Ongoing breeding efforts have focused on developing more drought tolerant crops. The early seedling stage is the most crucial for plant establishment ( 3 ).

When plants experience water deficit, they undergo a variety of molecular, biochemical and physiological changes to cope with drought stress through the regulation of gene expression. These mechanisms involve complicated biological processes that finely coordinate molecular signaling pathways to enhance survival and reproduction under drought stress by encouraging the development of more compact leaves and harnessing metabolism to produce protective compounds ( 4 - 6 ). The changes of the metabolism of plants in response to drought are complex and involve metabolic pathways and physiological and biochemical process. Because drought stress is the accumulation of a large number of reactive oxygen species (ROS) in plants. Hyperosmotic stress, in turn, leads to secondary oxidative stress. This is believed to disturb the balance between the formation and removal of ROS. To reduce high levels of ROS and accumulate osmolytes to regulate cellular osmolality, many plants have developed sophisticated enzymatic antioxidant defense systems that cause the accumulation of free amino acids such as proline (Pro) and malondialdehyde (MDA) and an increase in superoxide dismutase activity (SOD) ( 3 , 7 , 8 ). Proline, superoxide dismutase activity and malondialdehyde content affect drought tolerance levels, and the accumulation of osmolytes is a protective strategy against abiotic stress in plants. They are often measured as an indicator for drought tolerance ( 9 ).

With the development of biotechnology, the cost of sequencing the maize genome is getting lower and lower. The genome-wide association study (GWAS) is a powerful tool to dissect the genetic basis of the association of primary metabolites with trait-associated loci in maize ( 10 , 11 ). Compared by the traditional screening methods, genotyping, and sequencing technologies are more efficient; therefore, measuring the effects of genes on drought tolerance in maize is a convenient method for studying metabolic traits that vary during drought to unravel the genetic basis of drought tolerance. In maize, GWAS has been successfully utilized to identify numerous candidate genes that control a number of metabolic or morphological traits, such as husk traits ( 12 ), stalk cell walls components ( 13 ), oil biosynthesis ( 14 ), and agronomic traits related to drought stress ( 15 ).

Single nucleotide polymorphisms (SNPs) in inbred maize lines were mapped to the same reference sequence, and comparison with the reference genome B73 enabled the identification of SNPs. Associated genetic markers may be causal for the trait of interest or in linkage disequilibrium with a causal locus ( 16 ). In quantitative trait locus (QTL) studies, QTLs were typically localized to 10-20 cM intervals, and populations were constructed according to recombination mapping ( 17 ). Complex trait analysis can benefit from a number of statistical adjustments, which can decrease the noise from microenvironmental nongenetic sources in an experiment, making the results more robust and powerful ( 15 ).

2. Objectives

In this study, an association mapping panel with 68 diverse inbred maize lines was used to conduct an association analysis of drought tolerance at the seedling stages. The main goal was to identify drought-resistant genes involved in drought stress tolerance at the seedling stage. Drought-resistant genes can be developed into functional markers. This information helps for the marker-assisted selection of maize varieties for further use in maize breeding programme.

3. Materials and Methods

3.1. Plant Materials

To identify genes of response drought stress at the seeding stage in maize, a panel of 68 inbred maize lines was used for the whole-genome association analysis (Supplementary Table S1). The association panel included the parental lines of widely used commercial hybrids in China. The germplasm was used in this experiment to evaluate physiological and biochemical indicators in maize under drought stress in greenhouse in the summer. Young leaves from the three-leaf stage were collected in bulk and ground in liquid nitrogen. DNA was manually isolated from the leaves using a modified CTAB method ( 18 ). The DNA concentration and quality were evaluated by a Nano Drop 2000 (Thermo Scientific).

3.2. Trait Measurement and Phenotypic Data Analysis

For drought treatment, the seeds were sown in pots containing vermiculite and nutritional soil, and seedlings were watered with tap water until they reached the three-leaf stage. The seedlings were used as experimental materials. The control seedlings continued to be watered, while water was withdrawn from the drought treatment seedlings for 10 d. Three replicates were performed for each treatment. All experiments were carried out in a greenhouse; the seedlings were harvested prior to measuring the physiological and biochemical indicators. The measured physiological and biochemical traits included Pro, MDA, SOD, relative conductivity (REC) and leaf relative water content (LRWC). The drought resistance indices of the physiological and biochemical indicators were calculated as Bellows:

1) LRWC: drought resistance index = trait mean in drought condition/trait mean in control condition; 2) Pro and SOD: drought-resistance index = trait mean in drought condition - trait mean in control condition)/ trait mean in control condition; 3) REC and MDA: drought-resistance index = 1-(trait mean in drought condition - trait mean in control condition)/ trait mean in control condition. All data was conducted by analysis of variations.

3.3. Genotype, Quality Control, Population Structure and Estimation of the Kinship Matrix

The genotyping was performed for 68 inbred maize lines using the Illumina MaizeSNP55 K Bead Chip (Illumina). The maize 55 K SNP Affymetrix Axiom Genotyping Array contains 55,229 SNPs, and these SNPs are distributed throughout the whole maize genome. The SNP density ranged from 22 to 27 SNPs per Mb p. PLINK 1.07 software ( 19 ) was used to calculate the frequency. The SNP selection criteria were previously described ( 4 ). SNPs that were missing in more than 20% of the association-mapping panel and those with a minor allele frequency (MAF) <0.05 were excluded from association-mapping panel, and the reminder of the high-quality SNPs were used for the association analysis. A total of 40,580 SNPs with a MAF >0.05 and a missing rate <20% were used for GWAS analysis. A total of 2,642 SNPs was used for the population structure assessment, which was performed with STRUCTURE v2.3 software ( 20 ). Five runs were performed in STRUCTURE for each population (k) (1 to 10) ( 21 ). Using SPA Gedi software ( 22 ), the 2,642 SNPs were used to generate the relative kinship matrix (K). The kinship matrix measured the genetic similarity between individuals.

3.4. Association Mapping

A total of 40,580 SNPs was selected for the GWAS by combining the traits affected by drought stress. Linkage disequilibrium measurement was used to estimate the LD between each chromosome using TASSEL V3.0 software ( 23 ). A mixed linear model (MLM) was evaluated in TASSEL, and both the K and Q matrices were taken into account to avoid spurious associations. Quantile-quantile plots (QQ plots) show GWAS results of effect measures. The QQ plots and Manhattan plots were displayed using R software. In this study, the Bonferroni test (0.05/N, N=total SNPs used) criterion is typically used as a very strict threshold. The SNPs were considered according to the significance of the SNP-trait associations with the identified candidate genes. The extended region where the LD was decayed was searched for candidate genes.

3.5. Annotation of Candidate Genes

The SNPs with the most significant associated traits were selected. The physical positions of the SNPs were based on Maize B73 RefGen_V3. The B73 genome was used to identify the locations of the SNPs relative to those of the candidate genes. The corresponding genes were annotated by using the MaizeGDB database. To evaluate the candidate genes associated with various physiological and biochemical traits, we chose ZmEREB160 as a candidate gene associated with drought resistance.

3.6. ZmEREB160 Expression Under Drought Stress

To examine the expression patterns resulting from various drought treatments, seeds of the inbred maize line Ji853 were sown in vermiculite. Three kinds of treatments, dehydration, high salinity and abscisic acid (ABA), were separately applied to plants in the three-leaf stage. The roots of the seedlings were immersed in solutions containing 20% (w/v) PEG6000 and 100 mM ABA. Samples from a minimum of five seedlings were collected immediately and 3 h after ABA and PEG treatment. All collected samples were immediately frozen in liquid nitrogen and stored at -80 ℃ for prior to RNA extraction. Total RNA was extracted from the above described treatment samples by using a HiPure Plant RNA Mini kit (Magen, China), and first-strand cDNA synthesis was performed with M-MLV reverse transcriptase (Invitrogen) according to the kit protocol. Gene specific primer pairs were designed for the sequence of ZmEREB160 (Supplementary Table S2), and actin was used as the internal reference gene. qPCR was performed using a real-time PCR system (Light Cycler 480, Roche) with SYBR Green I mix (TaKaRa). The 2-ΔΔCt quantification method ( 24 ) was used and the variation in expression was estimated from three technical replicates.

4. Results

4.1. Analysis of the Phenotypes of the Inbred Lines Panel

To discover drought-resistant genes by association analysis, physiological and biochemical traits were tested in a subset of 68 inbred maize lines. Five traits related to drought tolerance were measured, including MDA, SOD, Pro, REC and LRWC at the seedling stage, the mean values of drought resistance indexes of MDA, SOD, Pro, REC and LRWC were calculated in elite inbred maize lines subjected to drought treatment. Extensive phenotypic variation was observed for MDA, SOD, Pro, REC and LRWC during drought stress. The descriptive statistics for the phenotypes related to drought stress are presented in Table 1. A wide range of variation was observed among the accessions in terms of the drought resistance indexes of physiological and biochemical traits.

Table 1

drought-resistance indexes in seedlings stage

Trait	Mean	SD	Minimum	Maximum	CV (%)	Skewness	Kurtosis
SOD	0.43	0.29	0.01	1.04	67.44	0.641	-0.669
MDA	0.67	0.19	0.31	0.99	28.36	-0.244	-0.846
Pro	0.62	0.49	0.03	1.92	79.03	0.973	0.190
REC	0.62	0.30	-0.06	1.14	48.39	-0.548	-0.499
LRWC	0.83	0.16	0.45	0.99	19.28	-0.818	-0.309

Trait: SOD, superoxide dismutase activity; MDA, malondialdehyde content; Pro, Proline content; REC,relative conductivity; LRWC, leaf relative water content.

SD, standard deviation; CV(%), coefficient of variation.

For traits under drought stress, the drought resistance index of MDA was 28.36. The drought resistance index of REC was 48.39. The drought resistance index of SOD was 67.44. The coefficients of variation (CV %) of the five drought resistance indexes of the traits ranged from 19.28 (LRWC) to 79.03 (Pro) under drought stress (Table 1). The drought resistance index of the MDA trait was evaluated and valued ranging from 0.31 to 0.99. The drought resistance index of the SOD trait ranged from 0.01 to 1.04. The results of the correlation analysis showed that there were significant positive correlations between SOD and MDA (p< 0.05) (Table 2). The drought resistance index of the Pro trait varied from 0.03 to 1.92. The drought resistance index of the REC trait ranged from -0.06 to 1.14. The drought resistance index of the LRWC trait varied between 0.45 and 0.99. Pro was highly negatively correlated with LRWC (p< 0.01) (Table 2). The results showed that highly significant genetic variation in the seeding stage was observed among the association panels under drought stress.

Table 2

drought-resistance indexes in seedlings stage

	MDA	SOD	REC	Pro	LRWC
MDA	-
SOD	0.242*	-
REC	0.127	0.191	-
Pro	0.082	0.120	-0.126	-
LRWC	0.033	0.131	0.075	-0.315**	-

,

: significant at the 0.05 and 0.01 probability level, respectively.

,

4.2. SNP Genotyping, Population Structure and Kinship

Using the maize 55 K SNP Affymetrix Axiom Genotyping Array, the genotypes were determined for the panel of 68 inbred lines. After excluding data with a MAF <0.05 and a missing data value >20%, data for 40,580 SNPs was used for further analysis. Based on the nucleotide polymorphisms, the population structure was assessed with STRUCTURE software for K values ranging from 1 to 10 for the entire panel using 2,642 SNPs among the 40,580 SNPs that were distributed throughout the whole maize genome. A significant increase was observed when the delta k value was changed. A sharp peak in the ΔK value was identified. As a result, the delta k reached its maximum value at k=5 (Fig. 1A). The largest value of the ΔK statistic was used as an indicator for evaluating the most likely number of subpopulations. Thus, the maize panel was divided into five subgroups based on the structure analysis (Fig. 1B).

Figure 1

Analysis of the population structure and relative kinships. A) Relative kinships: Δk value was based on 2642 SNPs; B) Population structure estimated by STRUCTURE (k =5).

4.3. Genome-Wide Association Studies and Candidate Gene Analysis

To dissect the genetic basis of the physiological and biochemical traits, we selected 40,580 SNPs across the entire maize genome for the GWAS of drought resistance indexes during drought stress. As shown in the Manhattan plots (Fig. 2), the criteria used for the Bonferroni multiple test correction for detecting associations related to the traits were very strict. The results of the GWAS revealed that 8 SNPs were significantly associated with the different traits. Those SNPs satisfied the threshold for significance applied in this study (-log10P value > [-log10 (0.05/40580)]). For the traits, four SNPs were significantly associated with MDA, three SNPs were associated with the SOD trait, and one SNP was associated with the REL trait (Table 3). Four SNPs were significantly associated with the MDA trait, and the 4 associated SNPs had –log10 (P) values that ranged from 14.9 to 31.6, and could explain 14.6-68.4% of the phenotype variation. For the SOD trait, the 3 associated SNPs had –log10 (P) values that ranged from 5.4 to 7.1, and could explain 83.1-85.2% of the phenotypic variation. For the REC trait, we detected one GWAS peak on chromosome 9. The SNP explained 85.2% of the phenotypic variation in the association panel. The MLM was then generated, and both population structure and kinship relationships were taken into account to avoid spurious associations. Compared by the MLM model, the GLM model does not consider the effect of kinship relatedness among individuals, while kinship relatedness is an indispensable effect included in the calculations for the MLM, the quantile-quantile plots of the GWAS results indicated that the model was well fitted to the data (Fig. 2). When we used the MLM model, Pro and LRWC were not significantly associated with any SNP as indicated by Bonferroni multiple test correction. If we included all a subpopulation, the Pro trait was associated with 12 very significant SNPs as indicated by Bonferroni multiple test correction (Supplementary Fig. 1A).

Figure 2

Manhattan plots (left) and quantile-quantile (QQ) plot (right) of GWAS for MDA, REC and SOD. A) Manhattan plot and QQ plot for MDA. B) Manhattan plot and QQ plot for SOD. C) Manhattan plot and QQ plot for REC.

Table 3

SNPs, chromosomal position and candidate genes significantly associated with drought-resistance indexes under drought stress

Trait	SNP	Chr.	Position (bp)	R2 (%)	P-valuea	Gene	Annotation
MDA	AX-116874493	3	82727904	68.4	2.67E-32	GRMZM2G081380	Ca2+-binding actin-bundling protein
	AX-90823123	3	113180327	61.7	6.20E-26	AC185252.4_FG002	Unknow
	AX-86278443	1	80073680	22.5	8.19E-21	GRMZM2G162690	trehalase
	AX-86307466	1	186655481	14.6	1.06E-15	AC208564.3_FG004	MADS transcription factor
REC	AX-123949274	9	11501257	85.2	6.32E-8	GRMZM2G171179	AP2/EREB160 transcription factor
SOD	AX-86272995	6	30931897	85.2	7.20E-8	GRMZM2G072606	Unknow
	AX-91220227	4	169065243	83.1	1.96E-7	GRMZM2G152278	Calcium-dependent lipid-binding protein
	AX-90538327	3	91677989	84.4	3.79E-6	GRMZM2G328374	glutathione S-transferase

R2,the percent of of phenotypic variance explained by the SNP.

Significant SNP-trait associations with a stringent threshold of –log (0.05/40,580).

Furthermore, we analyzed the drought resistance indexes of LRWC using a custom level, and 1 significantly associated SNP was identified (Supplementary Fig. 1B). These results indicate that the association panel provided sufficient statistical power to detect alleles associated with physiological and biochemical traits.

To identify candidate genes related to the drought resistance indexes of the physiological and biochemical traits, we searched the maize B73 reference genome based on the significant trait-SNP associations. We identified 8 candidate genes associated with stress-related traits (Table 3). The GWAS analysis revealed a site that was significantly associated with drought tolerance on chromosome 3 that contained three trait-associated-SNPs. The physical distance between these candidate genes and the significant SNPs varied from 8 to 41 kb. Among these genes, trehalase and the MADS transcription factors were associated with MDA under drought conditions. OsMADS18 responds to ABA signaling and the MADS transcription factor enables the recognition of similar DNA sequences as common DNA-binding domains to execute their regulatory roles in the transcription of genes in response to stress ( 25 ). Transgenic Arabidopsis plants overexpressing CaMADS showed enhanced tolerance to cold, salt and osmotic stress ( 26 ). OsTRE1 is a trehalase, the overexpression of which enhanced salt tolerance in rice (Islam, Kato et al. 2019). The SNPs explained approximately 14.6 to 68.4% of the phenotypic variation (Table 3). For SOD, three genes (including a calcium-dependent lipid-binding protein, glutathione S-transferase, and one gene with an unknown function) were identified. OsGSTU3, a known and thoroughly investigated gene, modulates a regulatory network that contributes to drought tolerance ( 27 ). These SNPs were found on chromosomes 3, 4, and 6, and explained approximately 83.1 to 85.2% of phenotypic variation. For REC under drought stress, one SNP on chromosome 9, AX-123949274, was found to be significantly associated with the REC trait. This SNP was near the gene AP2/EREB160 (ethylene responsive element binding protein), which was identified as encoding a putative stress-responsive NAC transcription factor. In the ABA-independent pathway, AP2/EREB transcription factors activate cis-elements that are present in the promoters of stress-induced genes. The expression of GmAP2/EREB genes was upregulated in beans under water-deficient conditions ( 28 ). Liu et al. found a significant association between the ZmDREB2.7 gene and drought responses in the seedling stage using a diverse population of maize consisting of 368 varieties from tropical and temperate regions ( 29 ). Out of the 8 identified SNPs, we chose one candidate gene associated with the REC trait that could serve as a potential target for future genetic study. A number of studies have reported that the family to which this gene belongs has been proven to be related to drought resistance. The results showed that there is wide genetic variation in seeding stage traits in maize germplasm. Therefore, we focused on the gene significantly associated with the REL trait.

4.4. ZmEREB160 Expression in Response to Abiotic Stress

To identify the function of the ZmEREB160 gene, we performed expression level analysis in plants under various abiotic treatments. As expected, ZmEREB160 expression levels rapidly increased in roots and during ABA treatment (Fig. 3), ZmEREB160 expression levels were remarkably increased by PEG treatment in roots and leaves. The results showed that ZmEREB160 is induced by ABA and drought stress, and the induction of ZmEREB160 expression in response to drought may be mediated by the ABA-dependent pathway.

Figure 3

Expression patterns of ZmEREB160 under abiotic stress in maize. Expression of ZmEREB160 in leaves (white bars) and roots (grey bars), respectively, under stress treatments with PEG (A) and ABA (B). The samples were collected for each time point per treatment.

5. Discussion

Drought may threat the survival of plant and ultimately reduce crop yields, inbred maize lines showed different drought resistance under drought stress because of remarkable genetic differences. These potential genetic variations can be used to breed maize with enhanced salt or drought tolerance. The seedling stage is the most vital for plant establishment. The specific mechanisms underlying drought and salinity tolerance during the seedling stage have not been elucidated.

For the purpose of the current study, it would have been ideal to test whether the chosen traits were appropriate for GWAS of quantitative traits. Moreover, this may have helped to save a large amount of money and time. As an alternative, we selected adaptive traits that are known to be linked to drought tolerance in plants. To determine the genetic basis and the underlying molecular pathways, we analyzed the associated drought resistance indexes of physiological and biochemical traits in an association panel, and a set of putative candidate genes were identified by performing GWAS analysis. For the GWAS analysis, we used 40,580 SNPs identified by the Illumina MaizeSNP55 K Bead Chip to explore the genetic architecture that forms the basis drought tolerance in maize and to survey desirable alleles for maize improvement. Because the population was small, there was the risk of missing real associations. We chose the Bonferroni test criterion as a threshold to avoid false positives in the GWAS results. A total of 8 SNPs was found to display significant associations with three important agronomic traits. For some of the traits, because of the costs and/or labour required for data collection, it would be costly and require a large number of fields to evaluate the chosen traits, so GWAS was effective in determining whether we should continue to measure the chosen traits in crops. The GWAS results indicated that MADS transcription factors may be associated with stress tolerance. The expression levels of ZMM7-L, a MADS transcription factor, were upregulated under cold-, NaCl- and PEG-induced stress and downregulated under ABA-induced stress in maize ( 33 ). Glutathione S-transferases (GSTs) were found to be regulated by various abiotic stresses, including dehydration and cold and salt stress. OsGSTU4 overexpression in transgenic Arabidopsis improved tolerance to salinity and oxidative stress ( 34 ). OsGSTU30 overexpression in plants improved drought tolerance and OsGSTL2 overexpression plants enhanced cold, osmotic and salt tolerance in Arabidopsis ( 27 , 35 ). Trehalase biosynthetic enzymes were involved in enhancing osmotic stress in plants, and AtTRE1 overexpression in Arabidopsis improved drought stress tolerance ( 36 ). OsTRE1 transgenic rice showed enhanced salt tolerance compared with wild type rice ( 37 ).

Accurate phenotyping is the most vital factor for GWAS of maize, different genetic and physiological mechanisms of drought tolerance exist in the diverse germplasms of maize, and drought tolerance mechanisms can vary among varieties ( 2 ). The effects of soil conditions on the seedling stage are extremely complex, and involve various physical and biochemical cues ( 8 ). Proline serves as a molecular chaperone, an osmolyte, and an ROS scavenger, all of which helps to maintain cell homeostasis during drought stress. Drought stress accelerates the accumulation of reactive oxygen species (ROS) in plant cells, and the imbalance between ROS production and detoxification accounts for oxidative stress. In response to oxidative stress, an effective antioxidant system can act to scavenges excessive ROS ( 3 , 38 ). We used GWAS to identify SNPs associated with the drought resistance indexes of physiological and biochemical traits, and GWAS provided much more precise positional estimates than the QTL method when mapping populations. The appropriate model was selected to obtain a higher level of confidence in the association results ( 39 ).

In support of this approach, we selected candidate genes based on the SNPs that have been reported to be tightly associated with the corresponding traits in gene homologue families. Among 368 maize varieties, DREB genes were found to be associated with resistance to drought stress ( 40 ). In maize, a family of AP2 genes expressed throughout plant development and in various physiological processes were associated with abiotic stress responses ( 41 ). Therefore, we regarded the AP2 gene as a candidate gene associated with drought resistance. It would be interesting to test whether the expression of this gene leads to increased drought stress. The expression levels of ZmEREB160 were significantly increased by ABA, and PEG treatment. The results indicated that ZmEREB160 responded to various stresses. These findings have important implications for understanding the genes associated with drought tolerance in the seedling stage and exploiting them to clone genes conferring drought tolerance. Therefore, understanding the genetic control of maize traits and applying that knowledge in maize breeding programs might be instrumental in developing improved germplasm in maize ( 42 ).