A bi-parental DH line population, consisting of 169 genotypes and the testcross hybrids, consisting of 190 genotypes was evaluated. The DH population was derived from an F1 cross between drought tolerant lines, La Posta Sequia C7-F64-2-6-2-1-B-B (LPSC7F64), and an elite inbred line from CIMMYT, CML495 (Cairns et al., 2013). The first parental line is traced to the La Posta Sequia Population (LPS), a white dent, Tuxpeño-related synthetic, well adapted to lowland tropics. Full sib recurrent selection was carried out under drought conditions from cycle 0 to cycle 7. The second parental line is a white flint line described as late-lowland with tolerance to rust, helminthosporium, good standability, belonging to heterotic group A. The lines were testcrossed to CML494 for the phenotypic evaluation.

Each of the populations (i.e., hybrids and lines) was evaluated under well-watered (WW) and drought stressed (DS) conditions in experiments conducted in different locations in Mexico during winter cycles of 2013, 2014, and 2016 (Table 1). A total of 11 experiments were conducted in Iguala (Guerrero, Mexico; IG), Agua Fria (Puebla, Mexico; AF) and Tlatizapan (Morelos, Mexico; TL); three for hybrids-WW (IG2013, AF2016, TL2016), four for hybrids-DS (IG2013, TL2013, TL2014, TL2016), two for lines-WW (IG2013, TL2013) and two for lines-DS (IG2013, TL2013).

For all the experiments, the design was an alpha-lattice (0, 1) replicated twice with incomplete blocks size of 5. Plots consisted of one row 4.5 m long at row spacing of 0.75 m. Plots were hand-seeded with two seeds per hill and thinned to one plant per hill (22 plants per plot; 6.6 plants/m²) three weeks after planting.

For the DS treatment, water deficit was induced by withholding irrigation 12-15 days (~190 GDD) before flowering with the aim of reaching the permanent wilting point at flowering at 40 cm soil depth. In the case of severe drought, as indicated by an ASI above 5 d on trial average, irrigation was applied 7 d after completion of anthesis, while a second irrigation was applied 3 wk after completion of anthesis. In the case of moderate DS, only one irrigation was applied during the grain filling period, 2 wk after completion of anthesis. For the WW experiments, evapotranspirated water was fully compensated for through weekly irrigations. Soil moisture content was measured at 10, 20, 30, 40, 60, and 100 cm soil depth three times weekly using Delta-t PR2/6 soil moisture probes (Delta-T devices, Cambridge, United Kingdom) to schedule irrigations in the drought stress treatment. Fertilization, insecticides and herbicides were applied as needed. Fertilizer quantities applied at each location are reported in Table 1.

Several phenotypic traits were measured in each plot throughout the growing season. Two, four, five, and six weeks after planting the NDVI was measured using an RT-505 Greenseeker (Trimble, Ukia, CA, USA). NDVI measurements were taken by running the sensor in the middle of each plot at a height of 80 cm above the canopy. NDVI was calculated per the following equation: R_NIR-R_Red/R_NIR+R_Red, where R_NIR is the reflectance of near infrared wavelength, and R_Red is the reflectance of red wavelength. The NDVI illustrates the part of red wavelength which is absorbed by the plant. At flowering, anthesis and silking dates were recorded when 50% of plants within a plot were shedding pollen and growing silks, respectively. The ASI was calculated as the difference between female and male flowering dates. Two, four and six weeks after flowering senescence was measured visually using a scale ranging from 1 (no senescence) to 9 (complete senescence) to approximate stay green (Trachsel et al., 2016). Leaf rolling was measured visually at flowering, and one and two weeks after flowering using a scale ranging from 1 (unrolled, turgid) to 5 (rolled, onion leaf). For NDVI, leaf rolling and senescence, the area under the curve was calculated by integrating a polynomial function of second degree fitted to individual measurements taken before (for NDVI) or after flowering (for senescence and leaf rolling). In this study, a new secondary trait indicative of early vigor, senescence and overall green leaf area and duration (GLAD) is proposed and evaluated. GLAD was calculated as:

