Genomic selection (GS), as an efficient breeding tool, was first implemented in animal genetic improvements (Schaeffer, 2006; Meuwissen et al., 2013). With the availability of genome-wide markers and low-cost genotyping technologies, GS is rapidly adopted in plant genetics and breeding (Heffner et al., 2009; Lin et al., 2014). In GS, a training population, which has been genotyped and phenotyped, is used to train a statistical model to predict individuals that have been genotyped but not phenotyped. The predicted value is termed as an estimated breeding value (EBV) if the model utilized a pedigree relationship between individuals and a genomic estimated breeding value (GEBV) if marker data is used in the analysis. The EBVs or GEBVs can be used to rank and select germplasm (Meuwissen et al., 2001; de Los Campos et al., 2013).

Implementing GS in the breeding program heavily depends on the prediction accuracy of GEBVs. Multivariate models showed higher prediction accuracy than univariate models in GS studies (Jia and Jannink, 2012; Sun et al., 2019). The additional information in genetically correlated traits is exploited in multivariate models, and the higher the correlation is, the greater the multivariate models would benefit (Zhao et al., 2022a). Rutkoski et al. (2016) included canopy temperature and normalized difference vegetation index in a multivariate model, which resulted in a 70% prediction accuracy improvement for grain yield in wheat. Important agronomy traits, such as days to flowering, plant height, etc., with malting quality traits, were used to assess the prediction accuracy of multivariate models, and the result showed a 76% higher predictive ability than univariate models (Bhatta et al., 2020).

GS can improve the genetic gain of target traits in crop breeding programs mainly by increasing selection accuracy, shortening breeding cycles, etc. (Falconer and Mackay, 1996; Xu et al., 2017). Studies to compare genetic gains achieved by different selection strategies have been conducted with real datasets and simulations (Woolliams et al., 2015; Marulanda et al., 2016). Simulation with ryegrass showed a 4-year reduction in cycle time, with genetic gain doubling or tripling when GS is incorporated into the breeding program (Lin et al., 2016). In the short breeding cycle, improving the selection accuracy would be important for higher genetic gain. Beyene et al. (2015) estimated genetic gains in grain yield in eight bi-parental maize populations and observed higher grain yield achieved from GS than conventional pedigree-based phenotypic selection (PS). A study with wheat indicated that the integrated PS + GS approach could result in an optimal genetic gain for yield (Lozada et al., 2020).

GS can simultaneously select multiple traits with the genomic selection index (GSI), a linear combination of ‘traits’ GEBVs. Both simulation and real data showed that GSI was more efficient than the phenotypic selection index (PSI) per unit of time (Ceron-Rojas et al., 2015; Habyarimana et al., 2020b). A study with maize confirmed that GSI is more efficient in obtaining sustainable gain than the phenotypic selection index (PSI) per unit of time in selection for grain yield, plant height, etc. (Ceron-Rojas et al., 2015). Various GSI were compared in wheat, and the results indicated the employed GSI mitigated the negative trade-off between grain yield and protein content and led to a substantial selection response for protein yield (Michel et al., 2019). In multitrait index selection, the choice of the selection index and the assignment of the weights to different traits strongly affected the selection gain (Marulanda et al., 2021).

Safflower is a multi-purpose crop that is grown worldwide in Africa, America, Europe, and Asia (Zohary, 1999). It is grown as a vegetable, cut flower, herbal medicine, bird feed, etc. However, the main growing interest has shifted to safflower seed oil because of its high oleic and linoleic acid content (Li and Mündel, 1996; Khalid et al., 2017). In 2020, the global safflower growing area was about 0.8 million ha, and it was around 40,000 ha in Australia (FAO, 2020). Safflower has the potential to be further incorporated into the farming system because of its drought tolerance. Two Australian safflower varieties released in 1998 carried resistance to leaf spot (Alternaria Carthami) and root rot (Macrophomina phaseolina), respectively, and most varieties grown today in Australia were developed overseas (Jochinke et al., 2008). Breeding for high-yield safflower varieties is under demand. In this study, we attempted to evaluate the genetic gain by adopting genomic selection to select elite crossing parents from a diverse Genebank collection. The aims of the study were 1) to estimate the genetic correlation between major safflower agronomy traits; 2) to evaluate the genomic prediction accuracy for safflower grain yield by combining other traits with multivariate prediction models; 3) to compare selection responses of grain yield to different safflower selection strategies.

Safflower is an underutilized oil seed crop. Breeding efforts for genetic improvement of target traits were modest. With limited genetic and genomic resources and limited research funds, the datasets we used in our study have been used for genetic diversity study, genetic characterization, and genomic parameters estimation of important agronomic traits in safflower. All those studies enhanced our knowledge of safflower and paved the way to transfer conventional safflower breeding into a highly targeted and more efficient modern crop breeding scheme. In this study, we further investigated the genomic prediction accuracy for grain yield with multivariate models and compared grain yield selection responses to different selection strategies. Combining GS into the conventional breeding program and optimizing the breeding strategies would facilitate and fast-track the genetic improvement for safflower.

Genetic parameters are important estimates for quantitative traits. In this study, we nearly doubled the SNPs number compared to our previous study (Zhao et al., 2021). A high marker density could affect the heritability estimation but not impact the prediction accuracy too much (Zhao et al., 2021). To increase the marker density, imputation of SNPs with missing rate up to 80% has been reported in the GS study for ryegrass and soybean (Jarquín et al., 2014; Faville et al., 2018). The estimated variance and heritability via the multivariate model in our study showed consistent trends with the univariate GBLUP models (Zhao et al., 2021). With the decomposition of covariance between traits in multivariate model, the h² could be estimated more accurately compared with the univariate model. Therefore, our study indicates that the multivariate models may perform as well as the univariate models in genetic parameters estimation.

