Genome-Wide Association Study Reveals Genetic Architecture and Candidate Genes for Yield and Related Traits under Terminal Drought, Combined Heat and Drought in Tropical Maize Germplasm

Maize (Zea mays L.) production is constrained by drought and heat stresses. The combination of these two stresses is likely to be more detrimental. To breed for maize cultivars tolerant of these stresses, 162 tropical maize inbred lines were evaluated under combined heat and drought (CHD) and terminal drought (TD) conditions. The mixed linear model was employed for the genome-wide association study using 7834 SNP markers and several phenotypic data including, days to 50% anthesis (AD) and silking (SD), husk cover (HUSKC), and grain yield (GY). In total, 66, 27, and 24 SNPs were associated with the traits evaluated under CHD, TD, and their combined effects, respectively. Of these, four single nucleotide polymorphism (SNP) markers (SNP_161703060 on Chr01, SNP_196800695 on Chr02, SNP_195454836 on Chr05, and SNP_51772182 on Chr07) had pleiotropic effects on both AD and SD under CHD conditions. Four SNPs (SNP_138825271 (Chr03), SNP_244895453 (Chr04), SNP_168561609 (Chr05), and SNP_62970998 (Chr06)) were associated with AD, SD, and HUSKC under TD. Twelve candidate genes containing phytohormone cis-acting regulating elements were implicated in the regulation of plant responses to multiple stress conditions including heat and drought. The SNPs and candidate genes identified in the study will provide invaluable information for breeding climate smart maize varieties under tropical conditions following validation of the SNP markers.


Introduction
The world's population is estimated to reach 10 billion by 2050 [1]. This, together with reduction in arable land and climate change will lead to a rise in abiotic and biotic stresses which are the major threats to food and nutritional security. Maize is an important food security crop and is crucial for economic development in many developing countries in Asia, Latin America, and sub-Saharan Africa [2]. The savannas of sub-Saharan Africa (SSA) contribute to high maize productivity through high incoming solar radiation, reduced incidence of pests and diseases due to low humidity, and low night temperature [3,4]. Additionally, climate change threatens the goal of achieving global food security and could have severe socio-economic consequences globally [5]. With the fast increasing world population, maize production and productivity are expected to be significantly reduced by the adverse impacts of climate change and lead to a global food crisis which is expected to have major bined environments (well-watered, drought and heat stress conditions). Of these SNPs, S1_140960558 on chromosome 1 was associated with days to anthesis under the three conditions while 76 SNPs were linked to grain yield under both drought and heat stresses. These suggested the complex nature of the traits. Drought and heat stresses are known to affect growth and yield of maize significantly as well as predispose the grains to pre-harvest aflatoxin contamination [48,49]. Therefore, understanding the genetic architecture of a trait of interest significantly influences the breeding strategy that will be most effective [50] for genetic enhancement.
Breeding for heat stress as well as drought tolerance in maize is crucial to the effort to overcome the challenges associated with these two key threats to sustainable maize production. Therefore, a total of 162 tropical maize inbred lines developed by the IITA were evaluated under two contrasting environments: heat combined with drought (CHD), and terminal drought (TD) conditions. Furthermore, the inbred lines were genotyped using the genotyping-by-sequencing (GBS) with the Diversity Array Technology (DArT) sequencing (DArTseq) [27][28][29]. The objectives of the study were to (i) map the genomic regions associated with GY and its related traits under CHD conditions, (ii) identify markers linked to the various traits under TD conditions, and (iii) predict the putative candidate genes underlying the genomic regions of the studied traits under stress conditions. The findings from this study will be useful for marker-assisted breeding for the development of resilient maize varieties. Additionally, the predicted putative candidate genes could form the foundation for future functional validations to unravel the molecular mechanisms involved in the response of maize to CHD and TD.

