Numerous genetic loci identified for drought tolerance in the maize nested association mapping populations

Analysis of phenotypes under the WW and WS conditions

The heritability of the seven drought-related traits and average phenotypic performance based on the BLUP values is shown in Additional file 1: Table S1 for CN-NAM and Additional file 2: Table S2 for US-NAM. Significant differences were observed for the means of all of the traits within the two NAM populations under the WW and WS environments using the F test. The results suggested that drought stress at the flowering stage had different effects on drought-related traits. The average grain yield per plant (GYPP) in the CN-NAM population decreased by 28 % under WS, whereas it decreased by 66 % decrease in the US-NAM population under WS. As expected, the average ASI in the two NAM populations was longer under WS than under WW because the ASI reflects the susceptibility of different genotypes to drought stress. The estimated heritability for the seven traits of the CN-NAM population under WS and WW ranged from 49.3 to 83.5 % and from 61.7 to 89.4 %,respectively. The heritability estimates of the seven traits were all higher under WW than under WS.

The phenotypic correlations among the seven traits are listed in Additional file 3: Table S3 and Additional file 4: Table S4 for the CN-NAM and US-NAM population, respectively. Except for ASI and EL under the two water regimes, significant phenotypic correlations were observed among all the traits, and significant negative correlations were observed between the ASI and the remaining traits under both water regimes.

Joint linkage mapping of drought-related traits

To validate the QTLs identified in the CN-NAM population, we also conducted a joint linkage analysis in the US-NAM population (Additional file 6). A total of 112 and 78 joint QTLs were detected for all of the traits in the US-NAM population under the WW and WS conditions, respectively. Among the 112 and 78 QTLs, 32 and 18 QTLs were identified under the WW and WS conditions in the CN-NAM population, respectively, and each trait under the different water conditions had consensus QTLs between the CN-NAM and US-NAM populations except for KNPR under WS (Additional file 7: Table S5).

GWAS of drought-related traits

The genome-wide association study (GWAS) identified 1075 SNPs under WW and 795 SNPs under WS with a bootstrap posterior probability (BPP) ≥0.05 for all of the traits in the CN-NAM population (Additional file 8), and it identified 700 SNPs under WW and 448 SNPs under WS in the US-NAM population (Additional file 9). In certain cases, clusters of associated SNPs were detected within less than 100 kb of each other. The SNPs within a cluster were all identified based on a strong linkage with the same causative genes. Therefore, 100 kb was selected as a non-overlapping sliding window, and the SNPs with the highest statistical significance were selected as representatives of this window. A total of 778 SNPs under WW and 649 SNPs under WS for the CN-NAM population and 586 SNPs under WW and 359 SNPs under WS for the US-NAM population were identified in this manner (Fig. 1).

To identify the most robust associated SNPs, we further selected a BPP ≥ 0.25 as a significance threshold, which detected a total of 255 and 146 strongly associated SNPs for CN-NAM and US-NAM, respectively (Fig. 1). The GWAS results for both NAM populations were compared with both NAM joint QTL intervals. In the WW treatment, 27 % of the strongly associated SNPs identified in both NAM populations overlapped with both NAM joint QTL intervals, whereas in the WS treatment, 29 % strongly associated SNPs overlapped. A total of 221 and 179 strongly associated SNPs were detected in both NAM populations under the WW and WS conditions, respectively. Among the SNPs detected under the different water regimes, 18 associated SNPs were detected in both water regimes, which suggests that these SNPs were constitutive loci in different water environments. Ultimately, 365 strongly associated SNPs were used to identify the candidate genes associated with drought tolerance.

Identification of candidate genes

