Dissection of Maize Drought Tolerance at the Flowering Stage Using Genome-Wide Association Studies

Drought is one of the most critical environmental factors constraining maize production. When it occurs at the flowering stage, serious yield losses are caused, and often, the damage is irretrievable. In this study, anthesis to silk interval (ASI), plant height (PH), and ear biomass at the silking date (EBM) of 279 inbred lines were studied under both water-stress (WS) and well-water (WW) field conditions, for three consecutive years. Averagely, ASI was extended by 25.96%, EBM was decreased by 17.54%, and the PH was reduced by 12.47% under drought stress. Genome-wide association studies were carried out using phenotypic values under WS, WW, and drought-tolerance index (WS-WW or WS/WW) and applying a mixed linear model that controls both population structure and relative kinship. In total, 71, 159, and 21 SNPs, located in 32, 59, and 12 genes, were significantly (P < 10−5) associated with ASI, EBM, and PH, respectively. Only a few overlapped candidate genes were found to be associated with the same drought-related traits under different environments, for example, ARABIDILLO 1, glycoprotein, Tic22-like, and zinc-finger family protein for ASI; 26S proteasome non-ATPase and pyridoxal phosphate transferase for EBM; 11-ß-hydroxysteroid dehydrogenase, uncharacterised, Leu-rich repeat protein kinase, and SF16 protein for PH. Furthermore, most candidate genes were revealed to be drought-responsive in an association panel. Meanwhile, the favourable alleles/key variations were identified with a haplotype analysis. These candidate genes and their key variations provide insight into the genetic basis of drought tolerance, especially for the female inflorescence, and will facilitate drought-tolerant maize breeding.


Introduction
Maize (Zea mays L.) is the most widely grown crop worldwide and has incredible importance for food, feed, and other industrial products [1]. It was estimated that maize production needs to be boosted by 67% to cope with the increased population growth and food demand in 2050 [2]. Thus, it is of great significance to increase maize yield productivity and reduce yield loss caused by biotic and abiotic stresses.
Drought is considered one of the most detrimental curbs of agriculture, which endangers maize production globally because of its erratic nature [3,4]. Maize is an openpollinating crop that is extremely sensitive to drought stress throughout its life span, especially one week before and three weeks after the flowering period, causing severe yield loss by 15-25% [5]. At this stage, drought inhibits ear development, causing abnormal differentiation in spikelets, and changes in hormone signalling involved in cell division, growth, and primordium development [6], leading to asynchronous development in tassel and ear, a tremendous extension of anthesis and silking interval (ASI), and reduction in silk This panel was obtained from an available association population for eQTL mapping on maize kernel development [25,26].
The field trials were conducted in Urumqi (Xinjiang province: 43 • 54 N, 87 • 28 E), from May to September. The flowering date of the association panel was recorded in 2016 under normal growing conditions and then divided into three patches according to growing degree day (GDD) value at −15 days before pollen shedding (−15 D). The daily temperature was recorded, and GDD was calculated for each line using the following formula: where L means daily lowest temperature ( • F), and H means daily highest temperature ( • F).
In order to obtain synchronised flowering, 279 lines were delayed sowing in patches. Each patch was in randomised complete block design, with three independent repeats during the years of 2017, 2018, and 2020. A total of 119 lines belonged to the 1st patch (No. 1-119, Table S1), with GDD ranging from 1612.  Table S1) with earlier flowering times (GDD = 962-1291. 6) were planted on 28 May 2017. Each row was 3.6 m in length, with planting space of 0.24 m and 1.1 m width across rows. Every line was planted side by side, one row for water stress treatment (WS) and the adjacent for a well-watered regime (WW, Figure 1A), with an independent valve to control the irrigation of each row. Drought stress was applied in the WS regime under the Protocol of Irrigation Management for Maize Drought Trials of Syngenta. Briefly, daily moisture level was recorded using the soil moisture sensors WATERMARK MONITOR 900 M (IRROMETER) installed in the field ( Figure 1B) and connected to the sensor in every patch at 60 cm and 90 cm depth in multiple locations according to instructions from the trial sponsor in both WS and WW regimes. For the intended drought stress management at the flowering stage, a drought stress trial was conducted by withdrawing irrigation initiated at the −21 D before anthesis estimated according to GDD. The soil moisture was maintained at <40 centibars before −15 D and at 80-120 centibars from −15 D to −7 D. At −7 D, the soil moisture reached 120-150 centibars, which is the trigger point for drought stress, and was kept at 150-200 centibars until 14 days after anthesis when irrigation was resumed. The fertilisers, herbicides, and insecticides were applied as per requirements according to the local recommendation practices. In the WW regime, normal irrigation was provided once a week.