Genetic Resources, Field Evaluation, and Phenotyping
In total, 162 early maturing inbred lines obtained from IITA, Ibadan, Nigeria were used for this study. The lines were evaluated at the Manga Station of the Savannah Agricultural Research Institute, Ghana (11 • 00.977 N, 000 • 15.912 W) and Kadawa, Nigeria (11 • Figure 1). Detailed descriptions of the panel of inbred lines used, treatments (CHD and TD), and traits measured have been published in our previous study [33]. Briefly, the maize plants were exposed to DH. Data were recorded on the following: days to 50% anthesis (AD), days to 50% silking (SD) as the number of days from planting to when 50% of the plant had extruded pollen and produced silks, respectively; anthesis-silking interval (ASI) was recorded as the difference between SD and AD; stalk lodging (SL) was recorded as percentage plants leaning more than 30 • from vertical; husk cover (HUSKC) was recorded 3 weeks post-flowering using a scale of 1-9 (see Table S1); plant aspect (PLASP) was measured based on the overall assessment of the architecture of the plant in a plot as they appeal to sight using a scale of 1-9 (see Table S1); leaf death (LD) was determined by estimating the percentage of leaves dead from the base of the plant to the flag leaf using a scale of 1-9 (see Table S1) at 70 days after planting (DAP); plant height (PLTH) and ear height (EARH) were measured on 10 randomly selected plant from each plot, as a distance from the bottom of the plant to the first tassel branch and to the height of node bearing the upper ear, respectively; root lodging (RL) was measured as the percentage plant bent below or at the highest ear node; ear rot (EARO) was measured by counting the number of ear showing signs of rot within a plot; ear aspect (EASP) and ears per plant (EPP) were taken at harvest. EASP was measured using the general appearance of the ear without the husk using a scale of 1-9 based on the ear size (see Table S1), uniformity of the size, texture, extent of grain filling, insect, and disease damage. Tassel blasting (TB) was measured as the total number of plants with tassel blast and leaf firing (LF) was measured as the total number of plants with leaf firing at the upper part of the plant, 1-2 weeks after silking and tasseling. Grain yield (GY kg/ha) was measured from the shelled grain weight per plot and adjusted to 15% moisture content as outlined by Longmei et al. [14] and Nelimor et al. [18]. The temperature and rainfall data for the experiments at the two test locations were recorded. the size, texture, extent of grain filling, insect, and disease damage. Tassel blasting (TB) was measured as the total number of plants with tassel blast and leaf firing (LF) was measured as the total number of plants with leaf firing at the upper part of the plant, 1-2 weeks after silking and tasseling. Grain yield (GY kg/ha) was measured from the shelled grain weight per plot and adjusted to 15% moisture content as outlined by Longmei et al. [14] and Nelimor et al. [18]. The temperature and rainfall data for the experiments at the two test locations were recorded.

Statistical Analysis
Descriptive statistics: mean, standard deviation (SD), range (minimum-maximum), and coefficient of variation (CV%) of the phenotypic data were computed using SAS version 9.3 (SAS Institute, 2010, Inc., Cary, NC, USA). The frequency distribution and correlation analyses were performed using Performance Analytics R package [51]. Combined analysis of variance was conducted for each of the traits under CHD and TD conditions using a generalized linear model as outlined by Karikari et al. [37]. The broad-sense heritability (H 2 ) was calculated following the formula by [52], = / + + , where is the genotypic variance, is the genotype by environment interaction variance, is the error variance, e is the number of environments, and r is the number of replications.

Genome-Wide Association Study Mapping
The total DArT SNP markers for this panel was 9684 as earlier reported in our study [33]; further filtering was undertaken in pLINK v1.07 with the indep-pairwise 50 10 0.5 command option [53], leaving 7834 SNPs of higher quality.
Association analysis among the SNPs and traits was performed using the mixed linear model (MLM) implemented in Traits Analysis by Association, Evolution and Linkage (TASSEL) version 3.2.3.1 (Figure 1) [54]. The MLM adopted was proposed by Yu et al. [55] with each molecular marker considered a fixed effect and assessed individually:

Statistical Analysis
Descriptive statistics: mean, standard deviation (SD), range (minimum-maximum), and coefficient of variation (CV%) of the phenotypic data were computed using SAS version 9.3 (SAS Institute, 2010, Inc., Cary, NC, USA). The frequency distribution and correlation analyses were performed using Performance Analytics R package [51]. Combined analysis of variance was conducted for each of the traits under CHD and TD conditions using a generalized linear model as outlined by Karikari et al. [37]. The broad-sense heritability (H 2 ) was calculated following the formula by [52], H 2 = σ 2 g / σ 2 g + σ 2 ge n + σ 2 e nr , where σ 2 g is the genotypic variance, σ 2 ge is the genotype by environment interaction variance, σ 2 e is the error variance, e is the number of environments, and r is the number of replications.