were PHT is plant height at flowering and NDVI and senescence are area under the curve of four NDVI readings and three senescence scores as described above. The area under the curve (AUC) for NDVI and senescence was calculated by integrating a polynomial function of second degree fitted to individual measurements taken before (for NDVI) or after flowering (for senescence).

After physiological maturity was reached, all ears of each plot were collected and shelled, grain moisture was recorded. Grain yield is reported at 12% moisture.

The mixed effect linear model used for the analysis of phenotypic data measured in multilocation trials was:

Where Y_hmlk is the trait value of the h_th genotype (h = 190 and 169 for hybrids and lines, respectively) for the l_th location (hybrids-WW: l = 3; hybrids-DS: l = 4, lines-WW: I = 2; lines-DS: I = 2), the m_th replication (m = 2); μ the overall mean, α_h the main effect of the genotype, E_ml the effect of the location, α_hEml the genotype-by-location interaction, r_m(E_ml) the replication within location effect and r_m(E_ml)δ_k the effect of blocks within replicates within locations and the random error term e_hmlk. All factors except μ were set as random. Best linear unbiased predictors (BLUPs) of genotypes, variance components, and broad sense heritability were obtained. Data for each population-by-treatment combination were analyzed separately. Analysis of the genotype-by-treatment (i.e., WW and DS) interactions was carried out separately for hybrids and lines. For that, terms for treatment and genotype-by-treatment interaction were added to the model described before. Location, genotype and treatment were set as fixed and the rest as random factors. Plant stand was used as a covariate for grain yield and NDVI calculations. Plots with less than 18 plants were removed from the analysis. Variance components were estimated by restricted maximum likelihood (REML) and heritability as the relationship between genetic and phenotypic variance, according to the formula:

where σG2 is the genotypic variance, σGxE2 the genotype-by-environment interaction variance, e the error term, l the number of environments and r the number of replications within environments.

BLUPs for genotypes effects are shrinkage predictors obtained as:

using matrix notation, where y is the vector of the response variable, G ˆ the matrix of variance covariance of the random effects, Z the design matrix for random effects in the model, V ˆ estimated variance of y, 1 a vector of ones and μ the overall mean, the only fixed parameter in the model. The cor.test function in R was used to calculate correlations among BLUPs. Genetic correlations among traits were estimated with a method described previously (Cooper and Delacy, 1994).

For all the maize lines tested in this study, leaf samples bulked from 12 plants of each line were used for DNA extraction with a Cetyltrimethylammonimum bromide (CTAB) procedure (CIMMYT, 2005). A genotyping-by-sequencing (GBS) protocol commonly used by the maize research community was applied at the Cornell University Biotechnology Resource Center in this study (Elshire et al., 2011; Wu et al., 2016). Briefly, the GBS libraries were constructed in 96-plex, and genomic DNA was digested with the restriction enzyme ApeK1. Each library was sequenced on a single lane of Illumina flow cell. SNP calling was performed using TASSEL 5.0 GBS Discovery Pipeline with B73 as the reference genome. Initially, 955,690 SNPs evenly distributed on maize chromosomes were called for each line; 955,120 of them were assigned to chromosomes 1–10, and 570 of them could not be anchored to any of the 10 maize chromosomes. A bin map was constructed by using 20,473 high quality filtered GBS SNPs, details on how to construct the linkage map were described previously (Cao et al., 2017). In brief, neighbor SNPs having high similarity haplotype information were clustered into one bin, and each bin was treated as single marker to construct the genetic map. The following steps were performed to reduce genotyping error and eliminate the low quality SNPs from the bin map: (1) DH lines with heterozygosity rate greater than 5% and/or missing rate greater than 20% were eliminated from further analysis; (2) unlinked SNPs were removed from further analysis, where the window size was 8, similarity rates of all the SNPs within each window were calculated to remove the unlinked SNPs, threshold of similarity rate was 95%; (3) the consecutive SNPs with high similarity rate, i.e., 95%, were merged into one bin; and (4) bins were treated as genetic markers to construct a genetic map. The genetic map was constructed with 191 bins in software QTL IciMapping Version 4.0 (www.isbreeding.net; Wang et al., 2014). The total genetic map length was 987.35 cM resulting in an average distance between markers of 5.15 cM.