Plant Materials and Experimental Design
An association panel for GWAS comprising 279 inbred lines with diverse genetic variation (Table S1) was selected from comprehensive geographical regions including tropical, subtropical, and temperate areas of the world, as described by Zhang et al. [24]. This panel was obtained from an available association population for eQTL mapping on maize kernel development [25,26].
The field trials were conducted in Urumqi (Xinjiang province: 43°54′ N, 87°28′ E), from May to September. The flowering date of the association panel was recorded in 2016 under normal growing conditions and then divided into three patches according to growing degree day (GDD) value at −15 days before pollen shedding (−15 D). The daily temperature was recorded, and GDD was calculated for each line using the following formula:

50
where L means daily lowest temperature (°F), and H means daily highest temperature (°F). In order to obtain synchronised flowering, 279 lines were delayed sowing in patches. Each patch was in randomised complete block design, with three independent repeats during the years of 2017, 2018, and 2020. A total of 119 lines belonged to the 1st patch (No. 1-119,  Table S1) with earlier flowering times (GDD = 962-1291.6) were planted on 28 May 2017. Each row was 3.6 m in length, with planting space of 0.24 m and 1.1 m width across rows. Every line was planted side by side, one row for water stress treatment (WS) and the adjacent for a well-watered regime (WW, Figure 1A), with an independent valve to control the irrigation of each row. Drought stress was applied in the WS regime under the Protocol of Irrigation Management for Maize Drought Trials of Syngenta. Briefly, daily moisture level was recorded using the soil moisture sensors WATER-MARK MONITOR 900 M (IRROMETER) installed in the field ( Figure 1B) and connected to the sensor in every patch at 60 cm and 90 cm depth in multiple locations according to instructions from the trial sponsor in both WS and WW regimes. For the intended drought stress management at the flowering stage, a drought stress trial was conducted by withdrawing irrigation initiated at the −21 D before anthesis estimated according to GDD. The soil moisture was maintained at <40 centibars before −15 D and at 80-120 centibars from −15 D to −7 D. At −7 D, the soil moisture reached 120-150 centibars, which is the trigger point for drought stress, and was kept at 150-200 centibars until 14 days after anthesis when irrigation was resumed. The fertilisers, herbicides, and insecticides were applied as per requirements according to the local recommendation practices. In the WW regime, normal irrigation was provided once a week.  (SD) were determined by the number of days from planting to pollen shedding of 50% of plants and 50% of plants having clearly visible silks, respectively. ASI was calculated as the interval between anthesis date and silking date (ASI = SD − AD) and counted in days [9]. PH was measured from the soil level to the lowest tassel branch of each plant and recorded in centimetres (cm). For ear biomass, the ear at its silking date ( Figure 1C) of each line was harvested and dried to a consistent weight at 72 • C. The ear biomasses were weighed with an electronic weighing balance and recorded in grams (g). The drought-tolerant index was calculated by dividing the mean values of studied traits in the WS regime by those of the WW regime. To estimate the random errors, each measurement per inbred line comprised 6-8 individual plants with 3 independent repeats.