Combining the correlated trait with high heritability would improve the prediction accuracy for the low heritability traits by the multivariate model. The correlation between traits could be a key factor determining the multivariate model’s advantage over the univariate model (Montesinos-López et al., 2021). In our study, the higher prediction accuracy for the multivariate models combining GY and PH confirmed that the correlation between traits is important for the success of multivariate models. A study with Sorghum also showed that genomic prediction for GY benefits mainly from using PH as a secondary trait (Velazco et al., 2019). The prediction accuracy increase varied across sites for GY, and the highest increase was at site LR, where the estimated heritability for GY and correlation between PH and GY were both low. This indicated that multivariate models could benefit more for the trait when estimated heritability is low. Up to 60% PA increase for low heritability traits was reported in a lentil GS study when implementing multivariate models (Haile et al., 2020). We also evaluated the prediction accuracy of different cross-validation scenarios. The scenario 2, PH only included in the training set (CV2), didn't improve the PA, indicating that adding information of correlated traits in the prediction dataset was essential for improving the prediction accuracy by the multivariate models. A similar study in barley reported that the correlated traits only included in the training population would not change the PA in validation and performed equivalently to the univariate model (Arojju et al., 2020). Plant height could be measured at the early stage, and high throughput phenotyping technology, such as the UAV-based RGB imaging system, has been reported to be used for plant height in wheat (Volpato et al., 2021). It indicated that early grain yield prediction by combining the high throughput PH measurements to improve the GS breeding efficiency could incorporate into future safflower GS breeding.

With the limited genetic and genomic resources for the safflower, taking advantage of the modern breeding tool, such as GS, could fast-track the genetic improvement of the crop. Tradition safflower breeding program depends on the phenotypic selection of the elite parental lines to initiate the crossing. In our study, we have chosen the top 20% of individuals according to the BLUEs and GEBVs, respectively, and the best gain achieved was 0.7 t/ha, about a 23% yield increase in the optimal site. This indicated the high yield potential for safflower genetic improvement. We expected to find the highest response by PS for within site selection because the selection differential was used to predict selection response in the next-generation. However, the selection responses in PS and GS were very similar, which could be associated with the fact that most of the top 20% of the accessions with the highest phenotypes had the highest GEBVs as well (Supplementary Table S4). Selection for rust resistance in wheat had similar observations that gain from GS was equal to PS when comparing the gain by two cycles of GS and one cycle of PS (Rutkoski et al., 2015). Bernardo (2020) recommended that selection based on breeding ‘lines’ performance (BLUEs) instead of GEBVs would be sufficient enough if the breeding goal is to pyramiding the elite alleles. However, with GS, selection decisions could be made in an early stage, for example, selection for grain yield before harvest. And when the breeding cycle increases, GS would benefit the breeding program significantly by reducing the breeding cycle. Those lines with high yield and high GEBVs in the safflower collection could be selected directly to achieve yield gain. Selection in the stressed environment, especially for drought resistance, has been reported to improve the selection in optimal conditions (Das et al., 2021; Kumar et al., 2021). However, we observed the selection in the RF and LR sites leads to the compromised gain for GY in IR. Although the gains from different strategies were high at the selected sites, the gap between the optimal and stressed environments was significant. The decreased genetic variation in the stress sites could be the reason for the lower selection response for GY in those sites (Xu et al., 2017). Falconer (1952) mentioned that the correlation between environments could impact the selection responses. In other words, g × E could reduce the rate of genetic gain achieved by breeding. The combined g × EGS selection increased gain slightly in RF and LR, indicating that combining g × E into the prediction model could benefit GS selection. Studies combined multi-environment trials in the genomic prediction model showed similar results (Gill et al., 2021; Khanna et al., 2022). In safflower, using suitable GS models for an optimized breeding pipeline, which selects genotypes that suit environments or genotypes with compromised performance but robust in all conditions, is essential.

Selection indices have been used efficiently for multi-trait breeding programs for animals and crops (Smith, 1936; Hazel, 1943). If selections were based mainly on GY, we observed that OL had limited positive gain in the optimal site and negative gains in stressed sites (Supplementary Table S2). The index selection achieved a high gain for OL and GY although the gain for GY was compromised. The loss of the gain for GY could be compensated by high OL, which would bring high economic values for the farmer. A study in Sorghum also indicated that GS with the optimal index selection to improve biofuel traits is the most promising strategy (Habyarimana et al., 2020a). An optimal breeding framework incorporating index selection has been proposed for different crops (Marulanda et al., 2016; Akdemir et al., 2019). As a minor crop, safflower faces similar challenges as other orphan crops, such as limited resources and limited research funds. GS with selection index is the fast way to pyramid elite alleles into the future safflower varieties. The optimal selection index for safflower GY and OL still requires further research. More information, e.g., biotic and abiotic tolerance, early vigour, seed protein content, etc., could be combined into the selection index. Further validation of different selection strategies, which would facilitate the optimization of breeding strategies by simulation, would also be required.

In summary, we estimated the heritability and genetic correlations as well as selection responses to different selection indices for a diverse safflower Genebank collection. We achieved higher genetic prediction accuracy for grain yield by combining the correlated trait into the multivariate model. High selection responses for GY were achieved at the selected sites with both PS and GS, but the gains varied across the sites. Combining the g×E model with the equal weighted selection index GS model could simultaneously improve GY and OL in all studied sites.