The predicted genes close to each of the 365 underlying genes might be identified using a publicly available maize genome database (http://www.maizesequence.org), and 354 candidate genes were identified. Among the 365 SNPs, 185 SNPs were located within the coding region of the candidate genes, whereas the remaining 180 SNPs were located closer to the candidate genes, with a physical distance between the SNPs and the candidate genes ranging from 1 to 172,496 bp based on the B73 reference genome v2 (Additional file 10).

Genomic predictions of drought-related traits

Figure 2a shows the accuracy of the genomic predictions for all of the target traits evaluated under the WW and WS conditions using GBS data in the CN-NAM population. The prediction accuracy differed among all of the predicted traits for both water regimes. The accuracy values under WW ranged from 0.51 to 0.80, with a mean of 0.64, whereas the values under WS ranged from 0.40 to 0.74, with a mean of 0.57. Except for the ASI, the accuracy values of other traits under WS were consistently lower than those under WW. Reflecting differences in the heritability and genetic architecture of the target traits, the GYPP and KNPR had small prediction accuracy under the same water regime compared with other agronomic traits. In addition, we conducted genomic prediction for all traits under the WW and WS conditions in the US-NAM population (Fig. 2b). Moreover, similar trends were observed between the two populations, although lower accuracy was observed for the US-NAM population compared with the CN-NAM population.

Plant materials and field environments

The CN-NAM and US-NAM panels were generated as previously described [20]. A total of 1972 CN-NAM lines were grown and measured in the spring of 2009 and 2010 in Urumqi of Xinjiang Province, where the institute of crop science belonging to the Chinese Academy of Agricultural Sciences has set up experimental field bases. The institute of crop science was approved for field experiments, and the field studies did not involve endangered or protected species. For each year, all of the plant materials were subjected to well-watered (WW) and water-stressed (WS) treatments. Trials of each treatment were planted in single-row plots of 11 plants, and there were two replications. Two irrigation regimes were applied using the furrow irrigation method starting at the seeding period. In the WW regime, irrigation was provided in 15-day intervals. In the WS regime, irrigation was given until 3 weeks prior to the expected anthesis date in each CN-NAM family. This stress condition was maintained until 4 weeks after 50 % of the RILs flowered. Nineteen of the 25 US-NAM families, including 2948 RILs, were grown in Urumqi, Xinjiang Province in the spring of 2013. Based on the previous flowering time data [20], four families (CML228, CML247, CML52, and IL14H) with late maturity were excluded from the drought tolerance identification. The P39 and Hp301 families belonging to the sweet and popcorn types were also excluded in the present study. Single-row plots of 11 plants with one replication were grown for each RIL under the WW and WS environments. Each plot was 3 m in length, and the rows were spaced 0.6 m apart. The RIL families were randomly arranged as previously described [39]. Within each family, incomplete blocks consisting of 40 random RILs, the B73 line and alternate parents of the family were planted in an alpha lattice. In the WS regime, the drip irrigation method was applied in 10-day intervals until 3 and 5 weeks prior to the expected anthesis date and after 50 % anthesis occurred. For the WW condition, the soil moisture was maintained at field capacity. All plant materials used in this study were conserved in our experiment lab and we declared that all plant materials used in this study comply with the “Convention on the Trade in Endangered Species of Wild Fauna and Flora”.

Phenotyping

A total of seven traits were measured for all of the lines under both water regimes. The anthesis-silking interval (ASI) was counted as the difference (in days) between male and female flowering times in each plot. The plant height (PH) was calculated as the average height of five random plants measured from the ground to the tassel tip in each plot. Five representative plants in each plot were harvested. The grain yield per plant (GYPP) was evaluated from an average of five plants. The ear length (EL), kernel number per row (KNPR), and ear weight (EW) were measured for five ears and averaged over the plot. The hundred kernel weight (HKW) was estimated from the average weight of 100 randomly selected seeds in three samples.

Genotypic data

A set of 0.95 million SNPs from the CN-NAM and US-NAM populations was generated using Genotyping-By-Sequencing technology [40]. The missing SNPs were imputed using the FILLIN method in TASSEL v.5.0 [41]. These marker sets are publicly available at http://www.panzea.org. The SNP sites exhibiting more than 20 % missing taxa, <5 % minor allelic frequency (MAF), and >20 % heterozygosity were excluded from the raw genotype datasets. Thus, a total of 333,577 and 404,543 SNPs were identified for the CN-NAM and US-NAM populations, respectively. These SNPs were subsequently used in the GWAS and genomic predictions. For the joint linkage mapping, two composite genetic maps that included 4932 and 5296 markers were constructed for the CN-NAM and US-NAM populations, respectively. Detailed information on the maps and genotypes scores has previously been described [20].

Phenotypic data analysis

For the CN-NAM and US-NAM populations, the best linear unbiased prediction (BLUP) for all of the traits of each line across environments (CN-NAM) and within a single environment (US-NAM) was calculated from a random effects model using PROC MIXED in SAS 9.2. In models across the environments of the CN-NAM population, the environment, family, family*environment and entry (family) were considered random effects. The mixed model was fitted within a single environment of the US-NAM population, and the family, RIL within family, blocks, rows and columns were included in the field design. Correlation coefficients were obtained based on the BLUP using Pearson’s statistic applied using the cor function of R software. The broad-sense heritability (h²) for each trait across environments in the CN-NAM population was calculated on a plot basis using the ANOVA tool in QTL IciMapping Version 3.3 [42].

Joint linkage mapping in CN-NAM and US-NAM

The joint linkage analysis for CN-NAM and US-NAM was conducted in SAS 9.2. The detailed information for joint linkage mapping has previously been described [21]. Briefly, PROC GLMSelect was implemented to select the significant marker effects in a family-nested QTL model. For all of the traits, the P-values for the entry and exit of the model were determined using permutation testing. The phenotypic variation explained by all of the QTLs was counted according to Li [43]. The joint linkage QTL support intervals were counted according to Tian [22].

GWAS in CN-NAM and US-NAM

The GWAS was performed in a single NAM population using the Fixed and random model Circulating Probability Unification (FarmCPU) method in R software [44]. To identify the SNPs with the most robust associations with phenotypes, a subsampling procedure was used in the GWA analysis [45]. For each subsampling, 80 % of the RILs of each NAM population were randomly sampled without replacement. This procedure was repeated 100 times between each trait and polymorphism. In each subsampling, we implemented a Bonferroni-corrected threshold probability of 0.05/N to verify the significance levels, where N is the number of individual trait-SNP combinations tested. The significance levels were used to control the false positives in the GWAS. The bootstrap posterior probability (BPP), which is defined as the proportion of times that SNPs were included in the 100 subsamples, was calculated for each significantly associated SNP. Only SNPs with BPP ≥ 0.05 were listed in the results. According to Valdar et al. [45], a BPP ≥ 0.25 was considered the most robust SNP association.

Validation of the candidate genes using RNA-seq data

Based on the maize B73 reference genome assembly V2, genes co-localizing with or adjacent to the associated SNPs were determined to be candidate genes for drought tolerance. Functional annotations of the candidate genes were conducted using blastp, conserved domain search tools, the Maize Genome Database and a literature-specific inspection for each gene. To validate the candidate genes for drought tolerance revealed in the GWAS, the available RNA-seq data of the inbred line B73 was used to evaluate the expression of candidate genes under the WW and WS environments [46]. The RNA-seq data from pollinated ovaries of drought sensitive inbred B73 under WW and WS environments were downloaded from the NCBI web site (http://www.ncbi.nlm.nih.gov/sra/). Briefly, after obtaining the RNA-seq data, the quality control of the raw reads were dealt with the FASTX toolkit (http://hannonlab.cshl.edu/fastx_toolkit/). High-quality RNA-seq reads were mapped to the maize B73 reference genome (B73 AGPv2; http://www.maizesequence.org) using the programs TopHat v2.0.4 [47]. Differential expression analysis was conducted using the HTSeq-DEseq workflow [48]. A false discovery rate (FDR) <0.05 after Benjamini-Hochberg correction for multiple tests was applied.

Cross-validated genomic prediction

We conducted genomic predictions using a mixed-model solver in the rrBLUP package [49] distributed by R software. The prediction accuracy of all traits was evaluated through cross-validation. To perform the cross-validation, a five-fold cross-validation scheme was applied and repeated 20 times. All of the lines in the NAM panel were randomly divided into five disjointed subsets. One of five subsets was selected as the validation population, and the remaining four subsets were used as the training population to estimate the SNP effects for predicting the lines’ values in the validation. The prediction accuracy was calculated according to Pearson’s correlation between the predicted and observed values averaged over all of the cross-validations and replicates.