A panel of 320 spring wheat genotypes [hereafter called wheat association mapping (WAM) panel], representing a world collection from 28 different countries, was used. This diverse panel was procured from the International Maize and Wheat Improvement Center (CIMMYT), Mexico [for details, see Gahlaut et al. (2019)]. The WAM panel was raised under irrigated (IR) and rainfed (RF) conditions each at two different locations: (i) Meerut [North India; (28°0.97′N 77°0.74′E)] and (ii) Powarkheda (Central India; 22°0.07′N 73°0.98′E) in crop season 2011–12 and 2012–13, respectively; thus providing following four environments: Meerut IR (E1), Meerut RF (E2), Powarkheda IR (E3), and Powarkheda RF (E4). Meerut and Powarkheda fall under mega environment (ME)-1 and (ME)-5 respectively, and the soil types of these two locations are deep clay soil and deep loam soil, respectively. At each location, evaluation of WAM panel under both IR and RF conditions minimized the effect of differences in soil types and other environmental conditions on the phenotypic variations at the two locations.

The panel was grown in a simple lattice design with two replications. Each plot, representing an individual genotype, consisted of three 150-cm-long rows with row to row distance of 25 cm. Five irrigations were given under IR condition, while only one irrigation (for sustaining the crop) was given under RF condition (21 days after sowing) to ensure water stress. The details of the experimental locations, sowing dates, meteorological data, and rainfall data are summarized in Supplementary Table 1.

Phenotypic data were recorded for the following four traits: (i) days to anthesis (DTA): recorded as the number of days from planting until anthers in 70% of the spike in a plot had extruded, (ii) grain filling duration (GFD); difference between the DTA and days to maturity (DTM) which was number of days from planting until 70% of the spikes in each plot had turned yellow), (iii) grain number per ear (GNPE); average number of grains per ear using five representative ears per plot, and (iv) grain weight per ear (GWPE); average weight (g) using five ears per plot.

