Natural variations in the non-coding region of ZmNAC080308 contributes maintaining grain yield under drought stress in maize

Background Natural variations derived from both evolutionary selection and genetic recombination, presume to have important functions to respond to various abiotic stresses, which could be used to improve drought tolerance via genomic selection. Results In the present study, the NAC-encoding gene of ZmNAC080308 was cloned and sequenced in 199 inbred lines in maize. Phylogenetic analysis showed that ZmNAC080308 is closely clusteredinto the same group with other well-known NAC genes responding to improve drought tolerance. In total, 86 SNPs and 47 InDels were identified in the generic region of ZmNAC080308, 19 of these variations were associated with GY (grain yield) in different environments. Nine variations in the 5’-UTR region of ZmNAC080308 are closely linked, they might regulate the gene expression and respond to improve GY under drought condition via Sp1-mediated transactivation. Two haplotypes (Hap1 and Hap2) identified in the, 5’-UTR region using the nine variations, and Hap2 containing insertion variants, exhibited 15.47 % higher GY under drought stress condition. Further, a functional marker was developed to predict the drought stress tolerance in a US maize inbred line panel. Lines carrying Hap2 exhibited > 10 % higher GY than those carrying Hap1 under drought stress condition. In Arabidopsis, overexpression ZmNAC080308 enhanced drought tolerance. Conclusions ZmNAC080308 is an important gene responding to drought tolerance, a functional marker is developed for improving maize drought tolerance by selecting this gene. Supplementary Information The online version contains supplementary material available at 10.1186/s12870-021-03072-9.


