Genomic Selection for Drought Tolerance Using Genome-Wide SNPs in Maize

Traditional breeding strategies for selecting superior genotypes depending on phenotypic traits have proven to be of limited success, as this direct selection is hindered by low heritability, genetic interactions such as epistasis, environmental-genotype interactions, and polygenic effects. With the advent of new genomic tools, breeders have paved a way for selecting superior breeds. Genomic selection (GS) has emerged as one of the most important approaches for predicting genotype performance. Here, we tested the breeding values of 240 maize subtropical lines phenotyped for drought at different environments using 29,619 cured SNPs. Prediction accuracies of seven genomic selection models (ridge regression, LASSO, elastic net, random forest, reproducing kernel Hilbert space, Bayes A and Bayes B) were tested for their agronomic traits. Though prediction accuracies of Bayes B, Bayes A and RKHS were comparable, Bayes B outperformed the other models by predicting highest Pearson correlation coefficient in all three environments. From Bayes B, a set of the top 1053 significant SNPs with higher marker effects was selected across all datasets to validate the genes and QTLs. Out of these 1053 SNPs, 77 SNPs associated with 10 drought-responsive transcription factors. These transcription factors were associated with different physiological and molecular functions (stomatal closure, root development, hormonal signaling and photosynthesis). Of several models, Bayes B has been shown to have the highest level of prediction accuracy for our data sets. Our experiments also highlighted several SNPs based on their performance and relative importance to drought tolerance. The result of our experiments is important for the selection of superior genotypes and candidate genes for breeding drought-tolerant maize hybrids.


INTRODUCTION
Traditional breeding strategies for selecting improved and resistant varieties of maize depending on the phenotypic trait have proven to be of limited success (Cushman and Bohnert, 2000). This direct selection is hindered by low heritability, and the existence of genetic interaction (e.g., epistasis), environmental-genotype interaction, and polygenic effects. This selection also takes a long period of time. Understanding the genetic basis of the plants' response to these various environments and the advent of new genomic technique and tools has allowed breeders to pave a way for selecting superior maize breeds (Tuberosa and Salvi, 2006). Drought stress has the most detrimental effect on maize, leading to reduced yields in maize production (Nepolean et al., 2014). Several QTLs have been identified for drought tolerance in maize and those QTLs have been used to improve stress tolerance through marker-assisted breeding. However, markerassisted breeding is limited to few major QTLs, thus minor QTLs are not part of the selection process, leading to a loss of genetic gain (Dekkers, 2004). To overcome this limitation, genomic selection (GS) has been proposed as a method to understand the effects of all the alleles across the genome to improve polygenic traits (Meuwissen et al., 2001). This method is advantageous over the traditional marker-assisted selection (MAS), as it addresses the effect of small genes which cannot be captured by the traditional MAS (Hayes et al., 2009).
GS is a form of MAS based on breeding values estimated from a genomic dataset that explores the genetic variances within each individual (Heffner et al., 2009). Current research in the area of genetic improvement explores GS as one of the approaches revolutionizing both animal and plant breeding (Hayes et al., 2009;Lorenzana and Bernardo, 2009). Genetic values of quantitative traits in maize and wheat datasets have been studied for the estimation of their higher predictive ability compared with molecular markers than pedigree information (Crossa et al., 2010).
GS reduced the selection time by almost half per cycle compared to the phenotypic selection for almost all traits in the different sets of maize, Arabidopsis and barley (Lorenzana and Bernardo, 2009). By replacing the phenotypic selection with the genomic estimated breeding value (GEBV), the gain for each unit cycle can be increased (Wong and Bernardo, 2008). GS can be appropriate even in the presence of modest molecular markers and diverse environmental conditions (Crossa et al., 2010). The prediction accuracy of breeding values in genomic selection has been found to be 0.58 for grain yield in maize (Zhao et al., 2012) and is estimated to be a better option than other methods considering the genetic gain each year (Lorenzana and Bernardo, 2009;Zhao et al., 2012).
Parametric (RR-Ridge Regression, LASSO-Least Absolute Shrinkage and Selection Operator, Elasticnet, Bayes A and Bayes B), semi parametric (RKHS-Reproducing Kernel Hilbert Space) and non-parametric (RF-Random Forest) models have been used to predict the genotype value, and machine learning programs (Long et al., 2007) have been proposed to develop prediction models for GS. These methods have been implemented in biparental (Lorenzana and Bernardo, 2009) and multi-parental populations (Heffner et al., 2011a) where the predictive ability using several models was compared among different datasets using Arabidopsis, wheat, maize, and barley. Different genomic selection models have been examined in diverse panels of maize and wheat germplasm (De Los Campos et al., 2009;Crossa et al., 2010). GS contributed appreciable genetic gain for grain yield and stover quality in bi-parental maize population (Massman et al., 2013) and drought stress tolerance in tropical maize germplasm (Beyene et al., 2015). In maize, prediction accuracy of GS among the full-sibs was more accurate than unrelated crosses (Riedelsheimer et al., 2013). Among the GS models, rrBLUP and BSSV were found equally efficient in identifying the Stenocarpella maydis resistant maize inbred lines using DArTseq markers (Pedroso et al., 2016).
Little information exists in comparing the efficiency of genomic models to select the better genotypes for drought tolerance in subtropical maize germplasm. The objectives of the present investigation were to predict the GEBVs of genotypes under drought stress using seven GS models, to compare the prediction accuracies of those GS models, and to validate the top selected SNPs from GS models with the SNPs identified through previous GWAS experiment.