Statistical analysis for all the four traits was conducted using software SPSS 17.0 (http://www.spss.com) to obtain values of range, mean, standard deviation, coefficient of variation (CV), and analysis of variance (ANOVA). To determine the normal distribution, skewness and kurtosis were also analyzed. Correlations between all pairs of traits were obtained separately for IR and RF conditions.

Phenotypic data were used to calculate the following seven indices to estimate the level of drought tolerance: (i) Drought resistance index (DI) = (Ys×(Ys/Yp)/Y¯s (Lan, 1998); (ii) Geometric mean productivity (GMP) = √(Ys × Yp) (Fernandez, 1992); (iii) Mean productivity index (MPI) = (Ys + Yp)/2 (Rosielle and Hamblin, 1981); (iv) Relative drought index (RDI) = (Ys/Yp)/(Y¯s/Y¯p) (Fischer and Maurer, 1978); (v) Stress tolerance index (STI) = (Ys×Yp)/(Y¯ p)ˆ2 (Fernandez, 1992); (vi) Yield index (YI) = Ys/(Y¯s) (Gavuzzi et al., 1997); (vii) Yield stability index (YSI) = Ys/Yp (Bouslama and Schapaugh, 1984), where Ys is the trait value under drought stress and Yp is the value for the same trait under normal condition.

SNP data was generated using the genotyping-by-sequencing (GBS) approach developed by DArT Pty. Ltd., Yarralumla, Australia. The details of the methodology are described by Sehgal et al. (2015). Population structure and linkage disequilibrium (LD) information were available from our earlier study (Gahlaut et al., 2019). Briefly, 9627 SNPs with missing data (<30%) and minor allele frequency (MAF) >5% were used for genotyping. These SNPs were randomly distributed across all the 21 wheat chromosomes spanning 5943.1 cM with an average of 17 SNPs per 10 cM genetic distance. Model-based cluster analysis was performed for population structure using Software STRUCTURE version 2.2 (Pritchard et al., 2000) assuming number of subpopulations (K) to range from 2 to 20, and burn-in and Markov Chain Monte–Carlo (MCMC) iteration were set to 50,000 and 100,000, respectively. The actual number of subpopulations was determined using the tool “STRUCTURE Harvestor” following delta K (ΔK) method (Evanno et al., 2005). LD was estimated using software TASSEL v. 5.0 (Bradbury et al., 2007); mean LD values were obtained for the whole genome as well as for individual chromosomes.

SNPs with no more than 30% missing data and >5% minor allele frequency were utilized for GWAS. Principle component analysis (PCA) was conducted using TASSEL v5.0, and first three components were incorporated as a covariate in association test model. Fixed and random model Circulating Probability Unification (FarmCPU), developed by Liu et al. (2016) was used for ML-GWAS. This method is believed to be the most efficient and eliminates confounding issues arising due to population structure, kinship, multiple testing, etc. The method utilizes both Fixed Effect Model (FEM) and a Random Effect Model (REM), iteratively. FEM allows identification of associated markers (MTAs) that are described as pseudo-quantitative trait nucleotides (pseudo-QTNs) and are used as covariates, in REM, which allows identification of QTNs; in this study, QTNs are described as MTAs. Bonferroni-correction was built in within FarmCPU, so that a default P-value threshold (0.01) was used to declare significant MTAs. Quantile–quantile (Q–Q) plots generated through FarmCPU were used to examine model fitting (account for population structure). Phenotypic differences between the two alleles of SNPs identified as MTA were tested using “t-test” analysis.

Wherever more than two SNPs were associated with the same trait, joint effects were estimated. This was done through linear regression performed using all desirable SNP alleles for the trait (independent variable) and corresponding trait values of the genotypes that contained more than one desirable SNP alleles (dependent variable). Contrasting genotypes were identified using phenotypic values for all the four traits under IR and RF conditions. For this purpose, average sum of ranks (ASR) of seven indices was calculated for each genotype. In the case of DTA, genotypes with higher ASR were considered superior, while for GFD, GNPE, and GWPE, genotypes with lower ASR were considered desirable.

CGs for the associated SNPs were identified by aligning the associated GBS sequences to wheat genome assembly IWGSC refSeq v1.1 available in the Ensemble database (http://www.ensembl.org/info/docs/tools/vep/index.html). The GO annotation (including molecular function and biological process) of each CG was extracted from the IWGSC website (http://www.wheatgenome.org/). Information about expression of CGs was collected using the online tool Genevestigator (Hruz et al., 2008).

The WAM panel exhibited wide range of phenotypic variation for each of the four traits, under both IR and RF conditions at each of the two locations (Figure 1A; Supplementary Table 2). Skewness and kurtosis were found within the range of normal distribution (i.e., ±2.0) for all the traits under all environmental conditions. The only exception was GFD in E4, where kurtosis was observed with a value of 2.2. It is also apparent that each of the four traits was adversely affected under RF condition at both the locations (Meerut and Powarkheda; Figure 1A; Supplementary Table 2).

ANOVA suggested that genotypes differed for the four traits. The performance of genotypes differed under IR and RF and at the two locations (Meerut and Powarkheda) for each of the traits. Genotype × location interactions were significant for DTA and GNPE, while genotype × environment interactions (IR vs. RF) were significant for DTA and GFD (Table 1). Correlations were largely significant under IR and RF conditions at both the locations, except those between GNPE and GFD under RF at Powarkheda (Figure 1B).

The details about population structure and LD are available in our earlier study (Gahlaut et al., 2019). In brief, three sub-populations involving 157 genotypes were recognized; the remaining 163 genotypes were admixed. The number of genotypes was 57 in sub-population I, 85 in sub-population II and 15 in sub-population III. Genome-wide LD decay was observed at 3 cM with a range of 2–20 cM in different genomic regions.

A total of 30 high confidence MTAs were detected for the four traits under one or more environmental conditions. These MTAs involved 27 SNPs distributed over 15 chromosomes (excluding 1D, 4D, 5A, 5D, 7B, and 7D) (Table 2; Figure 2). Manhattan and Q–Q plots showing appropriate model fitting for ML-GWAS are shown in Supplementary Figure 1. Only three SNPs (SNP_404, SNP_1555, and SNP_8047) had association with more than one trait. For individual trait, a maximum of nine MTAs were available for DTA while a minimum of six MTAs were available for GWPE. For DTA, out of nine MTAs, four SNPs (SNP_647, SNP_5369, SNP_7068, and SNP_8390) were uniquely identified under IR condition, four other SNPs (SNP_2283, SNP_2860, SNP_448, and SNP_5304) under RF condition, and one MTA (SNP_404) was identified under both IR and RF conditions. The effect size of these 9 MTAs ranged from 2.2 to 5.2. The effect size of alternate alleles of SNPs associated with DTA was significant in different environments (Supplementary Figure 2). Seven SNPs were associated with GFD, eight with GNPE, and six with GWPE. Effect sizes for associated SNPs were also estimated for individual traits (Table 2).

Based on the effect size, desirable alleles of associated SNPs were also selected (Table 2). Among the four traits, positive selection appeared desirable for GFD, GNPE, and GWPE, and negative selection for DTA. Higher absolute value of SNP effect size showed higher contribution of SNP on the phenotype. Frequency of desirable alleles ranged from 0.10 to 0.94. Two SNPs exhibited pleiotropic effect. SNP_404 was associated with DTA as well as GFD in E1. Similarly, SNP_1555 was associated with GNPE and GWPE in E3.

Joint effect of desirable alleles of multiple associated SNPs was determined using linear regression. For DTA, nine SNPs were each associated with the trait in one or more environments. An increase in the number of desirable SNP alleles from four to nine (but not less than four) led to a significant decrease in DTA (Figure 3). Interestingly, significant joint effect of nine SNPs on DTA was observed in all the four environments; however, strength of regression varied across environments and ranged from 0.33 (E4) to 0.50 (E1).

Significant joint effects of associated SNPs were also observed for GFD, GNPE, and GWPE. Trait values for each of these three traits increased with an increase in the number of desirable alleles in all the four environments (Figure 3). The regression coefficients ranged from 0.08–0.21 for GFD, 0.08–0.15 for GNPE, and 0.07–0.36 for GWPE.

Seven different stress-related indices were obtained for each of the four traits to better assess the genetics of drought tolerance at both locations (total of 28 index traits); this allowed the identification of 153 MTAs involving 85 SNPs (Supplementary Table 3). Manhattan and Q–Q plots showing appropriate model fitting for ML-GWAS tests are shown in Supplementary Figure 3. These SNPs were distributed over 18 of the 21 wheat chromosomes (except 4D, 5D, and 7D). A comparison of MTAs for 28 indices, including seven indices for each of the four main traits allowed the identification of 19 common SNPs associated with response to water stress (Table 2; Supplementary Table 3). As many as 19 SNPs for DTAs, 13 SNPs for GFD, 34 SNPs for GNPE and 28 SNPs for GWPE showed significant association with one or more indices (Supplementary Table 3).

Using ASR values of indices, two contrasting genotypes were selected, which included the superior genotype TX181, and the inferior genotype TX67 (Table 3). TX181 performed better under RF condition for all the four traits (Table 3). These two genotypes can be used for further studies involving crosses generating segregating populations for fine-mapping of QTLs leading to cloning.

The pattern of decline in trait values in RF condition was also examined to assess the sensitivity of these two genotypes to water stress. For this purpose, per cent decline in trait value under RF was examined. Interestingly, in case of TX181, the reduction in trait value under RF condition was relatively low [GFD, 5.11%; GNPE, 20.39%; and GWPE, 17.28%], when compared with those for TX67 [GFD, 10.62%; GNPE, 47.45%, GWPE, 56.45%]. In case of DTA, where lower value is desirable, larger decline was observed in TX181 (4.16%) relative to TX67 (1.24%) under RF condition. However, the DTA in TX67 was about one month longer than TX181. These observations revealed that TX181 is less sensitive (more tolerant) and TX67 is more sensitive under water stress conditions. We may therefore conclude that these varieties differ widely not only for absolute trait values, but more importantly, for tolerance to drought stress.

The above contrasting genotypes were also examined for the presence of desirable alleles of significant SNPs, assuming that desirable alleles for all SNPs are unlikely to be concentrated in TX181; similarly, undesirable alleles for all SNPs cannot be present in TX67. For some of the SNPs, the two genotypes may not differ. Interestingly, out of 27 significantly associated SNPs for four traits, desirable alleles of 24 SNPs were present in TX181. For one SNP, an undesirable allele was found and for two SNPs genotypic data were missing in TX181. In contrast, in TX67, undesirable alleles were present for 14 SNPs, desirables for 9 SNPs, and for the remaining four SNPs, genotypic data were missing (Table 3).

The 27 SNPs involved in 30 MTAs (as mentioned above) were used to identify candidate genes (CGs). As many as 10 of the 27 SNPs were co-located within protein-coding genes and were therefore treated as putative CGs. Details of CGs and their corresponding annotation information are provided in Table 4. Eight of these 10 CGs represent those MTAs that were detected either in IR or RF environments; two MTAs were detected in both the environments. GO annotations of the CGs showed their involvement in protein binding, innate immune response, microtubule cytoskeleton organization, protein self-association, protein phosphorylation and protein folding (Table 4).

Gene expression analysis for the 10 CGs is shown in Supplementary Figures 4–6. This analysis provides further support to their potential involvement in the trait phenotype in different wheat developmental stages and tissues under drought/water stress conditions. The results showed a wide range of expression. For instance, following four CGs had relatively higher expression in all wheat tissues/organ: TraesCS1A02G331000, TraesCS4B02G329500, TraesCS7A02G133300, and TraesCS7A02G176600. Some CGs had tissue-specific expression; for example, TraesCS6D02G394600 expresses in leaf tissues, TraesCS4B02G051200 in root tissues, Traes3B02G596100 in rachis while, TraesCS2D02G574400 in the shoots (Supplementary Figure 4). The CGs also showed varied expression during the different wheat development stages. Interestingly, most CGs had high expression during anthesis to ripening stages except one (Traes3B02G596100), demonstrating their possible role in regulating wheat yield (Supplementary Figure 5).

Under drought/water stress condition, only five of the 10 CGs showed differential expression (≥2-fold), either down-regulated or up-regulated (Supplementary Figure 6). For instance, TraesCS7A02G133300 (up to ~2-fold), TraesCS4B02G051200 (up to ~6-fold), and TraesCS1A02G331000 (up to ~2-fold) were down-regulated. Similarly, TraesCS4B02G329500 (up to ~2-fold) and TraesCS4A02G019800 (up to ~4-fold) were up-regulated during drought stress. Interestingly, out of these five, one CG (TraesCS4B02G051200) that encodes a microtubule-associated protein (MAP65), was associated with MTAs that were identified only in RF environment. Another CG (TraesCS1A02G331000) that encodes RNA helicase protein, belonged to MTAs that were identified in both IR and RF environments (Figure 4).