Association Analysis for ASI, EBM, and PH
For association analysis, the 279 lines were genotyped by 776,254 high-density SNP markers, with MAF > 0.05, which was excavated previously [26]. The genome-wide association analysis was conducted with a mixed linear model (MLM) embedded in TASSEL V5.0 [27]. Population structure (Q) and kinship (K) were estimated according to Pang et al. [26]. Briefly, kinship matrices and principal components were estimated based on 236,205 SNPs with MAF > 0.05. The top 3 principal components were selected as population structure, and for kinship estimation, the 'Normalised_IBS' method was used. The regression-based coefficient of determination values of all significantly associated SNPs was recorded to determine the variations explained by each SNP locus. Drought tolerance was controlled by several genes with minor effects and non-independence of SNPs in maize genome leading to strong LD. Therefore, a threshold of P < 10 −5 was applied in this study to eliminate false negatives. Multiple Manhattan plots were drawn in R (https://www.r-project.org/ (accessed on 8 November 2021)).

Drought Responsive, Linkage Disequilibrium, and Haplotype Analysis of Candidate Genes
In total, 197 inbred lines from a previous association mapping population were planted in pots under well-watered and drought stress treatments [6]. Drought stress conditions were relative to WW conditions, with the soil moisture levels of WS conditions dropping from 50% to 10% and soil moisture levels of WW conditions keeping values around 50% at the flowering stage. Candidate genes-Zm00001d013992, Zm00001d020506, Zm00001d029937, Zm00001d029938, and Zm00001d039319-were commonly identified under two environments. Therefore, to determine whether they were drought-responsive, their expression patterns were analysed according to the RNA profiling data of the abovementioned 197 lines (unpublished data from Dr. Mingqiu Dai, Huazhong Agricultural University). Haplotype analysis normally helps to understand the causal variation in a gene identified by GWAS [20,26], which in turn benefits future functional marker development. To confirm whether the SNPs identified for each gene were meaningful, haplotype analysis was conducted to check the allele effects of the most significant SNPs for the six overlapping candidate genes. A two-pair t-test was used to analyse their allelic effect in corresponding phenotypic performance.

Statistical Analysis
The phenotypic data were filtered by removing the suspicious value of each replicate by Q-test in excel [24]. The mean value of three replicates in each environment was used for association analysis. Best linear unbiased estimators (BLUEs) were estimated using the genotype and covariate as a fixed factor, and the rest as random factors. Analysis of variance (ANOVA), correlations analysis, and broad-sense heritability (H 2 ) were estimated by SPSS 25.0. The H 2 was estimated according to the following equation: where σ 2 G represents genetic variance, σ 2 GE shows the genetic and environmental interaction variation, σ 2 e shows residual error variance, n shows the number of environments, and r shows the number of replicates [28]. LD analysis of these candidate genes was conducted by SNPs within a gene using the 'LD heatmap' in the R software (https://www.r-project.org/ (accessed on 8 November 2021)).

Performance of Drought-Tolerant Phenotypes in the Association Panel
The test of normality revealed that the frequency distributions for all of the studied traits were near normal for most of the traits in the association panel ( Figure 2; Table 1). The existence of variation might be caused by fluctuating environmental conditions in the field or genotypic differences. The descriptive statistics, heritability analysis, and coefficient of variation (CV) for the phenotypic traits are listed in Table 1. Generally, drought stress significantly decreased the ear biomass and PH ( Figure 1B,C) while enlarging the ASI, suggesting that water stress at the flowering stage had diverse effects on drought-related traits.
for the six overlapping candidate genes. A two-pair t-test was used to analyse their allelic effect in corresponding phenotypic performance.

Statistical Analysis
The phenotypic data were filtered by removing the suspicious value of each replicate by Q-test in excel [24]. The mean value of three replicates in each environment was used for association analysis. Best linear unbiased estimators (BLUEs) were estimated using the genotype and covariate as a fixed factor, and the rest as random factors. Analysis of variance (ANOVA), correlations analysis, and broad-sense heritability (H 2 ) were estimated by SPSS 25.0. The H 2 was estimated according to the following equation: where σ 2 G represents genetic variance, σ 2 GE shows the genetic and environmental interaction variation, σ 2 e shows residual error variance, n shows the number of environments, and r shows the number of replicates [28]. LD analysis of these candidate genes was conducted by SNPs within a gene using the 'LD heatmap' in the R software (https://www.rproject.org/ (accessed on 08, 11, 2021)).

Performance of Drought-Tolerant Phenotypes in the Association Panel
The test of normality revealed that the frequency distributions for all of the studied traits were near normal for most of the traits in the association panel ( Figure 2; Table 1). The existence of variation might be caused by fluctuating environmental conditions in the field or genotypic differences. The descriptive statistics, heritability analysis, and coefficient of variation (CV) for the phenotypic traits are listed in Table 1. Generally, drought stress significantly decreased the ear biomass and PH ( Figure 1B,C) while enlarging the ASI, suggesting that water stress at the flowering stage had diverse effects on droughtrelated traits.     (Table 1). H 2 values of ASI, EBM, and PH were higher than 80% under the two water treatments across the three years (Table 1), indicating that these traits are highly heritable in the two water treatment conditions.

Correlations among Drought-Related Traits
The correlations among ASI, EBM, and PH values under WS and WW conditions are listed in Table 2. Significantly positive correlations (P < 0.05) were found between ASI and PH under WS (0.13) and WW (0.15) conditions. Additionally, significant correlations were found between ear biomass and PH under both WS (0.22) and WW (0.18) conditions. A weak and negative correlation was found between ASI and ear biomass.

GWAS for Maize Drought Tolerance Genes
A total of 71, 159, and 21 SNPs were significantly associated with ASI, EBM, and PH located in 36, 81, and 16 genes, respectively. These genes were scattered over 10 chromosomes, with R 2 ranging from 7.19% to 15.52%. The number of identified SNPs in each chromosome ranged from 9 to 31 on chromosomes 9 and 1, respectively. The information about the SNPs and candidate genes is shown in Tables S2-S4. A total of 17 and 48 SNPs located in 13 and 18 genes were identified for ASI-WS and ASI-WW, with R 2 ranging from 7.73% to 14.32% (Figure 3; Tables 3 and S2). As to the ASI delay, only five SNPs in four genes were identified in 2017 and 2018 (Table S2). In total, 49 and 93 SNPs located in 23 and 43 were associated with EBM-WS and EBM-WW, with R 2 ranging from 7.99% to 12.52% (Figure 4; Tables 3 and  S3). In addition, 17 SNPs in 15 genes were associated with the EBM drought tolerance index (WS/WW) across the three years (Table S3). For both conditions, nine SNPs were associated with PH-WS and PH-WW and located in seven and six genes, with R 2 ranging from 7.19% to 15.52% across the three-year environments ( Figure 5; Tables 3 and S4). Three genes in Chr. 2, 7, and 9 were identified for the PH drought tolerance index (Table S4).

Common Genes Identified for Ear Development across Multiple Years or Conditions
Several overlapping genes were identified among years and water treatments (Table 3). Four candidate genes encoding ARABIDILLO 1 protein (Zm00001d029938) and Glycoprotein (Zm00001d029937) on Chr. 1, Tic22-like family protein (Zm00001d039319), and zinc-finger family protein (Zm00001d042997) on Chr. 3, were identified for ASI under both WS and WW conditions in 2018 (Figure 3), suggesting that these genes might be promising candidates, as they can function in maize inflorescence development with or without water deficit. An SNP (S5_27121944), located in pyridoxal phosphate transferase encoding gene (Zm00001d013992), was significantly associated (P < 10 −5 ) with EBM under drought stress condition on Chr. 5 and was consistently detected for two years environments of 2018 and 2020 (Figure 4), suggesting the drought tolerance role of this gene in developing ear. A candidate gene (Zm00001d020506) encoding 26S proteasome non-ATPase regulatory subunit-9 was detected for EBM on Chr. 7 under both WS and WW regimes during the 2017 field trial (Figure 4), suggesting its possible role in maize ear development.

Common Genes Identified for Ear Development across Multiple Years or Conditions
Several overlapping genes were identified among years and water treatments ( Table  3). Four candidate genes encoding ARABIDILLO 1 protein (Zm00001d029938) and Glycoprotein (Zm00001d029937) on Chr. 1, Tic22-like family protein (Zm00001d039319), and zinc-finger family protein (Zm00001d042997) on Chr. 3, were identified for ASI under both WS and WW conditions in 2018 (Figure 3), suggesting that these genes might be promising candidates, as they can function in maize inflorescence development with or without water deficit. An SNP (S5_27121944), located in pyridoxal phosphate transferase encoding gene (Zm00001d013992), was significantly associated (P < 10 −5 ) with EBM under drought stress condition on Chr. 5 and was consistently detected for two years environments of 2018 and 2020 (Figure 4), suggesting the drought tolerance role of this gene in developing ear. A candidate gene (Zm00001d020506) encoding 26S proteasome non-ATPase regulatory subunit-9 was detected for EBM on Chr. 7 under both WS and WW regimes during the 2017 field trial (Figure 4), suggesting its possible role in maize ear development.

Candidate Genes Drought Responsive Pattern
Most quantitative traits' functional genes are responsive at the transcriptional level. Drought tolerance is complex and regulated by many quantitative trait loci (QTLs) with

Candidate Genes Drought Responsive Pattern
Most quantitative traits' functional genes are responsive at the transcriptional level. Drought tolerance is complex and regulated by many quantitative trait loci (QTLs) with minor effects. In order to determine whether these candidate genes were drought-regulated, we analysed their expression level using the expression data from 197 diverse inbred lines under WS and WW conditions. Unsurprisingly, significant differential expression existed in most candidate genes between WW and WS treatments. Zm00001d013992, Zm00001d029938, and Zm00001d039319 were increased by 45.44%, 17.46%, and 6.01%, while Zm00001d029937 was reduced by 30.26%, compared with their expression under the WW condition ( Figure 6). However, no obvious difference was observed for Zm00001d039319. Genes 2022, 12, x FOR PEER REVIEW 10 of 15

Allele Effects of Common Candidate Genes
SNPs within a gene might be critical to its function. Therefore, it was of great interest to find the true association between SNP variation and target traits in a population. Haplotype analysis is a useful strategy to extract more rare causal variants. In our study, inbred lines were grouped into HapA (A) and HapB (G) based on the most significant SNP (S5_27121944) variation in Zm00001d013992 (Table 3, Figure 7A), 226 (2018), and 213 (2020). The rare HapB group had significantly higher EBM (2.10 g) under WS conditions in both 2018 and 2020 ( Figure 7B, Table S5) than HapA. Five polymorphisms-namely, chr7.S_116288756, chr7.S_116288791, chr7.S_116288792, chr7.S_116285652, and chr7.S_116285655-were identified in Zm00001d020506 ( Figure 7C, Table S6). In total, 101 lines belonged to the HapA (AACCT) with average EBM values of 1.58 g and 1.79 g under WS and WW in 2017, while the other 107 lines carried the GTTTC haplotype with significantly lower EBM values of 1.35 g and 1.55 g under WS and WW ( Figure 7D). Three SNPs in Zm00001d029937 were significantly associated with ASI under the WW condition in 2018 (Table 3, Figure 7E), and accordingly, two haplotype AAC (HapA) and GGG (HapB) were obtained. Of the 193 lines harbouring rare haplotype HapA, 15 showed longer ASI (>10 days) than those harbouring HapB, whose ASI values were 5.6 days and 6.6 days under WW and WS conditions, respectively ( Figure 7F, Table S7). The same 15 lines belonging to the HapA group (TCGATAATC) presented longer ASI under both WW and WS conditions ( Figure 7G,H, Table S8), which was a similar case as that of the adjacent gene Zm00001d029938. PZE-103003226 was the only SNP identified in Zm00001d039319 ( Figure  7I). The 165 lines with shorter ASI, grouped into HapA (G), showed weaker sensitivity to drought stress in both WW and WS conditions ( Figure 7J, Table S9). For Zm00001d042997, 9 SNPs classified the 237 inbred lines into 3 groups ( Figure 7K, Table S10), among which 17 lines belonged to lower frequency of HapC (TGACTTAA) and had a larger ASI under both WS and WW conditions, compared with those belonging to HapA (AACTCCCG) and HapB (AACTCTAA), while no statistical difference was found between HapA and HapB ( Figure 7L).

Allele Effects of Common Candidate Genes
SNPs within a gene might be critical to its function. Therefore, it was of great interest to find the true association between SNP variation and target traits in a population. Haplotype analysis is a useful strategy to extract more rare causal variants. In our study, inbred lines were grouped into HapA (A) and HapB (G) based on the most significant SNP (S5_27121944) variation in Zm00001d013992 (Table 3, Figure 7A), 226 (2018), and 213 (2020). The rare HapB group had significantly higher EBM (2.10 g) under WS conditions in both 2018 and 2020 ( Figure 7B, Table S5) than HapA. Five polymorphisms-namely, chr7.S_116288756, chr7.S_116288791, chr7.S_116288792, chr7.S_116285652, and chr7.S_116285655-were identified in Zm00001d020506 ( Figure 7C, Table S6). In total, 101 lines belonged to the HapA (AACCT) with average EBM values of 1.58 g and 1.79 g under WS and WW in 2017, while the other 107 lines carried the GTTTC haplotype with significantly lower EBM values of 1.35 g and 1.55 g under WS and WW ( Figure 7D). Three SNPs in Zm00001d029937 were significantly associated with ASI under the WW condition in 2018 (Table 3, Figure 7E), and accordingly, two haplotype AAC (HapA) and GGG (HapB) were obtained. Of the 193 lines harbouring rare haplotype HapA, 15 showed longer ASI (>10 days) than those harbouring HapB, whose ASI values were 5.6 days and 6.6 days under WW and WS conditions, respectively ( Figure 7F, Table S7). The same 15 lines belonging to the HapA group (TCGATAATC) presented longer ASI under both WW and WS conditions ( Figure 7G,H, Table S8), which was a similar case as that of the adjacent gene Zm00001d029938. PZE-103003226 was the only SNP identified in Zm00001d039319 ( Figure 7I). The 165 lines with shorter ASI, grouped into HapA (G), showed weaker sensitivity to drought stress in both WW and WS conditions ( Figure 7J, Table S9). For Zm00001d042997, 9 SNPs classified the 237 inbred lines into 3 groups ( Figure 7K, Table S10), among which 17 lines belonged to lower frequency of HapC (TGACTTAA) and had a larger ASI under both WS and WW conditions, compared with those belonging to HapA (AACTCCCG) and HapB (AACTCTAA), while no statistical difference was found between HapA and HapB ( Figure 7L).

Discussion
Maize is highly sensitive to drought stress at the flowering stage. In this study, the silk extrusion was markedly delayed, and plant architecture was greatly affected under drought stress. These results are in line with previous reports that stated significant extension in ASI and a notable reduction in plant and ear height under WS conditions [29,30]. Water deficit at the flowering stage delays or inhibits plant growth and ear development, reducing ear biomass from 1.71 g to 1.41 g (Table 1). It was found that osmotic stress limits the dry matter accumulation by approximately 50% during serious water shortages [31,32]. Thus, considering the higher estimated heritability and vulnerability to drought stress, ear biomass could be an optimised option for improving maize selection under waterscarce conditions [33,34]. Plant height was significantly correlated with yield-related traits under both WW and WS conditions [35]. The average PH decreased less when the plant was exposed to drought, which ensures sufficient 'source' availability and exhibits better drought resistance. PH was positively correlated with EBM, suggesting the meaningful role of EBM in drought tolerance. A weak and negative correlation existed between the ASI and ear biomass, indicating that delayed silk extrusion has likely no correlation with ear development. Therefore, the three traits may not be tightly correlated with each other, or the correlation may be disturbed due to the variation of the fluctuating environmental conditions in the field.
Up to now, only few GWAS have been conducted on maize drought tolerance under complex field conditions [30,35]. Researchers tested a number of maize-nested association mapping populations under two contrasting water regimes for seven drought-related traits, including ASI-, PH-, and yield-related traits. Hundreds of promising QTLs and candidate genes were obtained through GWAS and linkage mapping [35]. In addition, many other candidate genes were detected to be associated with drought-tolerance-correlated yield and agronomic traits [30,[36][37][38]; however, none of them were found in this study, which may be due to the various growth conditions and drought treatments. Initially, we aimed to find candidate genes contributing to drought tolerance, especially those for the drought tolerance index. However, only a few genes were identified, and no overlapping genes were identified across different environments, which might be due to the inherence of drought tolerance controlled by multiple minor effect genes. Fortunately, few candidate genes identified in this study were colocalised with reported QTLs. The uncharacterised gene GRMZM2G173084, associated with ASI-WW in 2017, overlapped in a QTL for both ASI-WW and ASI-WS, which was detected by joint linkage analysis in a CN-NAM population [35]. Zm00001d003939, encoding 11-ß-hydroxysteroid dehydrogenase 2, is the candidate gene for PH-WS-2017 located in a consistent QTL for PH-WW, which was identified in both CN-NAM and US-NAM populations. In addition, several overlapping genes were identified under water treatments among different years, while only Zm00001d013992 was commonly identified in a two-year environment. This result implied that maize drought tolerance is a complex trait, highly affected by environment and treatment. Hu et al. [39] reported that grain yield in WS (qWS-GY7-1) and an ear setting percentage in WW (qWW-ESP7-1), located in Chr. 7: 132.2-135. 6 Mb, and seven other QTLs and one mQTL for drought-related traits also clustered in this region [40,41]. Coincidentally, GRMZM2G173084 encoding an uncharacterised protein was associated with ASI-WW in 2017. Additionally, Hu et al. [39] found that the drought tolerance allele of Chr. 7: 132.2-135.6 Mb improved GY in both WW and WS regions, suggesting that both regions influenced GY performance under water-limited conditions. GRMZM2G173084 colocalised in this interval, which indicated its promising role in drought tolerance. Zm00001d044411 for EBM-WS-18 fell in a QTL hotspot in Chr. 3: 219.8-223.7 Mb for GY, ESP, and ASI [39], and its extended region harbours several QTLs responsible for stay-green, leaf senescence, and chlorophyll content identified under normal growth and water stress conditions in maize [29,40,42,43]. Zm00001d032084 (chaperone protein dnaJ) for ASI-Delay-18 and Zm00001d026286 (40S ribosomal protein S11) for ASI-WS-18 colocalised in QTL clusters in bin 1.07 and 10.05-10.07, respectively, each harbouring at least three maize-flowering-time QTLs under different planting densities [44]. Zm00001d024783 (BHLH transcriptional factor 117) for PH-WW-18 and Zm00001d003939 (11-ß-hydroxysteroid dehydrogenase) for PH-WW-17 were localised in QTL (bin 10.07 and bin 2.05-2.07) for plant height of forage maize [45]. These results indicated that these colocalised genes might play important roles in maize development under different conditions, though further functional validations are needed.
Members of zinc-finger family protein play critical roles in plant growth and developmental processes, including flowering, senescence, and also abiotic stress responses [46,47]. A C 2 H 2 zinc-finger transcription factor determines stomatal closure by modulating H 2 O 2 -homeostasis related genes, for example, peroxidases, glutathione S-transferase, and cytochrome P450s [48]. Hence, drought response in rice is regulated. Zm00001d042997 was identified encoding a HIT-type zinc-finger family protein and was associated with ASI-WW and ASI-WS, indicating a potential conserved abiotic-stress-tolerant role of zinc-finger family protein. It was stated that the F-box domain-containing protein ARABIDILLO-1 is conserved in plants, involved in root architecture development, and functions during rice abiotic stress, mainly through regulating root branching and lateral root development [49,50]. It was reported that ARABIDILLO-1-mediated protein degradation, most likely through modulating the GA3 signalling pathway [49]. However, another study revealed that plants with ARABIDILLO-1 knockout and overexpression responded normally to auxin and abscisic acid [51]. Zm00001d029938, encoding ARABIDILLO-1 in maize, was associated with ASI in this study. However, whether it functions in maize drought tolerance was not clearly determined. For EBM, an SNP (S5_27121944), which is annotated as pyridoxal phosphate-dependent transferase (Zm00001d013992), was consistently associated for consecutive two years under drought regime (Table 3) and was upregulated by drought stress. The lead SNP S5_27121944 (A/G) in Zm00001d013992 separated the association panel into two groups, and only around 6% lines (including BY855 and BY4960, etc) carried the favourable haplotype, exhibiting higher ear biomass, 3 g under WS ( Figure 7A). Pyridoxal phosphate (PLP) is an active form of pyridoxine (vitamin B6), which functions as a coenzyme in several reactions such as decarboxylation, deamination, and transamination. The PLP dependent enzymes mainly perform in amino acid biosynthesis and the metabolism of its derived metabolites. Therefore, it is interesting to speculate that Zm00001d013992 might be involved in amino acid metabolism and promotes ear development in maize under drought stress conditions. Potential causal SNPs of candidate drought-tolerant genes could be used for drought-tolerant maize improvement through both genome selection and genome editing.

Conclusions
The findings of this study provide insights into the genetic basis of drought tolerance at the flowering stage, especially for the female inflorescence's development. The overlapping genes are proposed as candidate genes for drought tolerance in maize. Moreover, those lines carrying favourable alleles could be used for drought-tolerance marker development, which is of benefit for future marker-assisted or genome-wide selection for drought-tolerant maize breeding. Future investigation is needed to explore the candidate gene function using CRISPR/Cas9 mediated genome editing and the underlying molecular mechanism of maize ear and silk development under water deficit condition.
Supplementary Materials: The following information is available online at https://www.mdpi. com/article/10.3390/genes13040564/s1, Table S1: List of inbred lines used for GWAS analysis, Table S2: Annotation of SNPs associated with ASI under two water regimes, Table S3: Annotation of SNPs associated with ear biomass at silking date (EBM) under two water regimes, Table S4: Annotation of SNPs associated with plant height under two water regimes, Table S5: Haplotype analysis of Zm00001d013992 for EBM in 2018 and 2020, Table S6: Haplotype analysis of Zm00001d020506 for EBM in 2017, Table S7: Haplotype analysis of Zm00001d029937 for ASI in 2018, Table S8: Haplotype analysis of Zm00001d029938 for ASI in 2018, Table S9: Haplotype analysis of Zm00001d039319 for ASI in 2018, Table S10: Haplotype analysis for Zm00001d042997 for ASI in 2018.