Dataset
A set of 29619 cured SNPs, genotyped across a panel of 240 maize inbred lines from an earlier data set (Nepolean et al., 2013(Nepolean et al., , 2014 was used in this experiment. The curation of dataset was done on the basis of MAF <0.05, heterozygosity > 5%, removal of "no calls" (SNPs not included in any cluster were categorized as "no calls"), monomorphs and unmapped SNPs. Briefly, the total genomic DNA from 240 genotypes was isolated with a Nucleopore DNA Sure Plant Mini Kit (Genetix Biotech Asia). SNP detection was performed using the Infinium HD Assay Ultra (Illumina, San Diego, CA, USA). SNP chips were hybridized with 50 ng × 4 µl DNA per sample. The Maize SNP50 BeadChip was used to scan the 240 samples with 24 samples per Sentrix Array Matrix (SAM). All 240 genotypes were genotyped with an Infinium Maize SNP50 BeadChip (Illumina, San Diego, California, USA) containing 56110 SNPs.
Phenotypic data under drought stress in three different environments (IARI, New Delhi: 28 • N 77 • E; 229m AMSL), ANGRAU, Hyderabad: 17 • N 78 • E; 536m AMSL and RRS, Karimnagar: 18 • N 79 • E; 264m AMSL) during post-rainy seasons of 2010/11 and 2011/12 generated earlier (Nepolean et al., 2014) were used for predicting the GEBV in the current experiment. The "mean data" obtained by pooling datasets of these 3 locations was also included in the present experiment. All the drought experiments was followed the alpha lattice design consists of 16 incomplete blocks, and each block comprised of 15 plots with 3 replicates. Drought trials were phenotyped for the anthesis-tosilking interval (ASI, in days), the grain yield (GY, kilograms per plot), the number of kernels per row (KR), the number of kernel rows (KRN), the ear girth (EG, in centimeters), the ear length (EL, in centimeters), and the 100-kernel weight (HKW, in grams). Mean data across location for each genotype was calculated using the restricted maximum likelihood (ReML) approach and the best linear unbiased predictors (BLUPs) was used for further analysis.

Genomic Selection Models
Parametric models (RR, LASSO, EN, Bayes A, Bayes B), a nonparametric model (RF) and a semi-parametric model (RKHS) were used to estimate the genetic value of each genotype. Common variance was considered for all the markers by the ridge regression model (Meuwissen et al., 2001) and, therefore, for each marker effect it constricts uniformly. Estimation of RR βs was performed by minimizing the L 2 panelized residual sum of squares , where RR shrinks all marker effects toward zero rather than categorizing the markers as either significant or as having no effect (Breiman, 1995;Whittaker et al., 2000).
Another parametric model, LASSO, which estimates of the number of βs, was obtained by minimizing the residual sum of squares, and subjected to the constraint of L 1 -type penalty on regression coefficients . EN is a more generalized model that combines both the RR and LASSO penalties. EN's estimate of the number of βs was obtained by minimizing the residual sum of squares subjected to the constraints of both L 1 -and L 2 -type penalties on regression coefficients.
EN simplifies to RR when α = 1 and to LASSO when α = 0. For any other value of α (0 < α <1), EN is used. The L 1 part of EN performs automatic variable selection while L 2 executes grouped selection and stabilizes the solution paths with respect to random sampling. Predictive accuracies and significant SNPs for all traits were estimated through RR, LASSO, and EN with the use of an R package "glmnet" with penalty parameters optimized via tenfold cross validation (Friedman et al., 2010). The significant SNPs estimated on the basis of variable importance were compared with the previous genome-wide association (GWA) results from the water-stressed maize panels (Nepolean et al., 2014).
The other two parametric models, Bayes A and Bayes B (Meuwissen et al., 2001) do not consider the common variance across the effects of SNPs. The Gibbs sampler for 50,000 repetitions fitting the model was computed by discarding the first 5,000 samples as a burn-in and saving one of each of the ten samples for computing the posterior means for parameters. The Bayes A method assumes conditional distribution of each marker effect (given its variance) to follow a normal distribution. If the π value becomes zero, then the Bayes B model shrinks to Bayes A. We used the "BGLR" R package for the implementation of both Bayes A and Bayes B (De Los Campos et al., 2013).
A kernel function is used by the RKHS method to translate datasets of markers into a square matrix to be used in a linear model. There is a possibility that this method might capture nonadditive genetic effects because of its ability to perform non-linear regression in a higher dimensional space. RKHS prediction was performed using "BGLR" (Bernardo and Yu, 2007). The model can be formulated as follows: where ε can be defined as a vector of random residuals and µ as a vector of fixed effects. The parameters α and ε have independent distributions. The matrix K h depends on a kernel function with the smoothing parameter h, which measures the "genomic distance" between genotypes and can be interpreted as a correlation matrix. The "genomic distance" between genotypes is measured by K h , where h represents the smoothing parameter and can be elucidated as a correlation matrix. Here, the Gaussian kernel was used on the genetic distance. The decay rate of correlation between genotypes is regulated by the h parameter. Genomic prediction using parametric and semi-parametric models RR, LASSO, EN, Bayes A, Bayes B, and RKHS was based on 29,619 SNPs, while prediction using a non-parametric model Random Forest (RF) was performed using 5420 SNPs that were randomly selected from a set of 29,619 SNPs. The RF test uses a random subset of predictors to form a collection of regression trees on the basis of bootstrap samples of observations. This model was implemented using the R package "RandomForest" (Liaw and Wiener, 2002) where the number of trees was adjusted to 1,000, keeping mtry at p/3. Prediction accuracy for all agronomic traits was calculated through "Pearson correlation coefficient" between observed and predicted value in all seven GS models.

Validation
To compute predictive accuracies, a 10-fold cross-validation (CV) scheme was applied and iterated 10 times. In each iteration, 10 disjoint subsets of genotypes were formed randomly where one random subset was used as a validation set and the other nine subsets were used as a training population to estimate the parameters of the model used for prediction of excluded genotypes in the validation set. The Pearson correlation between observed and predicted values was calculated in each round. This procedure was iterated 10 times to obtain 100 crossvalidation runs. The predictive ability was calculated as the Pearson correlation coefficient between observations and crossvalidated GEBVs were thus referred for accuracy.

Mean Performance
Under drought condition, the performance of 7 agronomic traits-ASI, GY, KR, KRN, EG, EL, and HKW were recorded in all three locations. In summary, the mean data from all three locations explained that ASI under drought varied from 2 to 12 days with a mean value of 6 days and standard deviation of 2.53. GY had a range of 0.2-2.2 with an average of 1.7 and a standard deviation 1.92. For KRN, the mean and standard deviation was 31 and 2.57, respectively. Other agronomic traits i.e., EL, EG, KR, and HKW, the range varied from 7.8-17.5, 1.8-4.1, 10.9-18 and 15-32 with an average of 13.4, 3.3, 13, and 26 respectively (Nepolean et al., 2014).

Prediction Accuracy of GS Models
Accuracy of the seven GS models was predicted for all seven traits phenotyped for drought stress at the three locations (Table 1). While comparing the prediction accuracies among these traits and locations, we observed that the highest prediction accuracies of 0.93, 0.91, and 0.92 were identified for ASI, EG, and HKW, respectively, in Karimnagar, whereas for GY, Hyderabad and Karimnagar provided the best results, and Karimnagar and IARI showed better results than Hyderabad for EL. The highest prediction accuracy for KRN was identified in Hyderabad with a value of 0.91. KR was the only trait in which all the 3 locations provided consistent results. Across all traits, the maximum prediction accuracy was found for ASI and the minimum for KR. It was examined that among the 3 locations, Karimnagar provided the best results. Prediction accuracy and standard deviations ranged between 0.28-0.92 and 0.03-0.06 (Hyderabad), 0.53-0.92 and 0.02-0.06 (IARI), and 0.30-0.93 and 0.02-0.06 (Karimnagar), respectively, across all traits and models. RKHS, Bayes A, and Bayes B showed no difference in the prediction accuracy above 0.04, while RR, LASSO and EN showed prediction accuracies above 0.03 for all traits, except for KRN, which showed a difference of more than 0.07 at all locations. Bayes B estimated the highest prediction accuracy for all traits except for EG and KRN. Bayes A and RKHS provided the second best prediction accuracy with a drop of 0.01-0.02 (Hyderabad); 0.01-0.03 (IARI), except for GY; and 0.01-0.04 (Karimnagar), except for ASI, respectively, compared to Bayes B.
A great fall in the prediction accuracies for certain traits (ASI, GY, and KRN) under specific models over different locations was observed. Ridge, lasso and EN models predicted ASI with less accuracy (0.28) in Hyderabad. Similarly, in IARI, less prediction accuracy was found for GY (0.6), whereas Karimnagar, ASI (0.3), GY (0.56), and KRN (0.4) had less accuracy compared to the other traits and locations. For GY, the best results (0.9) for Hyderabad and Karimnagar were predicted by both Bayes A and Bayes B, while for IARI, RKHS predicted the best value (0.89). Overall, the prediction accuracy of the mean location was better than that for individual locations.
In location-wise comparison, it was noticed that though the results of RKHS, Bayes A and Bayes B were quite similar but the highest prediction accuracy was obtained from the Bayes B model in all 3 locations namely Hyderabad, IARI and Karimnagar while in the mean dataset RKHS was slightly better over the Bayes B.

SNPs Identified through Different Models
We had estimated the prediction accuracy for seven GS models. From these models, the Bayes B model provided the maximum accuracy for six of the seven traits across several environments. A set of the top 100 SNPs with the highest marker effect observed in each trait and environment was selected using the Bayes B since it produced highest accuracy in all three locations. From this exercise, a total of 2800 SNPs with the highest marker effect across several datasets (traits + environments) were identified. Out of these SNPs, 1053 SNPs were unique (Supplementary Table S1). These SNPs distributed across the genome, ranging from 52 SNPs in chromosome 2-150 SNPs in chromosome 1 (Figure 1).
Out of the 1053 SNPs, a set of 77 consistent SNPs identified across several traits and locations ( Table 2) were selected as a test for understanding their functional relationships with droughttolerant genes. The maize gene models explained that these 77 SNPs distributed across the genome were mapped 10 droughtresponsive TFs within their 150 Kb region. CAMTA mapped within 41 Kb region followed by bHLH (73 Kb The AP2-ERF TF family was mapped close to the maximum number of SNPs (30) on all chromosomes; whereas chromosome 9 had the most drought-responsive SNPs (6) and chromosome 3 contributed only one drought-responsive SNP. The MYB TF family was mapped to 17 SNPs located on all chromosomes except on chromosomes 1 and 10. The SNP PZE-109076471 on chromosome 9 was mapped 2 Kb from MYB TF. Both WRKY and GRAS TFs mapped 10 unique and seven SNPs in their vicinity, respectively. The BHLH TF family encompassed five SNPs on different chromosomes, including an SNP (PZE-110088632) on chromosome 10, which was located only 732 bp away from the BHLH TF. The NF-YA (3), bZIP (2), NAC (1), CAMTA (1), and NF-YB (1) TFs mapped 3, 2, 1, and 1 SNPs, respectively. In addition, SNPs were also mapped close to more than one TF family. The SNP SYN38859 on chromosome 2 was mapped close to MYB and AP2-ERF at a distance of 2 and 116 Kb away, respectively, and another SNP (PZE-107128846) on chromosome 7 was mapped at a distance of 54 and 108 Kb away from the GRAS and NF-YA TFs, respectively.
The set of 1053 SNPs detected using the Bayes B model were matched with the previously identified 67 significant SNPs from GWAS models of GenABEL and GAPIT (Nepolean et al., 2014). We found 10 SNPs which were commonly identified by GS as well as GWAS models (Table 3). These SNPs were mapped on different chromosomes i.e., chromosome 1 (7 SNPs), chromosome 3 (1 SNP), chromosome 4 (1 SNP), and chromosome 10 (1 SNP). All these SNPs were associated with 13 maize gene models and had drought-related functions. Six SNPs were annotated as transcription factors including MYB, bHLH, NF-YA, and FAR1 while rest of them as chaperone protein dnaj 49-like, duf231 domain containing family protein, tubulin beta-1 chain, glutathione peroxidase, and NADP-malic enzyme.

DISCUSSION
Many studies have implemented GS to test the gains in various genetic enhancement programs (Bernardo and Yu, 2007;Wong and Bernardo, 2008;Mayor and Bernardo, 2009;Shengqiang et al., 2009). A high level of correlation between true breeding values and the GEBV is found to be sufficient for genomic selection based on marker data (Heffner et al., 2009). Different approaches have been used to determine breeding values from GS models-penalized regressions (RR, LASSO, and EN), Bayesian approaches (Bayes A and Bayes B), and non-linear regressions (RKHS and RF) (Hayes et al., 2009;Heslot et al., 2012;Nepolean et al., 2013). Non-linear regression models are studied for higher prediction accuracy over penalized models. In our study, we found higher prediction accuracy for both nonlinear models compared to penalized models, with a maximum difference of 0.25 between non-linear models and penalized models across all seven traits. Among several different GS models, a better accuracy is found for non-linear models since they can capture non-additive genetic effects . However, if the additive genetic effects are solely included, using nonparametric models may not yield the expected level of accuracy.
The least prediction accuracy was observed for regularized linear models in this study. This model can be supported by the presence of epistatic interactions which may lower the performance of linear models (Ogutu et al., 2012). Among the penalized models, we observed better prediction accuracy for EN and RR than for LASSO. These results were in agreement with a previous study where EN outperformed LASSO in terms of consistency of model selection and prediction accuracy (Zou and Hastie, 2005).
Bayes B is a variable selection operator, and identifies a subset of markers with larger effects particularly those controlled    (VanRaden et al., 2009;Daetwyler et al., 2010;Jannink et al., 2010). We also observed that RKHS showed the highest prediction accuracy but was restricted to a few datasets.
The variation in prediction superiority for RKHS has also been observed in previous results (Shengqiang et al., 2009;Crossa et al., 2010). The Bayesian model incorporates additive genetic effects, while RKHS captures complex epistatic interactions (Gianola and Van Kaam, 2008). Therefore, one would expect the Bayesian method to perform well in traits where additive effects play a central role and RKHS to perform well in traits where epitasis is more relevant. This also implied that both additive and nonadditive components play significant role in trait expression in variable magnitude depending upon the genetic architecture of the traits (Crossa et al., 2010). Our results showed the presence of variation in predicting the breeding values in different locations which explained that breeding values are shaped up by the environment. The results were also in coherent with our previous GWAS results (Nepolean et al., 2014) where location-specific SNPs were identified. It is also interesting to note that several SNPs were consistent in across locations as well as across traits in the GWAS study.
In our experiment, the Bayesian models out-performed the RR and LASSO models, and this result may be because Bayesian models utilize marker-specific shrinkage of effects, while RR and LASSO equally penalize entire marker effects (Meuwissen et al., 2001). This effect was evident in this study where Bayesian models out-performed the BLUP model by a difference of 0.25 in prediction accuracies. The latter model considers equal variance in all markers, and does not require preliminary information on the variance of marker effects. However, this information is required in Bayesian approaches to estimate prediction accuracies. In addition, RR also incorporates familial relationships and is hence inferior to Bayesian method (Habier et al., 2007). Since the Bayes B method estimates higher prediction accuracies for six of seven traits, it was selected for further validation of SNPs associated with drought tolerance.

Functional Mechanisms of Selected Top SNPs
The SNPs selected based on their marker effects were found to be associated with 10 droughts responsive TFs. The collective role and the interaction of those SNPs with various stress-related mechanisms at a functional level are discussed below.
Gibberellic acid is another plant hormone that promotes growth and cellular elongation; however, its impact on drought stress is not completely understood. GA-responsive transcription factors such as SCARECROW (belonging to GRAS family) are studied for their involvement in cellular differentiation in the Arabidopsis root meristem (Sabatini et al., 2003;Ma et al., 2010). These plant hormones constitute a signaling network that involves various receptors, phosphatases, kinases, calcium-binding signaling molecules, and transcription factors.

Photosynthesis
Maintenance of photosynthesis is an adaptive trait that contributes to an improvement in grain yield in drought stressed plants. The transgenic maize with an increased NF-YB activity is studied for drought tolerance with a maintained photosynthetic rate and high yield (Nelson et al., 2007). These plants also showed good stomatal conductance and high chlorophyll index in water stressed conditions. In addition, the NF-YB transcription factor mapped close to an SNP on chromosome 5 and was co-localized with a QTL for GY, which was previously identified by Nikolić et al. (2012). Another transcription factor CAMTA is identified as a regulator of the photosynthetic machinery, where the T-DNA insertion line of AtCAMTAs has low photo system II efficiency under drought stress (Pandey et al., 2013).

Root Development
A deeper root system is capable of accessing all soil moisture in a drought-stressed plant system. Root development due to drought tolerance has been observed in SNAC1-overexpressing transgenic cotton plants (Liu et al., 2014), and root enlargement has been observed in root-overexpressing OsNAC10 droughtstressed rice plants (Jeong et al., 2010). However, the inhibition of lateral root development is an adaptive response to drought stress (Xiong et al., 2006), whereas reduced lateral root formation and improved drought tolerance have been found in Arabidopsis MYB96-overexpressing mutant plants (Seo and Park, 2009).

CONCLUSION
Our results showed that Bayes B is superior to the other GS models in predicting the genomic values of the studied genotypes. Using Bayes B, we found 77 SNPs that are significant by their marker effects and are related to drought-responsive TFs. We also identified common SNPs from current GS model and previous GWAS models. These significant SNPs are related to many functions, such as stomatal closure, root development, hormonal signaling and photosynthesis. As these SNPs are drought-related and involved in various molecular functions, they can be further used for the development of drought-tolerant hybrids.

AVAILABILITY OF SUPPORTING DATA
The raw SNP data (Submission No. 10.6070/H4BG2KX8) have been submitted to the website http://www.labarchives.com/. All other supporting data are included as additional files.

AUTHOR CONTRIBUTIONS
TN and HSG conceived and designed the experiments; MS, AK, ARR, MGM, and TN analyzed the data. TN, MS, and AK drafted the manuscript. All authors read and approved the final manuscript.

FUNDING
The experiment was funded by the National Agricultural Innovation Project (NAIP, Component IV), the ICAR Network Project on Transgenics in Crop Plants (Maize Functional Genomics Component), and Computational Biology and Agricultural Bioinformatics (Agril.Edn.14(44)/2014-A&P). The funding agencies had no role in the study design, data collection and analysis, the decision to publish, or the preparation of the manuscript.