Genomic prediction was implemented in rrBLUP package (Endelman, 2011) in DH population. SNPs in the genetic map were used for genomic prediction. Details of the implementation of rrBLUP were described earlier (Zhao et al., 2012). A five-fold cross-validation scheme with 100 replications was used to generate the training and validation sets and assess the prediction accuracy. The average value of the correlations between the phenotype and the genomic estimated breeding values was defined as genomic prediction accuracy (r_MG).

A population of DH lines was evaluated per se and in testcross combination under WW and DS conditions. Significant genotype-by-location interaction was detected when experiments were combined by population and irrigation treatment (i.e., hybrids-WW and hybrids-DS; lines-WW and lines-DS) for all traits, with exception of PHT for lines under DS (Table 2). For most traits, heritability of the combined analysis remained acceptably high. Heritability of PHT and anthesis were the highest in most experiments, with values above 0.75. The ASI had the lowest values, ranging from 0.17 to 0.76. Under DS conditions, grain yield, NDVI and GLAD had the highest heritability values, ranging from 0.45 to 0.81. Meanwhile, senescence and leaf rolling had the lowest heritability, ranging from 0.14 to 0.66.

Grain yield for lines and hybrids was on average reduced by 63% under DS relative to WW (Table 2). Drought stress equally increased ASI (1.94 and 0.38 days for hybrids and lines, respectively), decreased PHT (26 and 10 cm for hybrids and lines, respectively) and delayed anthesis (1.8 and 2.2 days for hybrids and lines, respectively). Differences in NDVI between WW and DS before the onset of drought, can be explained by differences in environments as a result of unbalanced experimental data. Drought trials were all carried out in the winterseason (with lower temperatures compared to the summerseason), while the non-stressed trials were grown in both the summer and winterseason. Since the crop typically develops slower in winter and NDVI readings were taken in calendar days after planting, plants were on average less developed in Winter, relative to plants grown in trials carried out in both summer and winter, resulting in lower NDVI values.

DS were measured at completely dry locations (~900–1,100 masl) whereas additional WW treatments were included in locations with higher precipitations.

Senescence, leaf rolling and GLAD were only recorded under DS conditions. Averaged across treatments, hybrids reached anthesis six days earlier, grew 101 cm taller and had 24% higher NDVI than lines. Leaf rolling, and GLAD were 30 and 177% higher in the hybrids than in the lines and senescence was 11% higher for the lines. The ASI did not differ among lines and hybrids (p > 0.05).

Correlations among grain yield and secondary traits differed across populations and irrigation treatments (Table 3). Grain yield was moderately correlated to NDVI (r_g = 0.84; r_p = 0.53) and GLAD (r_g = 0.70; r_p = 0.49), for hybrids-WW and hybrids-DS, respectively. Grain yield correlated moderately to weakly with PHT, with highest correlation coefficients (r_g = 0.67; r_p = 0.54) for hybrids-WW. Correlations with anthesis, ASI, leaf rolling and senescence were weak or even non-significant for some trait-by-treatment combinations.

The analysis revealed a total of 48 significant QTL for nine traits evaluated (Table 4). They included 13 QTL detected in hybrids-WW, 12 in hybrids-DS, 12 in lines-WW and nine in lines-DS. Thirteen QTL were detected for PHT, eight for grain yield, seven for anthesis, six for senescence, four for ASI, three for GLAD and two each for leaf rolling, NDVI and DSS. In most cases, the allele derived from CML495 increased trait values for anthesis, ASI, PHT, NDVI, DSS, and GLAD, while it reduced trait values for DSS, leaf rolling and senescence. The LOD scores for all detected QTL ranged from 2.0 (grain yield, hybrids-DS) to 7.2 (grain yield, lines-WW) explaining 4.4 (grain yield, lines-WW) to 19.4% (grain yield, lines-WW) of the observed phenotypic variance. Only one constitutive QTL for grain yield was detected, which mapped to bin 8.08 for lines in WW and DS. The trait increasing allele was derived from LPSC7F64 in both cases. None of the 39 QTL detected for secondary traits overlapped for hybrid and line or across treatments.

