Inheritance of rust resistance in wheat can be either qualitative or quantitative. Quantitative disease resistance is more durable but more difficult to evaluate because it is expressed in mature plants (adult plant resistance) (Rutkoski et al., 2011). Phenotyping adult plant resistance in large populations is expensive and labor intensive. Expertise is needed because expression is affected by, among other factors, the inoculum load and sequential infection (Hickey et al., 2012). Regarding financial costs, in CIMMYT the most basic field assay costs around US$30 to 40 per genotype (with two replications per location); however, if greenhouse screening is required, the cost increases significantly. In contrast, according to Lorenz et al. (2012), the current cost of genotyping for genomic selection (GS) is $20.
The DNA-based strategies such as marker-assisted selection and GS are an alternative for developing high-yielding wheat with high levels of polygenic adult plant resistance (Heffner et al., 2010; Mago et al., 2011). For example, genotyping can be used to select the best candidates before conducting field tests, thus reducing the number of evaluation plots (Heffner et al., 2010). The complex nature of traits, however, makes marker-assisted selection difficult to implement. For example, gene Sr25 confers a high level of stem rust resistance only when operating in some genetic backgrounds and especially when gene Sr2 is also present (Singh et al., 2011). Besides, the number of molecular markers linked to resistance genes is insufficient to conduct marker-assisted selection (Singh et al., 2011, 2005). Therefore, GS is a promising alternative for accumulating favorable alleles for rust resistance traits (Rutkoski et al., 2011; Xu et al., 2012).
Genetic values for quantitative traits can be predicted as the sum of all marker effects by regressing phenotypic values on all available markers, as first proposed by Meuwissen et al. (2001). There is a plethora of parametric and nonparametric models for genomic prediction. Ridge regression (RR) and the Bayesian least absolute shrinkage and selection operator (LASSO) (BL) (Park and Casella, 2008) are well-established parametric regression methods for GS. Bayesian LASSO has the advantage of being robust regarding prior distributions of the λ parameter controlling shrinkage of regression coefficients (Crossa et al., 2010; de los Campos et al., 2009; Long et al., 2011b). On the other hand, parametric models may not be flexible enough to handle the complexities of gene-to-gene and gene-to-gene-to-phenotype relations (e.g., epistasis, pleiotropy).
Several nonparametric approaches have been proposed for making phenotypic predictions considering these complexities (de los Campos et al., 2010; Long et al., 2011a; Ober et al., 2011). In such models, no assumption is made regarding the genotype–phenotype relationship. Rather, this relationship is described by a smoothing function driven primarily by the data. Support vector regression (SVR) is considered a state-of-the-art, nonparametric, machine-learning algorithm (Hastie et al., 2011; Smola and Schölkopf, 2004) that is designed to learn from finite training set (TRN) data sets; therefore, it seems to be well suited for GS applications. Furthermore, the incorporation of a nonlinear kernel allows modeling nonlinear relationships between phenotypes and marker genotypes (Long et al., 2011b).
This study aims to evaluate GS for predicting resistance to stem rust and yellow rust in five related wheat populations from CIMMYT’s Global Wheat Program. Stem rust was evaluated in two planting seasons, the main and off seasons in Njoro, Kenya, whereas yellow rust was evaluated in Kenya and Mexico (Toluca). Genomic prediction was performed using four different statistical models: BL, RR, and SVR with linear and radial basis function (RBF) kernels.
All five populations were molecularly characterized by a total of 1400 Diversity Arrays Technology (DArT) markers (Triticarte Pty. Ltd., Canberra, Australia). Prediction performance was analyzed using several cross-validation strategies, each assessing four different prediction problems: (i) prediction within populations and environments, (ii) prediction within populations when combining environments, (iii) prediction within populations and between environments, and (iv) prediction between populations.
Materials and Methods
Data sets used in this study come from five wheat populations, PBW343 × Juchi, PBW343 × Pavon76, PBW343 × Muu, PBW343 × Kingbird, and PBW343 × Kenya Nyangumi (K-Nyangumi), evaluated for resistance to both stem and yellow rusts. All populations were derived from crosses between resistant parents (Juchi, Pavon76, Muu, Kingbird, and K-Nyangumi) and PBW343, a moderately susceptible parent. The parents and recombinant inbred lines (RILs) were evaluated for reaction to stem rust at the Kenya Agriculture Research Institute, Njoro, as described by Njau et al. (2010). To analyze resistance to stem race TTKSK (Ug99) using different planting dates (Njau et al., 2010), lines were evaluated in two seasons, the main season (May to October) and the off-season (November to April), in Kenya. The off-season is also known as the wet season whereas the main season is also called the dry season. Percent rust severity for each plot was determined using the modified Cobb Scale (Peterson et al., 1948).
The F6 generation of RILs was also field tested for yellow rust resistance using artificial epidemics at CIMMYT’s research stations in Toluca (Mexico) and natural infection in Njoro (Kenya). The parents and RILs were sown on 75-cm-wide raised beds in paired-row plots, 1 m in length, with 20 cm between rows and a 50-cm pathway. Yellow rust epidemics were initiated about 4 wk after sowing by inoculating spreader rows of susceptible cultivar Morocco planted in hills adjacent to the pathway. To initiate the epidemics, Morocco was sprayed with a suspension of rust urediniospores in Soltrol 170, a lightweight mineral oil (Chevron Phillips Chemical Company, The Woodlands, TX). The yellow rust strains inoculated in Toluca were a pathotype mixture (MEX96.11, MEX07.16, and MEX09.01) whereas in Kenya, yellow rust severity was due to natural infection. The field design for each population was a randomized complete block design with two replicates. The length was of 1 m and the mean of all the plants measured (around 50) were recorded and use for the analyses.
Disease severity of stem rust and yellow rust of the six parents is shown in Table 1. Distribution of the severity for both traits in each population and environment or season is presented in Supplemental Fig. S1; infection levels of the different populations varied in the different environments or seasons.
|Parental line||Stem rust (Kenya)†||Yellow rust (Kenya)||Yellow rust (Toluca)|
The number of wheat lines characterized in each population varied depending on the trait and/or the environment where the populations were evaluated (Table 2). Phenotypic data were processed by square root transformation and standardized to mean zero and standard deviation of 1.
|Population||Off season||Main season||σg2†||σge2†||σe2†||H2†||rg†||Njoro (Kenya)||Toluca (Mexico)||σg2||σge2||σe2||H2||rg|
|PBW343 × Juchi||92||92||96.9||49.7||103.6||0.72||0.47||92||90||100.0||111.1||9.0||0.66||0.52|
|PBW343 × Kingbird||90||90||287.3||45.5||88.2||0.90||0.74||90||90||32.1*||74.8||10.0||0.45||0.35|
|PBW343 × K-Nyangumi||191||176||309.0||69.4||122.6||0.85||0.66||–‡||189||–||–||–||–||–|
|PBW343 × Muu||148||148||247.8||31.6||83.9||0.90||0.75||148||148||21.1*||75.5||16.3||0.66||0.21|
|PBW343 × Pavon76||180||176||181.23||38.5||95.6||0.83||0.64||147||180||76.8||107.9||18.3||0.56||0.48|
Broad-Sense Heritability and Genetic Correlation between Environments
Variance components and broad-sense heritability (H2) are shown in Table 2. Variance components were estimated for each population separately by restricted maximum likelihood (REML) using Proc Mixed of SAS version 9.0 (SAS Institute, 2004). Data from both trials, that is, Main plus Off Season for stem rust and Njoro plus Toluca for yellow rust, were analyzed using the standard linear mixed modelwhere μ is the overall mean, R(E)ijk is the effect of the kth replicate of the jth genotype in the ith environment, Ei is the effect of the ith environment, Gj is the effect of the jth genotype, and GEij is the interaction effect of the jth genotype with the ith environment. All effects were considered random. For both traits the number of replicates per environment was two. Data for yellow rust of population PBW343 × K-Nyangumi was missing in Njoro.
Broad-sense heritability (repeatability) was estimated as H2 = σg2/(σg2 + σge2/e + σe2/re), where σg2 is the genotypic variance, σge2 is the genotype × environment interaction, σe2 is the residual variance, e is the number of environments, and r is the number of replicates per environment.
Genetic correlations (ρg) between environments and seasons (ith and i′th) were calculated using equations from Cooper et al. (1996): ρg = , where is the arithmetic average of all pairwise genotypic covariances between environments i and i′ and is the arithmetic average of all pairwise geometric means among the genotypic variance components of the environments (Cooper et al., 1996).
Genotypic Data and the Genomic Relationship between Parents and Populations
The six parents and the five wheat populations were molecularly characterized using 1400 DArT markers. Molecular data were scored directly into a spreadsheet as present (1) or absent (0). Before each GS evaluation, noninformative markers (i.e., with an allele frequency of less than 5 or more than 95%) were removed for each individual or combined data sets.
The genetic relationship between lines was examined by means of a heat map of the realized genomic relationship matrix G. A common notation is G = t−1XX′ (Van Raden, 2008), where t = , pk is the frequency of allele A, qk is the frequency of allele a in the kth marker, X is the (n × k) matrix of marker data, and n is the number of individuals. The heat map of the G matrix of the six parents is depicted in Supplemental Fig. S2a, and the expected genetic relationship between populations using only molecular information of the parents is shown in Supplemental Fig. S2b. Supplemental Fig. S3 depicts the heat map of the G matrix for lines within a population as well as between populations; four populations are clearly separated in the heat map, but population PBW343 × Juchi does not seem to have lines that are closely related to one another or to lines from other populations. However, we observed that this relationship is completely consistent with expected genetic relationship between populations using only molecular information of the parents.
Linear Genetic Model
The standard linear genetic model considers that the phenotypic response of the ith individual (yi) is explained by a factor common to all individuals (μ), a genetic factor specific to that individual (gi), and the residual comprising all other nongenetic factors (ɛi), for example, environmental effects and effects described by the designs of the experiment. Then, the linear genetic model for n genotypes (i = 1,…,n) is represented as yi = μ + gi + ɛi. In this standard linear genetic model, the genetic factor gi can be described by using a summation of molecular marker effects or by using pedigree.
Ridge regression, one of the first methods proposed for GS (Meuwissen et al., 2001), is equivalent to best linear unbiased prediction (BLUP) in the context of mixed models (Meuwissen et al., 2001). We used the package rrBLUP (Endelman, 2011), which calculates maximum likelihood (ML or REML) solutions for mixed models of the formwhere XF is a full-rank design matrix for the fixed effects βF, Z is the design matrix for the random effects u ∼ N(0,Kσ2), K is a positive semidefinite matrix, and the residuals are normal with constant variance. The function kinship.BLUP performs genomic prediction based on the kinship between lines; the default method of this function is ridge regression, where K = XX′ is the realized genomic relationship matrix computed with the molecular data (Endelman, 2011).
Bayesian Least Absolute Shrinkage and Selection Operator
The model, as implemented in the BLR Package (Pérez et al., 2010), is described aswhere y is an n × 1 response vector, μ is an intercept, and XF, XR, XL, and Z are incidence matrices used to accommodate different types of effects. In our case, XF is a matrix of fixed effects indicating populations or environments, and XR is a matrix of marker genotypes (coded as 0, 1 for markers). The incidence matrix Z is used to accommodate random effects. Finally, ɛ is a vector of model residuals assumed to be distributed as N(0, Diag(σɛ2/ωi2)), σɛ2 is an (unknown) variance parameter, and ωi2 are (known) weights that allow for heterogeneous-residual variances. More details on the BLR (Bayesian linear regression) software are given in Pérez et al. (2010).
Support Vector Regression – Linear and Nonlinear Kernels
Conceptually, the support vector machine algorithm is different from BL and RR. It is based on the structural risk minimization (SRM) principle that aims to learn from finite TRN data sets taking into account the complexity of the hypothesis space and the quality of fitting the TRN data (Hastie et al., 2011). We investigated the ɛ-SVR or “ɛ-insensitive” SVR (Smola and Schölkopf, 2004) implemented in Workbench WEKA (Hall et al., 2009).
The general model can be described as y = β0 + h(x)′β, where y is a real-value response variable, for example, rust resistance, and h(x)T is a linear (or nonlinear) transformation of one (or more) predictors (x) additively combined with the vector of weights (β).
For the sake of clarity, we first describe the most basic linear function, which can take the formwhere 〈ω,x〉 denotes the dot product in χ, the space of the input patterns.
By the SRM principle, generalization is optimized by the flatness of the regression function; this is done by minimizingsubject to
One attractive feature of ɛ-SVR is that it assigns zero value to residuals smaller in absolute value than the constant (ɛ) and a linear loss function for larger residuals (ɛ-tube):
Usually nonlinear models are required to adequately represent and model the data. The simplest solution is to apply a nonlinear transformation to the data in a high dimensional feature space where linear regression is performed. A computationally efficient way to do this is to perform an implicit mapping using kernel functions (Hastie et al., 2011; Smola and Schölkopf, 2004).
We evaluated the performance of the linear kernel, which can be homogeneous, that is, k(xi,xj) = 〈xi,xj〉, or nonhomogeneous, k(xi,xj) = 〈xi,xj〉 + 1, and the RBF kernel, which can be expressed as k(xi,xj) = exp(−γ∥xi − xj∥2). Optimization of parameters was performed by a logarithmic grid search of base 2 over an extensive range of values. Each point on the grid was evaluated by internal fivefold cross-validation on the TRN using Pearson’s correlation as a criterion for success. In our preliminary research, we found no significant impact of ɛ tuning on the success of regression; therefore, given the computational cost involved in a grid search, optimization focused only on the C parameter for the linear kernel and C,γ for the RBF kernel.
Designs for Assessing Different Prediction Problems
Fitting Models Using the Full Data
In this study, we designed four different layouts of the TRN and the testing set (TST), each of them assessing four different prediction cases (Fig. 1). In the first prediction problem (Case A), accuracy is assessed within populations in each environment (season). Data were randomly partitioned 50 times using a TRN containing 90% of the lines and a TST having 10% of the lines. The second prediction problem (Case B) consisted of combining in the TRN the performance of one population evaluated in two environments (seasons or sites) and predicting only one environment (one season or site). For example, the matrix used for the TRN for yellow rust of population PBW343 × Pavon76 in Njoro plus Toluca was generated by combining the matrix used for the yellow rust data of the population evaluated in Njoro and the matrix used on the yellow rust data of the population evaluated in Toluca. This allows assessing the prediction ability of the models using the same TST under different TRN situations (Fig. 1).
The third prediction layout (Case C) involved training the algorithms with the complete information from one environment of each population and evaluating their performance using data from the complementary environment (i.e., the TRN for yellow rust of PBW343 × Pavon76 in Njoro and the TST for yellow rust of PBW343 × Pavon76 in Toluca). Finally, the fourth prediction design (Case D) involved training the algorithms with each of the populations using information from both environments, for example, Main plus Off Season or Njoro plus Toluca, and evaluating their performance separately in each of the remaining populations. It should be pointed out that all possible pairwise comparisons (direct and reciprocal) were performed; however, for simplicity, only direct predictions are reported.
Fitting Models Using Only Ninety Individuals
To examine the potential effect of different sample sizes in the TRN on the prediction of the TST, we examined three of the above-mentioned prediction problems (Cases A, B, and C) while fitting only 90 randomly selected individuals from the three populations with the largest sample size, that is, PBW343 × K-Nyangumi, PBW343 × Muu, and PBW343 × Pavon76 with 176, 148, and 180 lines, respectively. The procedure consisted of a simple random sampling scheme in which 90 lines were randomly selected 20 times from each of the three populations. Then, for each of the 20 simple random samples of 90 lines each, 50 random partitions were performed to split the data into a TRN (90% of the lines) and a TST (10% of the lines).
Software and Statistical Analyses
Scripts for SVR evaluation were implemented in Java code using Eclipse platform version 3.7 (Dexter, 2007). Bayesian LASSO predictions were made with the BLR package for R, version 1.2 (R Development Core Team, 2008; Pérez et al., 2010), and hyperparameters were chosen based on the guidelines of the author. Ridge regression was performed with the package rrBLUP, version 3.8, also for R (Endelman, 2011).
The Shapiro-Wilk test (Shapiro and Wilk, 1965) for normality of the correlations from the 50 random partitions was rejected for Cases A and B; the Wilcoxon test for paired samples was then performed. Statistical differences on the correlations for the cases presented in Supplemental Tables S1 and S2 were performed using the Mann–Whitney U test (Mann and Whitney, 1947). R project (R Development Core Team, 2008) and R Commander (Fox, 2005) were used to perform statistical or routine analyses.
Algorithms were evaluated in a Sun Fire X2250 Server, and the nodes were Quad Core Intel Xeon Processors, Model L5420 (2.50 GHz, 1333 MHz, and 50 W), all with 8 GB of RAM memory. Routine tasks and preliminary analyses were performed using a standard desktop computer, with an Intel Core Duo 2.80 GHz processor and 4 GB of RAM memory.
The total number of lines varied per population; PBW343 × Juchi and PBW343 × Kingbird had only half of the sample size as the other three populations (Table 2). The distribution of stem rust and yellow rust resistance of the five populations in environments varied (Supplemental Fig. S1). Incidence of stem rust ranged from 0 to 85% in the main season and from 0 to 90% in the off-season. The yellow rust incidence ranged from 1 to 75% in Toluca and from 0 to 70% in Njoro.
Variance components and H2 for both traits in each population are given in Table 2. Broad-sense heritability was higher for stem rust (ranging from 0.72 to 0.90) than for yellow rust (ranging from 0.45 to 0.88). Genetic correlations between the off and main seasons for stem rust were relatively high, varying across populations from 0.47 to 0.75, whereas genetic correlations between Njoro and Toluca for yellow rust ranged from 0.21 to 0.52.
In the next two sections, prediction ability results for the four different prediction problems are presented for stem rust and yellow rust.
The positive association between the heritability of the trait and prediction ability (correlation between predicted and observed values) of the BL model for the off and main seasons is depicted in Fig. 2 for the five populations. This result indicates the clear relationship between heritability and predictability of stem rust.
Predicting Resistance within Populations and within and between Environments – Fitting Models Using Full Data (Cases A, B, and C)
The ability of BL, RR, and SVR with linear or RBF kernels to predict stem rust incidence in five populations of wheat when evaluated in two seasons of in Kenya (Main plus Off Season) is given in Table 3. The performance of the algorithms was evaluated under three different prediction cases: Case A, Case B, and Case C.
|PBW343 × Juchi||0.38||0.26||0.41||0.41||0.55||0.55||0.52||0.56|
|PBW343 × Kingbird||0.75**||0.66||0.64||0.68||0.69**||0.61||0.59||0.60|
|PBW343 × K-Nyangumi||0.59**||0.52||0.55||0.56||0.39||0.39||0.32||0.37|
|PBW343 × Muu||0.57*||0.50||0.53||0.56||0.62*||0.53||0.57||0.61|
|PBW343 × Pavon76||0.68**||0.55||0.58||0.60||0.60**||0.46||0.50||0.52|
|PBW343 × Juchi||0.49 (28.9)||0.39 (50.0)||0.39 (−0.01)||0.53** (29.3)||0.58** (5.5)||0.45 (−18.2)||0.45 (−13.5)||0.47 (−16.1)|
|PBW343 × Kingbird||0.85** (13.3)||0.78 (18.2)||0.69 (7.8)||0.82 (20.6)||0.80 (15.9)||0.76 (24.6)||0.59 (0.0)||0.81* (35.0)|
|PBW343 × K-Nyangumi||0.74 (25.4)||0.66 (26.9)||0.78** (41.8)||0.65 (16.1)||0.70** (79.5)||0.67 (71.8)||0.62 (93.8)||0.62 (67.6)|
|PBW343 × Muu||0.72 (26.3)||0.72 (44.0)||0.55 (3.8)||0.75** (33.9)||0.78** (25.8)||0.75 (41.5)||0.60 (5.3)||0.72 (18.0)|
|PBW343 × Pavon76||0.78** (8.8)||0.66 (20.0)||0.71 (22.4)||0.72 (20.0)||0.70 (16.7)||0.67 (45.7)||0.64 (28.0)||0.70 (34.6)|
|PBW343 × Juchi||0.47||0.47||0.47||0.5||0.5||0.47||0.48||0.47|
|PBW343 × Kingbird||0.81||0.77||0.77||0.79||0.80||0.77||0.77||0.80|
|PBW343 × K-Nyangumi||0.61||0.64||0.62||0.62||0.62||0.67||0.67||0.61|
|PBW343 × Muu||0.77||0.77||0.77||0.80||0.80||0.77||0.80||0.77|
|PBW343 × Pavon76||0.70||0.68||0.68||0.68||0.70||0.69||0.63||0.68|
Concerning Case A, for resistance measured in the off-season, the Pearson’s correlation (ρ) ranged from 0.26 (PBW343 × Juchi trained with SVR with linear kernel) to 0.75 (PBW343 × Kingbird trained with BL). For the main season, correlations ranged from ρ = 0.32 (PBW343 × K-Nyangumi and SVR with RBF kernel) to ρ = 0.69 (PBW343 × Kingbird with BL). In most populations, in both the main and off seasons, best predicted values were obtained with BL. This superiority was statistically significant at 0.05 or 0.01 (Wilcoxon test for paired samples) when compared with the other models.
For Case B, the lowest correlation when predicting the off-season was given by SVR with linear kernel for population PBW343 × Juchi (ρ = 0.39), and the highest value was obtained with BL when predicting the resistance of population PBW343 × Kingbird (ρ = 0.85). When predicting the resistance of the individuals in the main season, the worst results were obtained with SVR with linear and RBF kernels for population PBW343 × Juchi (ρ = 0.45) whereas the best result was obtained with RR in population PBW343 × Kingbird (ρ = 0.81). Results indicated an important improvement in prediction accuracy in most of the populations for Case B over Case A (Table 3), except for the population derived from PBW343 × Juchi evaluated in the main season. Prediction ability of the different models differed when trained with data from both seasons as compared with training in one season. For example, BL increased its performance just 5.5% when predicting individuals derived from PBW343 × Juchi in the main season but showed 79.5% improvement when predicting individuals of PBW343 × K-Nyangumi, also in the main season. Support vector regression with RBF showed 93.8% improvement when predicting the same population in the same season (PBW343 × K-Nyangumi and the main season) but its performance decreased 13.5% when predicting lines derived from Juchi in the main season.
Finally, when algorithms were trained exclusively with data from one season and evaluated with data from the other season (Case C), success of prediction on main season (i.e., trained with off-season data) ranged from 0.43 (RR and PBW343 × Pavon76) to 0.80 (BL and RR for PBW343 × Kingbird and SVR with RBF kernel for PBW343 × Muu). When the off-season was predicted, correlations ranged from 0.47 (BL, SVR with linear kernel, and SVR with RBF kernel with PBW343 × Juchi) to 0.81 (BL and PBW343 × Kingbird). In general, training in one season and predicting in the other showed very good prediction ability. Overall, a strong consistency among the three training cases was observed. Population PBW343 × Kingbird always had, on average, the best predictions in both seasons and for the three prediction cases whereas population PBW343 × Juchi gave the lowest correlations in all prediction cases, except for Case A in the main season (the worst was PBW343 × Pavon76). These results can be explained by the higher genomic relationships among the lines of population PBW343 × Kingbird as compared with those comprising population PBW343 × Juchi (Supplemental Fig. S3). Regarding the prediction ability of the models, BL presented the best (mean) correlation in five of eight situations, but RR showed also superior prediction ability in several cases.
Predicting Resistance between Populations (Case D)
The four prediction problems consisted of training the models within each of the five populations using full data from both seasons that had been tested in each of the other populations. Results for the between-populations prediction problem are presented in Table 4 (for simplicity we present only one-way results because reciprocal prediction was, in all cases, very similar to the direct prediction presented here). Over all populations, the SVR models gave the best results, with ρ = 0.42 for the linear and ρ = 0.40 for the RBF kernel. Correlation for BL was 0.31 and for RR was 0.32. It seems that nonparametric models with linear and nonlinear kernels can better capture the relationship between some populations, while parametric models (RR and BL) can capture the relationship in others (Table 4). Interestingly, for stem rust a clear trend was observed between the relatedness of populations and the average correlation (Fig. 3) with the four models.
|Testing or training†
|PBW343 × Juchi||PBW343 × Kingbird||PBW343 × K-Nyangumi||PBW343 × Muu||PBW343 × Pavon76|
|Training or testing‡||PBW343 × Juchi||–||0.48||0.14||0.28||0.31||Bayes LASSO|
|PBW343 × Kingbird||0.53||–||0.29||0.25||0.54|
|PBW343 × K-Nyangumi||0.14||0.30||–||0.28||0.28|
|PBW343 × Muu||0.18||0.30||0.33||–||0.29|
|PBW343 × Pavon76||0.37||0.51||0.22||0.33||–|
|Testing or training§
|PBW343 × Juchi||PBW343 × Kingbird||PBW343 × K-Nyangumi||PBW343 × Muu||PBW343 × Pavon76|
|Training or testing||PBW343 × Juchi||–||0.28||0.32||0.48||0.32||SVR-linear|
|PBW343 × Kingbird||0.28||–||0.45||0.53||0.59|
|PBW343 × K-Nyangumi||0.35||0.31||–||0.25||0.39|
|PBW343 × Muu||0.14||0.58||0.26||–||0.55|
|PBW343 × Pavon76||0.35||0.64||0.41||0.67||–|
Predicting Resistance within Populations and within and between Environments – Fitting Models Using Ninety Randomly Selected Lines
To determine whether, besides its heritability, the size of the TRN had any effect on the predictive ability of the algorithms, we performed a second round of analyses similar to those presented in Table 3 but with the number of individuals of each population fitted to 90 lines. Results are presented in Supplemental Table S1, where the numbers in parentheses and italics represent the decrease in correlation as compared to the case where the full data were used (presented in Table 3). Populations PBW343 × Juchi and PBW343 × Kingbird were not included because they had only 92 and 90 individuals, respectively.
Clearly there is a decreasing effect on the prediction ability of the models for populations PBW343 × K-Nyangumi, PBW343 × Muu, and PBW343 × Pavon76 for the three prediction cases and for most of the models. For Case A, prediction ability decreases 0 to 42%, for Case B, correlations decrease 0 to 38%, and for Case C, the decrease ranged from 0 to 40%.
The association between heritability and prediction ability of the BL model for Njoro and Toluca is depicted in Fig. 2 for four populations (PBW343 × K-Nyangumi was not included). Note that the association is not as clear as in the case of stem rust.
Predicting Resistance within Populations and within and between Environments – Fitting Models Using Full Data
The average accuracy of BL, RR and SVR with linear and RBF kernels for predicting yellow rust severity in five populations of wheat in two environments, Njoro and Toluca, is given in Table 5. In general, prediction ability for yellow rust is lower than that for stem rust, reflecting the generally lower estimated heritability of yellow rust resistance compared to the estimated heritability of stem rust resistance (Table 2). For Case A, ρ values ranged from 0.12 (PBW343 × Muu and RR) to 0.42 (PBW343 × Juchi and SVR-linear) in Njoro whereas in Toluca, ρ ranged from 0.15 (PBW343 × K-Nyangumi and RR) to 0.63 (PBW343 × Pavon76 and BL). In Toluca, BL had the best correlations in two populations, PBW343 × K-Nyangumi and PBW343 × Pavon76, and the same prediction as RR in the other three populations. In Njoro, for PBW343 × Kingbird and PBW343 × Muu, BL gave the best prediction whereas for the other populations both RR and BL gave similar predictions.
|PBW343 × Juchi||0.33||0.42*||0.41||0.33||0.37||0.37||0.34||0.37|
|PBW343 × Kingbird||0.23**||0.16||0.14||0.20||0.49||0.42||0.48||0.49|
|PBW343 × K Nyangumi||–||–||–||–||0.30**||0.24||0.23||0.15|
|PBW343 × Muu||0.17||0.19||0.20*||0.12||0.49||0.45||0.46||0.49|
|PBW343 × Pavon76||0.26||0.20||0.33**||0.25||0.63**||0.54||0.58||0.59|
|PBW343 × Juchi||0.56 (69.7)||0.53 (26.2)||0.52 (26.9)||0.56 (69.7)||0.63 (70.3)||0.60 (62.2)||0.60 (44.1)||0.52 (40.5)|
|PBW343 × Kingbird||0.36** (60.9)||0.30 (100)||0.33 (135.7)||0.32 (60.0)||0.50* (2.0)||0.31 (−26.1)||0.47 (−14.6)||0.45 (−8.2)|
|PBW343 × Muu||0.19 (11.8)||0.15 (−21.1)||0.20 (0.0)||0.14 (16.7)||0.29 (−40.8)||0.34 (−24.0)||0.06 (−87.0)||0.26 (−47.0)|
|PBW343 × Pavon76||0.46 (76.9)||0.46 (130.0)||0.46 (39.4)||0.45 (80.0)||0.62** (−1.6)||0.47 (−13.0)||0.56 (−3.45)||0.60 (1.69)|
|PBW343 × Juchi||0.52||0.52||0.52||0.54||0.52||0.52||0.52||0.56|
|PBW343 × Kingbird||−0.09||0.30||0.31||−0.10||−0.16||0.34||0.42||0.10|
|PBW343 × Muu||0.18||0.13||0.13||0.18||0.17||0.15||0.15||0.17|
|PBW343 × Pavon76||0.48||0.47||0.47||0.48||0.47||0.44||0.41||0.45|
Concerning Case B in Toluca, the lowest correlation was obtained with SVR with RBF kernel PBW343 × Muu population (ρ = 0.06) and the highest value was obtained with BL when predicting the resistance of population PBW343 × Juchi (ρ = 0.63). When predicting the resistance of the individuals in Njoro, the worst results were also obtained with the PBW343 × Muu-derived population but with RR (0.14) followed by the two SVR models (0.15) whereas the best result was obtained with BL and RR and population PBW343 × Juchi (ρ = 0.56).
Prediction results when using both sites, Njoro plus Toluca, to predict one of them indicate a substantial increase with respect to Case A, with the algorithms showing a differential response in prediction accuracy between most populations (Table 5). For example, SVR with linear kernel improved 100 (PBW343 × Kingbird in Njoro) and 130% (PBW343 × Pavon76 in Njoro) and decreased 57.8% when predicting lines of PBW343 × Muu in Toluca. On most of the population, this response depended on the algorithms, that is, the performance of some algorithms increased while that of others decreased. However, population PBW343 × Muu showed a decrease in the performance of all the algorithms when predicting lines in Toluca using information from both trials.
For Case C, when algorithms were trained exclusively with information from Toluca to predict Njoro, correlations ranged from −0.09 (or −0.10) (RR or BL and population PBW343 × Kingbird) to 0.54 (RR and PBW343 × Juchi) whereas in the reciprocal, correlations ranged from −0.16 (BL and population PBW343 × Kingbird) to 0.56 (RR and PBW343 × Juchi). In general, training in one environment and predicting in the other gave good prediction accuracies for populations PBW343 × Juchi and PBW343 × Pavon76 but low correlations for PBW343 × Muu. Correlations for population PBW343 × Kingbird were negligible.
Predicting Resistance between Populations
The other prediction problem is concerned with predicting yellow rust resistance in one population with the aim of predicting the resistance of another. Correlations shown in Table 6 were lower than those presented in Table 3 for stem rust. This was expected because stem rust was generated with only one pathotype. Results show that BL and RR gave relatively higher predictions than the SVR models.
|Testing or training
|PBW343 × Juchi||PBW343 × Kingbird||PBW343 × Muu||PBW343 × Pavon76|
|Training or testing||PBW343 × Juchi||–||0.12||0.06||0.09||Bayesian LASSO†|
|PBW343 × Kingbird||0.06||–||0.16||0.14|
|PBW343 × Muu||0.05||0.16||–||0.07|
|PBW343 × Pavon76||0.10||0.14||0.09||–|
|Testing or training
|PBW343 × Juchi||PBW343 × Kingbird||PBW343 × Muu||PBW343 × Pavon76|
|Training or testing||PBW343 × Juchi||–||0.09||0.05||0.07||SVR-linear§|
|PBW343 × Kingbird||0.09||–||0.04||0.21|
|PBW343 × Muu||0.09||0.11||–||0.09|
|PBW343 × Pavon76||0.12||0.18||0.12||–|
Predicting Resistance within Populations and within and between Environments – Fitting Models Using Ninety Randomly Selected Lines
Results are presented in Supplemental Table S2. Regarding the evaluation of individual data sets (Case A), prediction ability using adjusted data sets decreased between 5 and 61%. On average, algorithms decreased their performance by 28%. In the combined data sets (Case B), performance decreased 0 to 39%, averaging −21% in performance variation. Finally, using only information from one environment to predict the other one the prediction ability of the models decreased, on average, 27% on data sets fitted to 90 individuals compared to the case where the full data were used.
Comparing Model Prediction and Different Prediction Problems
We compared the performance of RR, BL and SVR, with linear or RBF kernel, for predicting stem rust (race TTKSK) and yellow rust resistance of five wheat populations. For stem rust resistance, we analyzed the performance of the algorithms on 10 individual data sets (representing 10 population × environment combinations), each including from 90 to 180 individuals. Performance was evaluated under different training–testing scenarios.
The first scenario, consisting of the random partitioning of the populations in each environment, was proposed considering the problems associated with rust phenotyping—for example, if a program has to evaluate 200 lines in one environment but only has enough resources to phenotype 100 individuals, the genetic value of the remaining candidates can be inferred with relatively high prediction ability as long as the lines in the environment where the model is trained are related to the unobserved lines to be predicted. The second scenario was proposed considering the multienvironments trials commonly used at CIMMYT for evaluating grain yield and diseases (Singh et al., 2012). As shown in Burgueño et al. (2012), combining information from correlated environments can improve prediction accuracy.
The third case allows testing prediction accuracy when the TRN and TST data sets are independent. The algorithm is trained with data from one environment and evaluated using information from the other environment, that is, TRN (Population 1 in Environment A) and TST (Population 1 in Environment B). This situation gives an idea of how to perform GS in other environments (site–year combinations). Finally, in the fourth prediction problem, the TRN and TST data sets are also independent, but now allele effects are obtained by training in one population and used to predict the allele effects in another population, that is, TRN (Population 1 in Environment A) and TST (Population2 in Environment A).
Predictions of stem rust in all prediction cases, populations, and models were higher than those for yellow rust. This may be because these populations were originally designed for stem rust quantitative trait loci mapping (Bhavani et al., 2011) and also because only stem rust race TTKSK was used to inoculate and cause stem rust infection in both seasons whereas yellow rust infection was achieved using natural populations in Kenya or a mixture of pathotypes in Mexico. These and other unknown genotype × environment effects resulted in a much lower genetic variance for yellow rust.
The scenarios presented here attempted to reproduce real situations faced by plant breeders. As reported by Heffner et al. (2010), if net merit (i.e., overall performance) exceeds 0.50, GS could greatly outperform conventional marker-assisted selection in terms of gain per unit time and cost. The results of the three prediction problems (Case A, Case B, and Case C) examined in Table 3 for stem rust may be due to the high heritability of the trait, the high genetic correlation between the off and main seasons (inoculation was performed using only race TTKSK), and the high genomic relationship among the lines comprising the populations, except PBW343 × Juchi. Regarding this population, the low genomic relationship observed among the lines can be attributed to the distance between the parents, that is, PBW343 and Juchi. This observation was confirmed by calculating the expected matrix G based on parent data (Supplemental Fig. S2) and also the G matrix of the five populations (Supplemental Fig. S3).
Additive effects due to multiple minor slow-rusting genes reported in the literature (Singh et al., 2008) may also explain the good results obtained with BL and RR, which in most data sets outperformed the two SVR models. This is not the case for more complex traits; for example, Long et al. (2011b) found that SVR with RBF kernel compared favorably to BL when evaluating a data set for wheat grain yield.
Clearly, when environments were correlated (as the off and main seasons), this favored the prediction of one environment while training in the other (Case B; Tables 3 and 5). In relation to this, differential prediction ability of the models is observed when trained with data from both seasons as compared with training in one season. For example, BL increased its performance just 5.5% when predicting individuals derived from PBW343 × Juchi in the main season (Table 3) but showed 79.5% improvement when predicting individuals of PBW343 × K-Nyangumi, also in the main season. These differences can be attributed to the size of the populations, that is, Juchi has 92 individuals and K-Nyangumi, 186. However, when we fitted the individuals of K-Nyangumi with 90 individuals, a 27.8% improvement in prediction ability was obtained (Supplemental Table S1). Support vector regression with RBF showed a 93.8% improvement when predicting the same population in the main season, but SVR with linear kernel decreased 13.5% when predicting lines derived from Juchi in the main season.
Finally, the same algorithm (SVR with RBF kernel) showed a 20.3% increase in its prediction ability when trained with population PBW343 × K-Nyangumi population but fitted to 90 individuals (Supplemental Table S1). As previously pointed out by Jannink et al. (2010), different models capture different aspects of the relationship between marker genotype and resistance, and thus they could complement each other.
Reducing the number of individuals to 90 lines significantly influenced the success of the prediction, indicating the importance of having large TRN, especially for traits with low heritability.
Prediction between Populations, the Size of the Training Population, and the Two Traits
A difficult prediction problem arises when attempting to predict one population using another. However, as shown in Table 4, some populations are excellent predictors of others when the SVR model is used. For example, data from PBW343 × Kingbird were used to predict PBW343 × Muu (0.73) and also PBW343 × Pavon76 (0.70). Figure 2 shows that the cross products of these pairs of populations in the G matrix had average values of 0.97, 0.95, and 0.84, respectively. Figure 3 shows a strong association between the relatedness of the populations and average accuracy.
Results show that the decrease in size of the TRN population produced a decrease in the models’ prediction ability (Supplemental Table S1) compared with prediction using full data (Table 2), indicating the importance of having large TRN. We should point out that the size of all the populations was only adjusted to 90 to 92 individuals to compare the different populations without the effect of TRN size.
We evaluated nine individual data sets for yellow rust resistance, and predictions were, in general, not as high as those for stem rust resistance. This lower prediction ability may be due to the lower heritability of yellow rust resistance as compared to stem rust resistance. Natural infection of yellow rust in Njoro could contribute to the lower heritability and thus a decrease in GS accuracies. With natural infections there is less uniformity of disease pressure throughout the field. Lorenz et al. (2012) obtained similar results (ρ = 0.72) when predicting Fusarium head blight of barley (Hordeum vulgare L.) using GS. In Rutkoski et al. (2012), ρ reached almost 0.55 when predicting incidence, severity, and kernel quality index of Fusarium head blight in wheat using a similar number of individuals (170).
Although the yellow rust races operating in the two international environments (Njoro and Toluca) were different and heritability was low, the results for yellow rust improved when the algorithms were trained using simultaneous information from Njoro plus Toluca. Correlations of 0.63 for BL were reached when predicting population PBW343 × Juchi in Toluca (Table 5). Also for PBW343 × Juchi and PBW343 × Pavon76, prediction in Njoro while TRN in Toluca or vice versa produced good prediction ability for most models, indicating a good genetic correlation. These results, especially those obtained with stem rust, indicate that is possible to use GS to predict the performance of rust resistance in future environments.
When predicting one population based on another, results show that, in general, most populations were well predicted by any other and with any of the models used. As shown in Table 4, some populations (for example, PBW343 × Kingbird and PBW343 × Juchi) were excellent predictors of others when BL or RR models were used.
In summary, the main conclusions drawn from this study may serve to delineate the implementation of genomics in stem rust and yellow rust resistance work in wheat breeding programs. The additive effects of the two traits resulted in the superior prediction ability of the two linear models RR and BL over the SVR with RBF kernel model. Despite the fact that rust resistance is not governed by as many genes as grain yield, correlation values show that is possible and practical to implement a GS strategy in a rust resistance breeding program.
Results indicate that using data from one environment to predict another correlated environment can be useful in the implementation of a practical GS breeding program and that, even under population structure, good prediction accuracy between populations may be achieved when the TRN and the TST are genetically related. This prediction ability increases when the heritability of the trait increases. Finally, for both stem and yellow rust traits, prediction ability decreased when the sample size of the TRN population decreased. Further research is needed to study genomic-enabled prediction of rust diseases in the presence of genotype × environment interaction and determine whether prediction patterns follow trends similar to those observed for grain yield (Burgueño et al., 2012).
Supplemental Information Available
Supplemental material is included with this article.