Background
Food security in terms of grain yield will become more difficult to achieve with limited arable land as human populations grow. Crop yields are suffered from different abiotic stresses and it will continue to be challenged due to climate changes. The latest report from FAO indicates billions of US dollars were lost in agricultur from 2005 to 2015 due to environmental disasters. More than 30 % of these losses were caused by drought, which alone resulted in 29 billion US dollars in agricultural losses (http://www.fao.org/news/story/en/ item/1106977/icode/). Therefore, improveing the crop tolerance to abiotic stresses, especially drought, becomes very urgent.
Transcription factors (TFs) play important roles in molecular mechanisms of stress tolerance in plants [1]. In particular, NAC is a member of a large plant-specific family of TFs that includes NAM, ATAF, and CUC with important diverse functions. The first NAC gene, i.e. NAM (no apical meristem), was identified in petunia [2]. Subsequently, CUC1 and CUC2 were identified in Arabidopsis,because the N-terminus of CUC2 was highly conserved with those of NAM and ATAF, this Nterminal region was then named the NAC domain [3]. As shown in the TF database (http://planttfdb.cbi.pku. edu.cn/), more than 500 NAC genes have been identified in Arabidopsis, rice, and maize, and these genes have been categorized into different subgroups [4][5][6]. Therefore, NACs involved in the plant stress-responsive subgroups are known as SNACs (stress-responsive NACs) [7][8][9][10]. These SNACs are involved in many stressresponsive metabolic pathways in plants including ROS generation, hormone function, osmotic adjustment, Ca 2+ signaling, secondary metabolic processes, etc. [11], where activated NACs can regulate the expression of downstream elements.
Recently, researchers have focused more on the function of NACs in abiotic stress responses, especially in crops. So far, most of the drought-related NACs reported in major crops are positive regulators. Overexpression of ZmNAC55 in Arabidopsis can enhance drought tolerance compared to wild-type lines [12]. TaNAC2 is a NAC TF in wheat that could enhance resistance to drought, salt, and freezing stresses in TaNAC2-transgenic Arabidopsis [13]. The root-specific expressed rice gene OsNAC10 can improve resistance to abiotic stress when it is overexpressed. Moreover, overexpression OsNAC10 in rice increase 25-42 % GY under drought stress conditions in the field [14]. Several other NAC TFs, such as SNAC1 [15], ONAC045 [16], SNAC2 [17], OsNAC5 [18], OsNAC6 [19], OsNAP [20], ZmNAC111 [21], and TaNAC29 [22], have also been identified, these genes have important functions responding to abiotic stresses tolerance, which are useful for developing stress-tolerance crops.
Identification natural variations in genes associated with target traits becomes more important, especially for quantitative traits including drought tolerance, as some cis-or trans-regulatory functions involved in responses to abiotic stress tolerance can be lost or gained, due to natural variations in these genes [23]. In maize, many natural variants associated with drought tolerance have been identified. Two InDel variations in the upstream region of ZmVPP1 were identified, which were significantly associated with drought resistance in maize. The haplotype of the insertion variants containing three MYB-binding domains, has greater drought tolerance than the haplotype carrying deletion variants [24]. Another natural variant in the 5'-UTR of ZmPP2C-A10 deletes an ERSE (endoplasmic reticulum stress response element), the loss of this variant prohibits the expression of ZmPP2C-A10 during ER stress. Interestingly, disabling the ZmPP2C-A10 expression results in drought tolerance improvement in maize [25]. Moreover, some natural DREB (dehydration responsive element binding protein) variants in maize have also been reported, they are significantly associated with drought stress tolerance, which is also useful for improving drought tolerance in maize [26].
So far, natural variations associated with drought tolerance have only been reported in ZmNAC111 [21], variants in other NAC genes associated with drought tolerance have not been discovered. In the present study, we used maize RNA-Seq data to identify a gene encoding a NAC TF, which has relatively high homologies to other known drought response-related NACs. Association analysis revealed that the variants in ZmNAC080308 are significantly associated with GY in maize, especially under droughtstress conditions. Haplotype of ZmNAC080308 variants associated with drought stress tolerance were identified, A functional marker was developed and validated in an association mapping panel, which can be used for improving drought stress tolerance in maize.

Results
ZmNAC080308 is a nuclear-localized TF in the abiotic stress-responsive SNAC TF subfamily As compared with other NAC TFs, ZmNAC080308 contained the conserved NAC domain including five complete consensus NAC subdomains (A, B, C, D, E) ( Fig. 1 A). Phylogenetic analysis showed that genes with similar functions tended to cluster together, the stressresponsive SNAC TFs clustered into one branch with ATAF, while NAM and CUC, which were reported involving in regulation the plant growth and development, appeared in the other branches (Fig. 1B). ZmNAC080308 was more closely related to OsNAC3, ZmNAC55, ZmSNAC1, and SNAC1, according to our phylogenetic analysis. These four genes are associated with drought stress responses in plants [5,12,15,27], indicating that ZmNAC080308 might also be involved in drought stress responses in maize.
We subsequently transformed ZmNAC080308-GFP into maize protoplasts and detected strong GFP signals in the nuclei of the transformed cells, indicating nuclear localization of the ZmNAC080308-GFP fusion protein (Fig. 2 A). Next, we used two-hybrid transactivation assays to investigate whether the ZmNAC080308 protein could activate the expression of reporter genes AH109. The full-length CDS of ZmNAC080308 was cloned into the pGBKT7 vector, then the constructed or the empty vector was transformed into yeast strain AH109, respectively. As shown in Fig. 2B, both transformants grew well on the SD/-Trp medium, but only the strain containing pGBKT7-ZmNAC030803 grew on the SD/-Trp/-His medium and showed galactosidase activity. Therefore,  ZmNAC080308 transcription factor was a nuclearlocalized protein with transcriptional activity in maize.

Expression of ZmNAC080308 in tassel and seedling
Firstly, gene expression in young tassels harvested from ten inbred lines under well-watered (WW) and droughtstressed (DS) conditions were analyzed using RNA-Seq data (Wang et al., 2018). Interestingly, results showed that the transcript expression of ZmNAC080308 under DS conditions differed significantly from that under WW conditions in five maize lines (Tie7922, X178, Dan340, Ji81162 and CA339) at the thresholds of false discovery rate (FDR) ≤ 0.001, and an absolute value of the Log 2 Ratio ≥ 1. In the inbred lines of Tie7922, X178, and Dan340 showing better drought tolerance, the average abundance of ZmNAC080308 transcripts increased 2.14-fold after drought treatment, while it decreased 1.68-fold after drought treatment in inbred lines of Ji81162 and CA339 ( Table 1), showing that the differences in transcript expression could be guided by the signals of drought stress.
qRT-PCR (quantitative Reverse-Transcription Polymerase Chain Reaction) was used to analyze the relationship between ZmNAC080308 expression and drought stress tolerance in maize seedlings. ZmNAC080308 highly expressed in root and the lower portions of the stem (Fig. 3 A). After drought stress treatment, the transcript expression of ZmNAC080308 increased significantly in root, a peak was identified at 1 h after stress treatment. However, after exogenous ABA treatment, ZmNAC080308 transcript expression increased by more than 80-fold in stem ( Fig. 3B and C). This result presumes that ZmNAC080308 has a function regulating abiotic stresses in ABA pathway.

Genetic characterization of ZmNAC080308
We then sequenced the ZmNAC080308 gene in 199 NLs. A total of 3065 bp of sequence data was obtained from each line, including three exons, two introns, a 645-bp untranslated region and a 1232-bp promoter sequence (Table 2). After aligning the sequence of ZmNAC080308, 86 SNPs and 47 InDels were identified among these 199 inbred lines. Among all these SNPs and InDels, 72 variations were identified in the promoter region. Except for the promoter and the 5'-UTR regions, the exon III contained the highest SNP density, with an average of one SNP per 15 bp, while the intron I contained the most InDel density with an average of one InDel per 26 bp. No InDels were identified in the coding region. Among the 32 SNPs identified in exons, eight SNPs were non-synonymous mutations, and the other 24 SNPs were synonymous mutations. Average minor allele frequency (MAF) value was 0.22 L across all the 133 variations, with the lowest average MAF value of 0.28 in the promoter region and the highest average MAF value of 0.50 in the Exon I region.
Linkage disequilibrium (LD) analysis was performed on all the 133 variants using Haploview. Three LD blocks (r 2 > 0.6) were identified in the 5'-UTR region, part of Exon I and Intron I region, and the Exon III, respectively ( Fig. 4 C). Several variants in the promoter region exhibited strong LD, but most of other variants were not tightly linked with each other.
Variations in the`ZmNAC080308 gene are associated with GY under WW and DS conditions in maize A candidate gene association study was used to detect the variants in ZmNAC080308 that were significantly associated with drought stress tolerance in maize. Firstly, the association analysis between the variants and GY was performed, considering the importance of GY in evaluating drought stress tolerance. Using the Mixed Linear Model (MLM) method, 64 variations were detected that were significantly associated with GY in all ten environments. Among these 64 variations, only four variations in Intron I region were identified under WW condition, the rest of the 60 variations were identified only under DS conditions ( Fig. 4A and B). Seven variations had p values less than 0.01, including four markers in the promoter region and three markers in the 5'-UTR region, each marker explained an average of 6.10 % of the variation of GY. A total of 19 variations including five InDels and 14 SNPs were detected that were significantly associated with GY in at least two environments (p < 0.05) ( Fig. 4A and B), each variation explained an average of 4.45 % of the phenotypic variation of GY. All of the variations significantly associated with GY werein the non-coding regions of the ZmNAC080308 gene (five in the promoter, nine in the 5' UTR, four in Intron, and one in the 3' UTR). Strong LD (r 2 > 0.8) was observed among the variations located in the 5'-UTR, InDel-600, and SNP-586 in the promoter region, as well as between the four variations within the Intron I region (Fig. 4 C).
Variations in the 5'-UTR of ZmNAC080308 significantly associated with GY under the drought stress condition ZmNAC080308 was significantly associated with GY, especially under DS conditions, and the variations in the non-coding regions of the gene play an important role in regulating GY under drought stress conditions in maize. Among all the variations in ZmNAC080308 associated with GY under DS, nine variations in the 5'-UTR, including seven SNPs and two InDels, were in high LD (Table 3; Fig. 4). The 5'-UTR was important in regulating the gene expression, all of the inbred lines were divided into two haplotypes built with the nine variations in the 5'-UTR of ZmNAC080308. The two haplotypes had similar GY values under WW conditions, the GY of Hap2 was 421.17 g/plot under DS conditions, which was 15.47 % higher than that of Hap1 under the same conditions (p = 0.028) (Fig. 5). Then the cis-acting elements in the 5'-UTRs of the two haplotypes were analyzed to identify the differences in these elements between the variants. As expected, several cis-acting elements related to stress responses were present in the 5'-UTR region of ZmNAC080308 (Supplemental Fig. S1), including an ABRE (ABA-responsive element), a DRE (dehydrationresponsive element), a CGTCA sequence motif (involved in MeJA response), and a G-box (involved in light response) (Supplemental Fig. S1A). The presence of these motifs suggested that ZmNAC080308 might be regulated during drought stress responses. Moreover, an additional cis-acting element, Sp1 (light-responsive element), was caused by a nucleotide change from T to C at 42 bp upstream of the start codon of ZmNAC080308, which was identified in Hap2 genotypes (Supplemental Fig. S1B). Interestingly, relative luciferase activity was significantly higher under control of the 5' -UTR from Hap2 than that under the 5'-UTR from Hap1 in both controlled (p = 1.69E-07) and ABA-treated (p = 2.99E-06) conditions (Supplemental Fig. S1C). These results demonstrated that the mutations in the 5'-UTR of ZmNAC080308 might functionally regulate the gene expression between these two haplotypes.

Development of functional marker for prediction and selection of drought tolerant lines
To distinguish the two haplotypes of ZmNAC080308, a functional marker was developed based on the two InDel variations in the 5'-UTR region. This marker was then tested in another maize inbred line panel containing 187 US off-PVP inbred lines. These inbred lines were then classified into two haplotypes by screening all genotypes on SDS-PAGE (Fig. 6). In total, Hap1 existed in 116 inbred lines carrying the deletion variation, Hap2 existed in 70 inbred lines carrying the insertion variation. Hap2 had an average GY value of 63.22 g/plot under WS condition, which was higher than that of Hap1 under the same condition, but the differences were not significant (p > 0.05). However,, the average GY of Hap2 was   (Table 4). These results were consistent with those of the previous experiment, which proved that the variations in the 5'-UTR of ZmNAC080308 were related to maize drought resistance. The functional marker developed by the present study could be used to select the drought stress tolerance inbred lines.
Overexpression of ZmNAC080308 enhances drought resistance in transgenic Arabidopsis Three independent transgenic lines (OE1, OE3, and OE4) of ZmNAC080308 overexpression in Arabidopsis were chosen for further analysis. During the process of drought stress, it was observed that most wild-type (WT) plants withered at earlier stage, while the 35 S:: ZmNAC080308 transgenic plants remained in green condition (Fig. 7 A). After withholding water for 3 weeks, most leaves of both WT plants and the transgenic plants began to dry for rehydration. Further statistical analysis of the survival rate of the Arabidopsis after drought stress revealed that~50 % of transgenic plants were still alive, whereas~85 % of WT plants died (Fig. 7B). These results showed that overexpression of ZmNAC080308 can improve drought tolerance of the transgenic Arabidopsis.

Discussion
ZmNAC080308 encodes a NAC transcription factor affecting drought stress tolerance in maize NAC TFs and their SNAC subfamily comprise an important group of genes involved in plant abiotic stress responses. Previous phylogenetic analysis showed that, 16 NAC genes belong to the SNAC subfamily in maize [6]. Among these 16 SNAC genes, three NAC genes including ZmNAC010312 (GRMZM2G347043, ZmSNAC1) [27], ZmNAC050436 (GRMZM2G336533, ZmNAC55) [12], and ZmNAC100475 (GRMZM2G127379, ZmNAC111)  [21] have also been reported responding to drought stress tolerance in maize. In the present study, we cloned another maize SNAC subfamily gene of ZmNAC080308 based on RNA-Seq data analysis, this gene presumed to respond to drought stress resistance. After DS treatment, the transcript abundance of ZmNAC080308 was higher in the young maize tassels of inbred lines showing moderate to high drought stress,, while it was relatively low in the susceptible inbred lines ( Table 1). The 5'-UTR of ZmNAC080308 was analyzed, the stress responsive cis-elements, such as ABRE and DRE, were identified (Supplemental Fig. S1). Further, just like that of other SNAC genes such as ZmNAC55 [12], OsNAC3 [5], and ZmSNAC1 [27]. qRT-PCR analysis also showed that expression of ZmNAC080308 could be induced by both ABA and drought stresses in roots and stems (Fig. 3), However, the expression of ZmNAC080308 under drought and ABA stresses differed between roots and stems. In roots, the expression of ZmNAC080308 was more sensitive to drought, whereas the expression of ZmNAC080308 rose more rapidly in stems after ABA treatment (Fig. 3). Results show that the expression of the ZmNAC080308 gene presumed responding to drought stress through both ABA-dependent pathway in stem and ABA-independent pathway in root, which was consistent with the results of our previous study, the ZmNAC080308 gene contains both ABRE and DRE in its promoter [6].
Overexpression of ZmNAC080308 enhances drought tolerance in transgenic Arabidopsis. The survival rate of the transgenic Arabidopsis was not as high as expected, which might cause by the severe drought stress, although it was still much higher than the survival rate of WT. All the results in the present study indicate that ZmNAC080308 variants play an important role responding to drought stress through both ABA-dependent and ABAindependent pathways in maize.
A total of 32 SNPs in the ZmNAC080308 coding region were identified, eight of them were nonsynonymous mutations, which located in Exon III region. Association analysis showed that two of the eight variants were significantly associated with GY evaluated under DS conditions. However, the effects of both variants were only identified in one of five tested environments (location × year) i.e. SNP-118 was detected in 2008 in Xinjiang and SNP120 was detected in 2009 in Xinjiang, this result indicates that the effects of these two associations were not stable. Moreover, TFs can be involved in multiple pathways, the variants identified in  the present study might also be related to other traits, such as flowering time. Therefore, further analyses should be conducted to detect the associations between ZmNAC080308 and other important agronomic traits.
Variations in the non-coding regions of ZmNAC080308 affecting the drought stress tolerance Candidate gene-based association analysis was used to detect the associations between genetic polymorphisms and the target traits. Compared to genome-wide association (GWAS) strategies, the candidate gene approach requires fewer genetic markers, it can minimize the statistical issues caused by multiple-testing across the entire genome [28]. Candidate gene association analysis was firstly used to analyze the association between the dwarf8 gene and flowering time in maize [29], where nine flowering time-associated polymorphisms were identified. In a previous study, association analysis showed that a deletion variant in the 5'-UTR of ZmPP2C-A10 decreased the gene expressionlevel, which resulted in maize drought tolerance improvement [25]. Another drought tolerance-related gene, ZmVPP1, was also identified by GWAS and candidate gene association analysis, two InDel variants on the upstream of its start codon and two variants within its CDS were detected, which were significantly associated with drought tolerance in maize. Furthermore, the drought tolerancerelated haplotype of ZmVPP1 was identified [24]. An insertion variant in the promoter region of ZmNAC111 was also associated with drought tolerance in maize. This variant correlated with lower expression of the ZmNAC111 gene, which affected DNA and histone methylation and caused maize seedlings to be more sensitive to drought stress [21]. In the present study, nine highly linked polymorphisms including two InDel variants in the 5'-UTR of ZmNAC080308 were identified and associated with GY under drought stress conditions. Two haplotypes of ZmNAC080308 were built using these variations, Hap2 exhibited higher GY than Hap1 under drought stress conditions. In addition, a functional marker in the 5'-UTR of ZmNAC080308 (Fig. 6) was developed based on the InDel variants and an off-PVP maize inbred line panel was genotyped with this functional marker. Results showed that the GY of the inbred lines carrying the favorable genotypes of this functional marker was 11.30 % higher, indicateing that this functional marker could be used in MAS for improving droughttolerance in temperate maize. The 5'-UTR is an important region that regulates the expression of eukaryotic genes (Barrett et al., 2012;Cenik et al., 2011). This study detected a SNP located at 42 bp upstream of the start codon of ZmNAC080308, which creates the cis-acting element Sp1 in Hap2 genotypes (Supplemental Fig. S1B). A previous study showed that promoters from Sp1 mutant lines had reduced GUS expression compared with promoters isolated from wildtype lines [30]. Moreover, the Sp1 cis-element is one of the cis-elements that are presented in the promoter region of expressed drought responsive genes [31]. Therefore, we hypothesized that this kind of mutation in ZmNAC080308 could modulate drought-stress tolerance in maize. Further, the 5'-UTR of the ZmNAC080308 Hap2 genotype showed higher transcriptional activity than that of the Hap1 genotype in the dual-luciferase reporter assay system (Supplemental Fig. S1D). Therefore, we tested the expression level of the ZmNAC080308 gene in different haplotypes using RNA-Seq analysis. Among the five lines we tested, only X178 carries Hap2, the other four lines carry Hap1 ( Table 1). The RPKM (Reads Per Kilobase per Million mapped reads) value of ZmNAC080308 in X178 was higher than that of the other four lines under both WW and DS conditions, indicating that the addition of the Sp1 motif can increase the expression of ZmNAC080308, which is related with drought stress tolerance in maize. Finally, the drought tolerance ability of the 35 S::ZmNAC080308 were tested in transgenic Arabidopsis. The survival rate of the transgenic plants was much higher than that of the WT plants (Fig. 7). These results confirm that ZmNAC080308 functions as a positive regulator of drought stress, it can be used for maintaining the GY under drought stress conditions in maize.

Conclusions
ZmNAC080308 is an important gene responding to drought tolerance and its marker from variations could be used as a functional identifier for improving maize GY under drought stress.

Collection of GY data under well-watered (WW) and drought-stressed (DS) conditions
In the present study, GY data for 389 lines in two association panels were collected and analyzed. More detailed information regarding these data can be found in our previous studies [32,33].
The other association panel was composed of 186 US off-PVP inbred lines (ALs; US lines no longer under patent). This panel was planted in a α-lattice design in Hainan and Xinjiang, the same locations as the NLs in 2013 with two watering treatments (WW and DS conditions). The experiment design and irrigation methods were the same as those used for the NLs, which have been described in our previous study [32]. In brief, each line was planted two replications in single-row 4-m long plot with 20 cm between plants, and 21 plants per row. The GY data of ALs under each treatment was collected three weeks after harvest.

Patterns of mRNA expression of ZmNAC080308
The expression levels of ZmNAC080308 in young tassels under normal and drought-stressed conditions were from a previous RNA-Seq study [34]. Expression data of five inbred lines including Ji81162, CA339, Dan340, X178, and Tie7922, exhibited different degrees of drought stress tolerance (Table 1).
In the present study, qRT-PCR was conducted to characterize the expression of ZmNAC080308 after drought or ABA treatment at seedling stages. Seeds of inbred line Tie7922 were first planted in quartz sand in the greenhouse and watered every other day. Endosperms were removed at the V1 stage, and then seedlings of similar sizes were transferred into buckets filled with Hoagland's nutrient solution for continued growth. The solution was refreshed every 2 d until the V2 stage. The seedlings were then transferred into Hoagland's solution containing 20 % PEG-6000 for the drought treatment or 100 µM ABA for the ABA treatment. The roots, leaves, and stems of each treated material and controls were sampled separately after treatment at 0 h, 1 h, 3 h, 6 h, and 12 h. Five samples from each treatment at the same sample collection were mixed equally as a replication, and two replications were used in this study.
Total RNA was isolated from samples using TransZol Up (Transgen, China) according to the manufacturer's instructions. A FastQuant RT Kit (Tiangen, China) was used to reverse transcribe the first-strand cDNA from 1000 ng of total RNA. qRT-PCR was then performed using a SuperReal PreMix Plus Kit (SYBR Green) (Tiangen, China) with IQ5 (Bio-Rad, US). The relative expression level of ZmNAC080308 was calculated relative to the expression level of the maize NADPH gene using the 2 −ΔΔCt method [35].

Sequencing of full-length ZmNAC080308 in NLs
The seeds of all NLs were planted in three 72-cell germination boxes in the greenhouse. One leaf was sampled from each line at the V2 stage, and genomic DNA was extracted from each line using the CTAB method. The gene sequence of ZmNAC080308 was downloaded from the National Center for Biotechnology Information database (NCBI, https://www.ncbi.nlm.nih.gov/). The gene structure of ZmNAC080308 was predicted using the online tool FGENESH (http://www.softberry.com) and its 5'-UTR was searched for cis-acting regulatory elements using the PlantCARE database [36]. We then cloned the full-length ZmNAC080308 genomic sequence and about 1.5 kb upstream before start site as its promoter region from each of the 199 NLs by PCR using Prime STAR GXL DNA Polymerase (Takara, China) and four pairs of overlapping primers (Supplementary Table 1). The amplified gene fragments were sequenced and assembled by Tsingke Biological Technology (China).
Multiple alignments of DNA sequences and amino acid sequences were performed separately using DNAM AN and MEGA 7.0 software [37]. SNPs and InDels with minor allele frequency greater than 0.05 were identified visually. Phylogenetic analysis of ZmNAC080308 with other NAC proteins was performed by constructing a tree in MEGA 7.0 using the neighbor-joining (NJ) method with a bootstrap of 1000 replications. Evolutionary distances were computed using the p-distance method and all ambiguous positions were removed from each sequence pair [38,39].

Linkage disequilibrium (LD) and association analysis of ZmNAC080308
After identifying variants in the ZmNAC080308 gene, LD analysis was performed using Haploview software [40]. Candidate gene association analysis was conducted using a Mixed Linear Model in Tassel 5.0 [41]. Population structure (Q) and kinship (K) data from previous studies [33,42] were used to improve the fidelity of GYassociated SNP detection [43]. Markers with p-value < 0.05 (-log 10 (p) = 1.30) were considered to be significantly associated with GY.

Subcellular localization of the ZmNAC080308 protein
The vector pRTL2 was digested with the restriction enzymes SacI and KpnI (NEB, USA) for > 3 h. The coding region of ZmNAC080308 was cloned from Qi319, an NL inbred line, without the terminator codon (TGA) by PCR to include the corresponding 15-bp extensions homologous to the vector ends (sense: 5'-CCGGGAAT TCCATGG-3'; antisense: 5'-CTCTAGAATGGTGAG-3'). The digested vector and the PCR product were recovered using a Zymoclean™ Gel DNA Recovery Kit (Zymo Research, USA). We then constructed the recombinant plasmid using an In-Fusion HD Cloning Kit (Takara, China) according to the manufacturer's instructions. The recombinant plasmid and the empty vector plasmid were transformed into maize leaf protoplasts using the PEG method [44]. After a 16-h incubation in darkness, GFP fluorescence was observed using an LSM880 fluorescence microscope (Zeiss, Germany).

Transactivation assay
To analyze whether ZmNAC080308 possessed transactivation activity, the yeast strain AH109 containing HIS3 and lacZ reporter genes was used to express the ZmNAC080308 gene fused to the GAL4 DNA-binding domain. The fulllength coding sequence of ZmNAC080308 was obtained by PCR using primers: 5'-ATGGCCATGGAGGCCGAATTCA TGGCAATGGTGGCGGCG-3' and 5'-ATGCGGCCGC TGCAGGTCGACTCAGAACGGAGGCAAGAT-3'. The resulting PCR product was cloned downstream of the GAL4 binding domain in the pGBTK7 bait vector digested with SalI (NEB, USA). The recombinant plasmid was then transformed into yeast strain AH109 using Matchmaker™ Gold Yeast Two-Hybrid System (Takara, China) according to the manufacturer's protocol, the empty vector was separately transformed into AH109 as a negative control. The transformed strains were diluted 1:10 0 , 1:10 1 , 1:10 2 , 1:10 3 . The dilutions were plated on SD/-Trp and SD/-Trp/-His (x-gal) solid medium, respectively.
To compare the transcriptional activation abilities of the the two ZmNAC080308 haplotypes at the 5'-UTR region, the Dual-Luciferase Reporter Assay System (Promega, USA) was applied. The maize ubiquitin (Ubi) promoter and the two ZmNAC080308 5'-UTR haplotypes were inserted into the Dual-Luciferase Reporter vector pGreen II 0800, this vector contained a Renilla luciferase (REN) gene under the constitutive CaMV 35 S promoter as a control (Supplemental Fig. S1C). The three different recombinant plasmids were cotransformed into maize protoplasts separated using a PEG-mediated method as mentioned above. For ABA treatments, the transfected protoplasts were first incubated in solution without ABA [45]. After a 16-h incubation, LUC (luciferase) and REN activities were measured using a fluorescence microscope (Zeiss, Germany).

Transgenic Arabidopsis plants construction and stress treatment
The coding region of ZmNAC080308 was amplified from the maize inbred line of Tie7922 by PCR. The PCR products were then inserted into the pCAMBIA1303 vector under the CaMV 35 S promoter. The reconstructed vector was introduced into Agrobacterium tumefaciens train after sequencing, and then it was transferred into Arabidopsis by floral dip method. The T0 generation seeds of transgenic plants were screened on 1/2 Murashige and Skoog (MS) medium with kanamycin. Then the gained positive plants were transplanted to pots filled with a 2:1 mixture of vermiculite and nutrition soil.
For the drought tolerance assay, the T2 plants were transferred from 1/2 MS medium into pots filled with a 2:1 mixture of vermiculite and nutrition soil at the twoleaves stage. After three weeks of recovery and growth under the well-watered condition, the plants were started to be treated under drought stress by withholding watering for three weeks. Then the plants were recovered by rehydration. Three days later, the survival rate of the plants was recorded.