A collocation of QTL for grain yield (hybrids-WW), anthesis (lines-DS), PHT (hybrids-WW), NDVI (hybrids-WW), and GLAD (hybrids-DS) was detected in bin 1.03 (Table 4). The phenotypic variance explained by each QTL ranged from 4.8 (GLAD, hybrids-DS) to 13.8% (NDVI, hybrids-WW). The allele derived from CML495 delayed anthesis by 1.07 d, increased GLAD by 0.16, grain yield by 0.18 t/ha, NDVI by 0.13 and PHT by 2.89 cm. Another collocation was identified in bin 1.02, where QTL for DSS (hybrids), anthesis (hybrids-DS), PHT (hybrids-WW) and GLAD (hybrids-DS) were detected. The phenotypic variance explained by those QTL ranged from 4.7 (GLAD, hybrids-DS) to 10.9% (anthesis, hybrids-DS). The allele derived from CML495 delayed flowering by 0.26 d, increased GLAD by 0.16, PHT by 1.55 cm and DSS by 1.12%. Another collocation of QTL for PHT (hybrids-WW), leaf rolling (hybrids-DS) and senescence (hybrids-DS) detected in bin 7.04 is indicative of a beneficial effect of early vigor when stress occurs during the post flowering period. The allele derived from CML495 increased PHT and decreased leaf rolling and senescence.

The positive effect of a short ASI on grain yield was confirmed by a collocation of repulsive QTL for ASI and grain yield in bin 4.08 (Table 4). A collocation of repulsive QTL for anthesis and grain yield in bin 4.10 is indicative of the contribution of this chromosomal region to drought escape by early flowering. The grain yield QTL in this bin had the largest phenotypic variance explained among all detected QTL in this study (19.4%), with the trait increasing allele derived from LPSC7F64.

The variation (R²) explained by all QTL for a single trait-by-experiment combination was moderate for grain yield (25.1%, lines-WW), PHT (37.3%, lines-WW) and senescence (25.2%, lines-DS), indicating that genetic control was well captured and is potentially usable in QTL-MAS (Table 5). Lower R² values for the rest of the trait-by-experiment combinations indicate that traits are controlled by many minor effect QTL and genotype-by-environment interactions are high, which are not suitable for QTL-MAS. The prediction accuracy of GS models was larger than the R² values for grain yield and secondary traits, for hybrids and lines both under WW and DS conditions; except for grain yield and ASI in lines-WW and for leaf rolling in lines-DS. The advantages of the GS-MAS over the QTL-MAS approach were larger under DS than under WW conditions for grain yield, anthesis and PHT as indicated in differences between R² (QTL-MAS) and prediction accuracies (GS-MAS); for instance, prediction accuracy of GS and R² values for grain yield were 16.9 vs. 0 and 22.3 vs. 0 for hybrids-DS and lines-DS, respectively. Moreover, the prediction accuracy of grain yield under WW was better than under DS (23.5 and 19.6 for the average of hybrids and lines under WW and DS, respectively) and the prediction of most of the secondary traits were better than for GY except for lines-WW. A similar trend was also observed for the secondary traits. A positive correlation was observed between the genomic prediction accuracy and trait heritability for hybrids-WW (R² = 0.83; p < 0.02, Figure 1). For the other treatment-by-population combinations the correlation between genomic prediction accuracy and trait heritability was low.

The best ten hybrids for WW and DS conditions out-yielded the commercial check (DK357) and the trial mean by 12 and 13% under WW, and by 18 and 11% under DS conditions, respectively (Table 6). Although the genotype-by-water treatment interaction was significant (data not shown) three hybrids with outstanding yield potential and stability (i.e., good performance under WW and DS) were identified (Table 6). The hybrids (LPSC7F64/CML495)DH220/CML494, (LPSC7F64/CML495)DH290/CML494 and (LPSC7F64/ CML495)DH30/CML494, ranked 1st, 2nd, and 3rd under WW conditions, and 22nd, 18^th, and 13th, out of 190 under DS conditions, respectively. On average, the three hybrids with high yield potential and yield stability across environments out yielded the commercial check and the trial mean by 12 and 13% under WW and by 12 and 6% under DS, respectively.