Genome-Wide Association Study Mapping
The total DArT SNP markers for this panel was 9684 as earlier reported in our study [33]; further filtering was undertaken in pLINK v1.07 with the indep-pairwise 50 10 0.5 command option [53], leaving 7834 SNPs of higher quality.
Association analysis among the SNPs and traits was performed using the mixed linear model (MLM) implemented in Traits Analysis by Association, Evolution and Linkage (TASSEL) version 3.2.3.1 (Figure 1) [54]. The MLM adopted was proposed by Yu et al. [55] with each molecular marker considered a fixed effect and assessed individually: where Y is the observed vector of means; β is the fixed effect vector (p × 1) other than molecular marker effects and population structure; α is the fixed effect vector of the molecular markers; ν is the fixed effect vector from the population structure; u is the random effect vector from the polygenic background effect; X, W, and z are the incidence matrices from the associated β, α, ν, and u parameters; and ε is the residual effect vector. This model takes into consideration population structure (Q) and relatedness (K), hence, the Q matrix obtained in the population structure analysis was the same as the one reported earlier by Osuman et al. [33] while the K matrix computed in TASSEL was used. Additionally, to re-duce false positive associations, the Bonferroni correction −logP (1/m), where m = number of SNPs, was adopted as the threshold, −logP (1/7834) ≈ 3.89 and this was equivalent to p-value = 1.27 × 10 −4 for declaring the significance of marker-trait association (MTA). Some selected trait association results were visualized using the Manhattan and quantilequantile (QQ) plots in R with package CMplot (https://github.com/YinLiLin/R-CMplot, accessed on 25 October 2020).

Candidate Gene Prediction and In Silico Analyses
Model genes within ±120 kb of stable SNP positions (SNPs detected under either CHD and TD or across (combined) effects ( Figure 1) were downloaded from B73 reference genome (version 4) via the online platform qTeller MaizeMGB (https://qteller.maizegdb.org/, accessed on 14 December 2021) [56]. The expression of model genes was retrieved from qTeller MaizeMGB with control and drought leaf deposited by Forestan et al. [57] and, control and heat treated seedlings by Waters et al. [58], and the expression profiles retrieved were heatmapped via TBtool [59]. The functional annotations such as gene ontology (GO) [60], protein families (Pfam) [61], and Kyoto Encyclopedia of Genes and Genomes (KEGG) [62] of the model genes ±120 kb of stable SNP positions were retrieved from Phytozome version 13 available on https://phytozome-next.jgi.doe.gov/, accessed on 14 December 2021) [63]. Promoter regions of the model genes were analyzed for Cisacting regulatory elements (CAREs) that may be involved in modulating plant response under CHD and TD conditions with 1.5 kb upstream of the start codon (ATG) in PlantCare database [64]. Based on the expression profiles, annotations, and CAREs retrieved, potential candidate genes were predicted. Gene structure was visualized via the online platform Gene Structure Database Server (http://gsds.gao-lab.org/, accessed on 14 December 2021).

Phenotypic Variation, Correlation, and Heritability among the Tropical Maize Germplasm
The rainfall and temperature recorded during the experiments are presented in Figures S1 and S2. The frequency distribution, scatter plot, and correlation among the quantitative traits evaluated/computed are shown in Supplementary Figures S3 and S4 under CHD and TD, respectively. All the agronomic traits measured under the test environments showed continuous and normal distributions which is typical of quantitative traits. For example, the average of the two years phenotypic data from the two countries under CHD experiments had: AD, SD, and GY mean ± standard error (range) values of 66 Table S2a,b) were recorded. These results suggest that CHD and TD stresses alter the phenological, growth, and yield related traits in maize.
Positive and significant correlations (p < 0.05, p < 0.01, p < 0.001) were observed between AD and SD (r = 0.97), HUSKC and PLASP (r = 0.59), and PLASP and EASP (r = 0.54) under CHD ( Figure S3). Similarly, under TD conditions, AD and SD, AD and HUSKC, SD and HUSKC, and EASP and EAROT had a strong correlation (r < 0.60) ( Figure S4). These revealed that these traits as well as others with positive correlations could be selected without any negative effect on the other traits. The majority of the measured traits differed significantly among the lines (genotypes), environment, and genotype by environment interaction (GxE) (Table S3a-c). Average H 2 of the traits was ≈45.54% with the minimum and maximum values of 4.70% (in Manga, Ghana) and 82.00% (in Kadawa, Nigeria), respectively, both under CHD conditions (Table S3a-c). These results indicated that the measured traits may have been influenced by either genetic or environmental factors or both.

Marker-Trait Associations under Terminal Drought Condition
Ten traits viz., AD, EARH, EARO, EASP, EPP, HUSKC, RL, SD, SL, and GY were used for mapping for MTAs employing Q + K in TASSEL. A total of 27 SNPs were linked to these traits with −logP and R 2 ranges of 3.90-24.11% and 12.53-64.79%, respectively,

Mapping Using Data Averaged across Combined Heat and Drought Conditions plus the Terminal Drought Conditions
Eleven common traits selected across CHD and TD conditions, i.e., AD, ASI, EARO, EASP, EPP, LD, PLASP, RL, SD, SL, and GY, with means across these conditions (CHD and TD) were used for MTA mapping. A total of 24 SNPs with a mean number of ≈2.8 SNPs, minimum = 1 SNP (on either Chr07 or Chr08), and maximum = 5 SNPs (on Chr04) were associated with the 11 traits ( Figure 2c; Table 2). These SNPs peaked in the range of 3.94-15.57% and caused 13.05-47.77% phenotypic variation (Table 2). Three SNPs (SNP_203396231, SNP_265278865, and SNP_269120178) on Chr01 were linked to GY, whereas two SNPs (SNP_112491022 and SNP_99059711) were associated with EPP (Figures 3 and 4; Table 2). The hotspot for EASP was on Chr04 with 4 SNPs followed by 3 and 2 SNPs for GY on Chr02 and Chr01, respectively, and 2 SNPs on Chr02 associated with ASI ( Figure 2c; Table 2). One repetitive SNP (SNP_11473743) on Chr01 associated with both AD and SD (Table 2).     (SNP_203396231, SNP_265278865, and SNP_269120178) on Chr01 were linked to GY, whereas two SNPs (SNP_112491022 and SNP_99059711) were associated with EPP (Figures 3 and 4; Table 2). The hotspot for EASP was on Chr04 with 4 SNPs followed by 3 and 2 SNPs for GY on Chr02 and Chr01, respectively, and 2 SNPs on Chr02 associated with ASI ( Figure 2c; Table 2). One repetitive SNP (SNP_11473743) on Chr01 associated with both AD and SD (

Candidate Gene Prediction and In Silico Analysis
With availability of reference genomes and user-friendly bioinformatics tools for insilico prediction, further downstream analyses have become feasible in recent years [65]. Candidate genes are of relevance for functional validation to unearth the regulatory mechanisms underlying each trait of interest. Therefore, we used the 12 common SNPs to retrieve models within the interval position of the SNPs (±120 kb) ( Table 3) via reference genome (B73 v4) through online platform qTeller MaizeMGB (https://qteller.maizegdb.org/, accessed on 15 December 2021) [56]. In all, we identified 55 model genes within 12 SNP positions ±120 kb with varied expression under drought and optimal growing conditions [57] or heat and control conditions [58] (Supplementary Figure S5); out of these, 12 potential candidate genes were predicted (Table 4). With the exception of Zm00001d041124 located at 99054551-99065304 bp on Chr03, the remaining 11 genes have annotations related to either CHD/TD or both conditions (Table 4). For instance, Zm00001d048531 located 117.36 kb (downstream) from SNP_158056460 encodes RNA helicase which has been implicated in improving abiotic stress tolerance in crop plants [66,67]. Most of the predicted candidate genes contain CAREs essential to regulate gene expression under CHD and TD conditions (CGTCA/TGACG (MeJA-responsiveness); ABRE (Abscisic acid); circadian, P-box (Gibberellin responsiveness); AuxRR-core/TGA (auxin); TCA (Salicylic acid); TC-rich repeats (defense and stress responsiveness); (Table S5). Candidate gene structure analyses revealed that Zm0001d002374 possessed 1 exon with no intron while the Zm00001d046434 possessed 17 exons and 16 introns ( Figure 6).  gene expression under CHD and TD conditions (CGTCA/TGACG (MeJA-responsiveness); ABRE (Abscisic acid); circadian, P-box (Gibberellin responsiveness); AuxRR-core/TGA (auxin); TCA (Salicylic acid); TC-rich repeats (defense and stress responsiveness); (Table S5). Candidate gene structure analyses revealed that Zm0001d002374 possessed 1 exon with no intron while the Zm00001d046434 possessed 17 exons and 16 introns (Figure 6). Figure 6. Gene structure analysis of predicted candidate genes. Figure 6. Gene structure analysis of predicted candidate genes.

Discussion
During the last century, conventional breeding approaches have contributed to the development of several high-yielding and superior quality maize breeds [68][69][70][71]. However, the progress achieved so far via conventional approaches is not enough to feed the growing population in the midst of climate change, reduction in arable land, and declining soil fertility. Breeding via conventional approaches is time consuming, expensive, and also phenotypic selection alone does not always result in the development of donor materials with the best potential for trait diversification and improvement. Hence, the present study was conducted to identify SNPs that are linked to the numerous agronomic traits under either CHD or TD conditions or both conditions for marker-assisted breeding.
Under CHD conditions, a total 66 SNPs associated with the 15 measured traits, i.e., AD, ASI, EARO, EASP, EPP, HUSKC, LD, LF, PLASP, PLTH, RL, SL, SD, TB, and YIELD (at Bonferroni-correction −log10(P) ≥ 3.89) were found across the 10 chromosomes in the maize genome (Figure 2a-d; Table S4). Either of these SNPs caused 8.37-49.23% variation in phenotype. This indicated that the 15 traits were complex in nature with numerous loci affecting the phenotypes. The hotspot for YIELD in this study was on Chr04 (SNP_42991074, SNP_49826316, SNP_168609387, SNP_213750918, and SNP_231684692). This is not surprising as Yuan et al. [23] previously reported 4 SNPs for grain yield under heat stress condition. Similarly, chromosomes 2 and 7 were hotpot regions for SL and Chr06 for RL (Figure 2a-d; Table S4). These pinpoint the important roles of these chromosomes in regulating maize response under CHD conditions. A number of earlier studies have reported that several traits such as AD and ASI usually exhibit strong genetic correlations with grain yield [43,72,73] but in this study weak and positive significant correlations (r = 0.25) were observed (Supplementary Figure S3). In this study, four SNPs, SNP_161703060, SNP_196800695, SNP_195454836, and SNP_51772182 on chromosomes 1, 2, 5, and 7, respectively, showed pleiotropic effects on both AD and SD (Table S4). This implied that the four SNPs regulated the two important traits simultaneously and could be targeted for MAS to improve these traits concurrently. This could be the basis for the strong correlation (r = 0.97) coefficient between AD and SD. Flowering in maize (denoted as AD and SD in this study) is an important trait in breeding for drought/heat tolerance. In this study, these two traits were concurrently identified to be linked together, emphasizing their importance in selecting for maize tolerant to drought and combined heat and drought. For example, early flowering days resulted in shorter anthesis-silking interval, increased plant and ear height, increased EPP, and delayed senescence [74]. Additionally, when maize flowers under drought conditions, there is a delay in silking and the interval between anthesis and silking increases thereby giving rise to a longer anthesis-silking interval (ASI) [75].
Terminal drought stress is one of the most vital environmental stress factors that cause a significant reduction in maize productivity. The present study identified markers linked to 10 traits of agronomic significance viz, AD, EARH, EARO, EASP, EPP, HUSKC, RL, SD, SL, and YIELD. Twenty-seven SNPs were associated with these 10 traits at the peak of 3.90-24.11 with 12.53-64.79% phenotypic variation in all chromosomes except chromosome 10 ( Table 1). This implies that breeding by traditional conventional approaches may result in lower genetic gain since these traits are regulated by both major and minor loci [76]. The most prominent SNP (SNP_269120178) on chromosome 1 with −logP (24.11) and R 2 (64.11%) plus two other SNPs (SNP_104170418 and SNP_162101608) on chromosome 2 were linked to YIELD (Table 1. Two SNPs (SNPs: SNP_51109816 and SNP_99059711) were linked to EPP (Table 1). These five SNPs could be targeted for MAS to improve YIELD and ear aspect-related traits under terminal drought condition since grain yield is the key output of production in agriculture. To the best of our knowledge, no GWAS has been reported for grain yield and other secondary traits in terminal drought. Therefore, the information provided in the present study could be useful for breeding for tolerance to terminal drought.
Stability of QTLs/SNPs across multiple environments and genetic backgrounds is a requisite for their utilization in practical plant breeding [77,78]. Hence, the means of the traits (AD, ASI, EARO, EASP, EPP, LD, PLASP, RL, SD, SL, and YIELD) under CHD and TD were used for the mapping of the genes. In all, 24 SNPs associated with these traits across nine chromosomes were used in the present study. Three SNPs on Chr01 (SNP_203396231, SNP_265278865, and SNP_269120178) and 2 SNPs on Chr03 (SNP_112491022 and SNP_99059711) were linked to YIELD and EASP, respectively (Figures 3 and 4; Table 2). These five SNPs could be targeted for allele pyramiding aimed at improving these traits simultaneously. Additionally, the results of the mapping with the means across CHD and TD conditions revealed that some SNPs were detected across either CHD or TD and combined conditions (means) ( Table 3). These SNPs are considered as stable and are recommended for further validation and use in practical maize breeding programs. In all, 105 SNPs were associated with the traits evaluated under the three scenarios. These SNPs could be useful in future genomic selection breeding programs targeted at all the traits at once.
In an effort to select potential candidate genes for cloning and functional verification using gene overexpression and CRISPR/Cas technology, a number of candidate genes were predicted based on propelling evidence from in silico and could be validated to unravel their actual regulatory role in the studied traits. For instance, recent evidence suggest that Auxin induced-like protein (Zm00001d002374) predicted in the present study is reported to play a significant role under drought or heat and combined heat and drought stress in Arabidopsis [79,80]. Circadian clock and CARES are known to regulate gene expression under different conditions [81]. It is of great interest to identify the candidate genes underlying the genomic region for practical plant breeding to unravel the molecular mechanism underlying any trait of interest [38], therefore, the 12 candidate genes predicted together using the CAREs, expression profiles, and gene structure analyses (Tables 4 and S5; Figure 6) provided valuable insight for future functional verification experiments/projects.

Conclusions
Multiple traits from 162 inbred lines and 7834 high quality SNPs were used to conduct MTAs. In all, 117 SNPs associated with the measured traits were evaluated in the present study. Specifically, 66, 27, and 24 SNPs associated with the traits were evaluated under CHD, TD, and combined CHD and TD conditions, respectively. Of these, three SNPs were repetitive under CHD and combined CHD and TD conditions and nine were found concurrently between TD and combined conditions. The highest number of SNPs detected under CHD conditions elucidated the complex nature of the genetic factors that regulated the maize responses under the test conditions. The stable SNPs and those linked to more than one trait could be targeted for MAS. Twelve candidate genes as well as their CAREs, gene structure, and expression profiles were predicted by bioinformatics analyses. The findings from this study offer useful clues about the genetic architecture of these traits under multiple abiotic stresses. The results of this study provided essential information that could be used to improve maize breeding programs in SSA.  Figure S3. Correlation analysis of traits evaluated across the combined heat and drought conditions with average values from Manga (Ghana) and Kadawa (Nigeria). *, ** and *** are significant correlations at p < 0.05, p < 0.01, p < 0.001, respectively. Days to 50% anthesis (AD); days to 50% silking (SD); anthesis-silking interval (ASI); grain yield (GY), plant height (PLTH); root lodging (RL); plant aspect (PLASP); ear aspect (EASP); ear rot (EARO); leaf firing (LF); tassel blasting (TB); ear per plant (EPP); and leaf death (LD). Figure S4. Correlation analysis of traits evaluated under terminal drought in average values from Manga (Ghana) in 2018 and 2019. *, ** and *** are significant correlation at p < 0.05, p < 0.01, p < 0.001, respectively. Days to 50% anthesis (AD); days to 50% silking (SD); grain yield (GY); stalk lodging (SL); ear aspect (EASP); ear rot (EARO); leaf firing (LF); tassel blasting (TB); ear per plant (EPP); and leaf death (LD). Figure S5. Expression of the 55 model genes with the SNP positions ± 120 kb of the markers linked to common traits listed in Table 3. (a). Expression of the genes in leaf tissue under drought and control conditions with T7 maize cultivar used by Forestan et al. [57]. (a). Expression of the genes in leaf tissue under heat and control condition with B73 maize cultivar used by Waters et al. [58]. The genes with dotted red lines are predicted candidate genes based on their annotations retrieved from https://phytozome-next.jgi.doe.gov/ (version 13) (accessed on 14 December 2021). The color legend is shown on right side of the Figure (green and red mean down-and upregulated, respectively, while the black indicates zero expression). Table S1, Description of traits evaluated in this study. Table S2a Table S4. SNPs associated with the 15 traits evaluated under combined heat and drought condition. Table S5. Cis-acting regulating elements identified at 1.50 kb upstream of the start codon of the predicted candidate genes obtained from PlantCARE Database (those with red fonts have been reported to play a role in abiotic stress tol-erance in crop plants).