A collocation of QTL for anthesis (hybrids-DS), PHT (hybrids-WW), GLAD (hybrids-DS), and DSS (hybrids) is indicative of the importance of bin 1.02 for the genetic control of grain yield and early vigor. Detection of QTL for anthesis and for DSS in this region indicates that the physiological mechanism conveying grain yield under drought stressed conditions is drought escape achieved through an earlier flowering. The importance of this bin is further supported by collocations with QTL detected for anthesis (Salvi et al., 2011) and PHT (Schön et al., 1993) in other genetic backgrounds. Collocations of QTL were detected for PHT, NDVI and senescence, which are all components of GLAD, in the same bin (Trachsel et al., 2016). Two candidate genes potentially accounting for the detected QTL are ct2 (compact plant2) and cfr1 (coupling factor reduction1). The first gene is involved in CLAVATA signaling, which controls shoot meristem size and shoot growth (Bommert et al., 2013), while cfr1 affects chloroplast function and seedling vigor (Echt et al., 1987).

A collocation of QTL for anthesis (lines-DS), GLAD (hybrids-DS), grain yield (hybrids-WW), PHT (hybrids-WW), and NDVI (hybrids-WW) was detected in bin 1.03, indicative of the importance of this bin for the genetic control of early and general vigor (i.e. NDVI and PHT) and grain yield. Candidate genes for the response related to this chromosome region are a set of genes related to chlorophyll fluorescence (hcf3,hcf31, hcf6) and leaf color (pg^*-N484A, pg^*-N484B, and pg^*-N526C), which may relate to seedling vigor.

A collocation of QTL for PHT (hybrids-WW), leaf rolling (hybrids-DS) and senescence (hybrids-DS) detected in bin 7.04 may suggests that general vigor confers stress avoidance later in the season, since the QTL for PHT was in repulsion with the QTL for leaf rolling and senescence. While PHT reflects general shoot vigor, it may also indicate root vigor (Richner et al., 1996; Hammer et al., 2009; Grieder et al., 2013) as a result of allometric root-shoot relations. Shoot vigor indicative of a vigorous root system, would indirectly allow for greater water and nutrient uptake from deeper soil layers resulting in lower stress levels and reduced leaf rolling and senescence under drought stress. This hypothesis is supported by QTL previously detected for PHT (Sibov et al., 2003) and root architecture in this bin (Tuberosa et al., 2003; Trachsel et al., 2009; Cai et al., 2012). Potential candidate genes underlying the observed response are hcf101, hcf103, and hcf104, which are related to chlorophyll fluorescence conveying sufficient assimilates and plant vigor.

Although a total of 48 QTL were detected for grain yield and secondary traits, none of them was consistently detected in hybrids and lines as a result of the low correlation found among lines and hybrids and across treatments as a result of genotype-by-environment interaction, epistasis and heterosis (Mei et al., 2005; Mihaljevic et al., 2005; Szalma et al., 2007; Hallauer et al., 2010). These results highlight the need to use the testcross' phenotype in mapping studies rather than lines' as done previously (Trachsel et al., 2009, 2010), when aiming to identify QTL to be deployed in hybrids. Only one QTL detected for grain yield in lines was constitutive across treatments. Since there was low correlation of grain yield among lines and hybrids, its usefulness in breeding programs is limited. No QTL for any other trait was detected across treatments, as observed previously (Edmeades et al., 1999). Nevertheless high correlations were found across treatments for anthesis (hybrids and lines) and for PHT (lines).

One constitutive QTL was identified when QTL reported here were compared to results from another study evaluating the same population under nitrogen deficient conditions (DHpop1; Liu personal communication) and an advanced backcross population with a common parent (LPSC7F64; Trachsel et al., 2016) under DS and WW. A senescence QTL under drought (Trachsel et al., 2016) and nitrogen deficient conditions (Liu personal communication) was also detected in bin 4.05; only for lines-DS here, for two populations of hybrids under low N stress (Liu personal communication). These findings are in agreement with two QTL related to senescence detected in this bin by Belícuas et al. (2014) under rain-fed conditions. This QTL has great value for breeding as it could bring yield advantages under two common stresses occurring in tropics (i.e., drought and low soil nitrogen) through improved stay-green. Two candidate genes related to senescence have been reported in this bin. One is SWEET15a, which regulates sucrose translocation in the plant (Chen, 2014). The second is nnr1, which regulates nitrate reductase, a crucial enzyme in nitrite assimilation in plants (Rockel et al., 2002).

In this study, the results showed that the prediction accuracy of grain yield under WW conditions was better than that under DS conditions. The prediction accuracy of the secondary traits were generally higher than the prediction accuracy of GY under almost all the conditions. However, a positive correlation was only observed between the genomic prediction accuracy and trait heritability for hybrids WW. Low correlation between genomic prediction accuracy and trait heritability was observed for all other treatment-by-population combinations. Since the training population was of the same size for all traits the lower prediction accuracy irrespective of the heritability could potentially be attributed to reduced phenotypic variation or large genotype-by-environment interaction (Zhang et al., 2017). Since Combs and Bernardo (2013) additionally show that prediction accuracy may also be dependent on the genetic architecture of a specific trait it is conceivable that genetic and physiological mechanisms acting under stressed conditions are responsible for the observed lack of correlation between trait heritability and prediction accuracy.

Prediction accuracy of GS-MAS was higher than the overall variance explained by all QTL for a trait (R²) in QTL-MAS for grain yield as observed previously (Meuwissen et al., 2001; Bernardo and Yu, 2007; Lorenzana and Bernardo, 2009; Mayor and Bernardo, 2009; Heffner et al., 2010; Guo et al., 2012; Zhao et al., 2012). A similar pattern was observed for secondary traits (i.e., anthesis, PHT, ASI, NDVI, senescence, leaf rolling and GLAD). From a practical point of view, strong QTL remain important in QTL-MAS, as suggested by Heffner et al. (2010). While GS-MAS requires several hundred markers, only flanking markers of target QTL are needed in QTL-MAS. Detected QTL with beneficial effects on early vigor, drought escape, grain yield and stay-green, such as the ones detected in bins 1.02, 1.03, 7.04, and 4.05 could be used in forward breeding to enrich alleles for these traits in a breeding program or for line conversions, while GS-MAS could be used in more mature breeding programs to additionally capture alleles with smaller additive effects (Heffner et al., 2010; Cao et al., 2017). Ideally selection could be carried out for major and minor additive effects simultaneously by using major QTLs as fixed factors in GS-MAS as described by Bernardo (2014).

To be commercially successful, a hybrid needs to perform well under non-stressed and stressed conditions. The fact that no hybrid reached the top ten under both WW and DS conditions is indicative of the difficulty to achieve high grain yield across environmental conditions due to potential physiological tradeoffs between optimal and stressed conditions. However, hybrids (LPSC7F64/CML495)DH220/CML494, (LPSC7F64/CML495)DH290/CML494 and (LPSC7F64/CML495)DH30/CML494 performed well under WW conditions (all of them in the top ten) and drought stressed conditions (all of them were within the best 22 out of 190). Their superior yield potential and stability was reflected by 12% higher grain yield relative to the commercial check (DK357) under both WW and DS conditions. Also, the hybrid (LPSC7F64/CML495)DH109/CML494, ranking 14th and 27th out of 190 under WW and DS conditions, respectively (data not shown) ranked fourth in a study where the same set of hybrids was grown under low nitrogen (Liu et al., personal communication). After further evaluations across sites, in combination with multiple testers, lines DH220, DH290 and DH30 could be released as CIMMYT maize lines for deployment in drought prone environments, while line DH190 could potentially be used in environments prone to drought and low nitrogen.

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.