Combining Crop Growth Modeling With Trait-Assisted Prediction Improved the Prediction of Genotype by Environment Interactions

Plant breeders evaluate their selection candidates in multi-environment trials to estimate their performance in contrasted environments. The number of genotype/environment combinations that can be evaluated is strongly constrained by phenotyping costs and by the necessity to limit the evaluation to a few years. Genomic prediction models taking the genotype by environment interactions (GEI) into account can help breeders identify combination of (possibly unphenotyped) genotypes and target environments optimizing the traits under selection. We propose a new prediction approach in which a secondary trait available on both the calibration and the test sets is introduced as an environment specific covariate in the prediction model (trait-assisted prediction, TAP). The originality of this approach is that the phenotyping of the test set for the secondary trait is replaced by crop-growth model (CGM) predictions. So there is no need to sow and phenotype the test set in each environment which is a clear advantage over the classical trait-assisted prediction models. The interest of this approach, called CGM-TAP, is highest if the secondary trait is easy to predict with CGM and strongly related to the target trait in each environment (and thus capturing GEI). We tested CGM-TAP on bread wheat with heading date as secondary trait and grain yield as target trait. Simple CGM-TAP model with a linear effect of heading date resulted in high predictive abilities in three prediction scenarios (sparse testing, or prediction of new genotypes or of new environments). It increased predictive abilities of all reference GEI models, even those involving sophisticated environmental covariates.


INTRODUCTION
The objective of plant breeders is to develop varieties well adapted to target environments. For this purpose, they evaluate each year candidate varieties in multi-environment trials (MET). Given that the speed of the process is a key element of genetic progress and since phenotyping is expensive, most varieties are evaluated in a small number of environments, considered as a combination of year × site × condition. Consequently, the environments in which the varieties are evaluated can be quite different from the target environments, because of the significant variation between years. In addition, only a limited number of varieties are evaluated each year to control the phenotyping costs. All these constraints are reducing the chance of success as they limit the number of lines that can be evaluated in the target environments.
One way to raise this strong constraint is to predict the performance of candidate varieties using molecular information with genomic selection (GS) models (Whittaker et al., 2000;Meuwissen et al., 2001). In GS, a phenotyped and genotyped calibration set is used to estimate DNA marker effects. Once the model is calibrated, new candidate varieties can be predicted, as long as their genotypic information is available. Specific GS models were proposed to predict the performance of varieties in different environments, taking the genotype by environment interaction (GEI) into account. It was first proposed to adapt the reference GS models to the GEI context by attributing environment specific effects to the markers (Schulz-Streeck et al., 2013;Crossa et al., 2016), or by modeling genetic covariances between environments (Burgueño et al., 2012). In other studies, environmental covariates (EC) were introduced in the GS model (Heslot et al., 2014;Jarquín et al., 2014;Malosetti et al., 2016;Ly et al., 2018), which allows predicting the performance of varieties in new environments. Crop growth models (CGM) were sometimes used to adjust the EC estimates to phenological stages, or to derive EC estimating the stress experienced by the plants (Ly et al., 2017;Rincent et al., 2019), instead of directly using pedoclimatic data.
One efficient way to improve the prediction accuracy of GEI is to introduce secondary traits measured in each environment on both the calibration set and the test set in multi-trait GS models. It was indeed shown that if the secondary trait is sufficiently heritable and correlated to the target trait, the approach called "trait-assisted prediction" (TAP) or "phenotype imputation", can be very efficient (Henderson and Quaas, 1976;Calus and Veerkamp, 2011;Jia and Jannink, 2012;Pszczola et al., 2013;Fernandes et al., 2018;Bustos-Korts et al., 2019). It resulted in successful applications in wheat (Rutkoski et al., 2016;Sun et al., 2017Sun et al., , 2019Crain et al., 2018;Lado et al., 2018;Michel et al., 2018;Schulthess et al., 2018). In particular Rutkoski et al. (2016) showed that secondary physiological traits measured in each environment on both predicted and test sets resulted in increased predictive ability for grain yield using a TAP model. TAP is particularly useful when the secondary trait is easy to obtain, highly heritable and strongly correlated to the target trait. In MET, the secondary trait measured in each environment can capture GEI and serves as a proxy of the target trait in each environment (Bustos-Korts et al., 2019). But the applicability of this approach in the GEI context is considerably constrained by the necessity to phenotype all varieties (calibration and test sets) in all environments for the secondary trait, which means that the test set has to be sown and phenotyped in each environment. The interest of this approach would be considerably stronger if the secondary trait would be accessible without phenotyping the test set.
Crop Growth Models (CGM) are also powerful tools to predict GEI. They model plant development using mechanistic relationships with genetic characteristics (genetic parameters) and environmental variables as input. These genetic parameters characterize the varieties (e.g., sensitivity to photoperiod) independently from the environment, and so are supposed to be constant from one environment to another for a given variety (Reymond et al., 2003(Reymond et al., , 2004. A method combining CGM and GS, called gene-based modeling (GBM), can be used to predict GEI for unobserved varieties. If the genetic parameters of varieties of the calibration set are known, a GS model can be calibrated and used to predict the genetic parameters of the test varieties. These predictions can then be used as input for the CGM together with environmental variables to predict the performance of the test set in various environments.
This strategy has proven to be effective on simple traits such as leaf elongation rate in maize (Reymond et al., 2003;Chenu et al., 2008), fruit quality (Quilot et al., 2005;Prudent et al., 2011), and phenology of various species (White and Hoogenboom, 1996;Nakagawa et al., 2005;Yin, 2005;Messina et al., 2006;White et al., 2008;Uptmoor et al., 2011;Zheng et al., 2013;Onogi et al., 2016) including wheat (Bogard et al., 2014;Rincent et al., 2017). The major limit of this approach is that it remains difficult to efficiently predict complex traits such as yield, although promising results have been obtained (Technow et al., 2015;Cooper et al., 2016;Messina et al., 2018). In any case, GBM can be very efficient to predict secondary traits, for example related to phenology. If we were able to use CGM (or GBM for new varieties) to accurately predict a secondary trait correlated to the final target trait in each environment, then we would be able to combine CGM and trait assisted prediction to overcome the limits of both approaches. The idea here is to use a TAP approach in which the secondary trait is not phenotyped for the test set (so there is no need to sow the test set in any experiment) but predicted in each environment of interest thanks to CGM.
We propose to test this new approach (CGM-TAP) on winter wheat with grain yield (GY) as target trait (the trait the breeder is interested in), and heading date (HD) as secondary trait. HD is very highly correlated to flowering time, which is indeed an important adaptive trait with an optimal date depending on the environment. Phenology is a major trait for the plants to benefit from the most possible resources while avoiding stressing conditions at key stages (Richards, 1991;Lopes et al., 2014;Semenov et al., 2014;Flohr et al., 2017). In addition, it was shown that GBM was efficient to predict HD in winter wheat (Bogard et al., 2014;Rincent et al., 2017). Therefore, this trait seems to be a particularly good candidate to test the CGM-TAP approach with GY as target trait. For this purpose, we used a MET composed of 220 varieties and 42 managed environments, that was used to compare the accuracy of CGM-TAP to reference GEI models in different prediction scenarios: the prediction of new varieties, of new environments, or of the missing phenotypes of an incomplete design (sparse testing).

MATERIALS AND METHODS
Genetic Material, Genotyping, and Estimation of a Genetic Covariance Matrix (Kinship) The plant material has been previously described in Rincent et al. (2018Rincent et al. ( , 2019 and Touzy et al. (2019). The genetic panel is composed of 220 European elite varieties of winter wheat. It was genotyped with the TaBW280K high-throughput genotyping array described in Rimbert et al. (2018). This array was designed to cover both genic and intergenic regions of the three bread wheat subgenomes. Markers with a minor allele frequency (MAF) below 5%, or with heterozygosity or missing data rate above 5% were removed. Markers in strong Linkage Disequilibrium (LD) were filtered out using the pruning function of Plink (Purcell et al., 2007) with a window of size 100 SNPs, a step of 5 SNPs and a LD threshold of 0.8, as proposed in Charmet et al. (2020). Eventually, we obtained 26,116 polymorphic high resolution SNPs, with an average missing data rate of 1.0%. Missing values were imputed as the marker observed frequency.
Genotype of individual i at marker l (M i,l ) was coded as 1, 0.5, or 0 for homozygote for an arbitrarily chosen allele, heterozygote, and the other homozygote, respectively. Genomic relatedness (kinship) between individuals was estimated following (VanRaden, 2008): p l × 1 − p l , p l being the allelic frequency of the reference allele in the corresponding diversity panel, L the number of markers. The matrix of genomic relatedness coefficients will be denoted K in the manuscript.

Phenotypic Data Measured in the Multi-Environment Trial
The same phenotypes as in Rincent et al. (2019) were used. Briefly, the panel was phenotyped for grain yield (GY) and heading date (HD) in a multi-environment trial composed of 42 managed environments located in France between 2012 and 2016 (Supplementary Table S1). These 42 managed environments correspond to 26 combinations of years and locations, with two treatments for 16 of them. Among these 16 combinations of years and locations, three had an irrigated (WW) and a rainfed (WD) treatment, one had a well-watered (WW) and a rainout shelter (RO) treatments, and 12 had a high (HN) and a low (LN) nitrogen fertilizer treatments. These 42 combinations of location, year and management will be called environments in the rest of the document. Experimental designs, estimation of adjusted means and heritabilities were presented in Rincent et al. (2019). Briefly summarized, in each environment the varieties were grouped in six to eight blocks according to their earliness. The designs were two-replicate designs or augmented designs. Adjusted means and heritabilities were estimated with SpATS (Rodríguez-Álvarez et al., 2018) for the two-replicate designs to take spatial trends into account and with block effect only for the augmented designs. On this dataset, the use of weights in the second step of the analysis to take into account the difference of precision of the adjusted means (Damesa et al., 2017) did not improve the results, and so a basic analysis without weights was applied here.

Environmental Characterization and Estimation of an Environmental Covariance Matrix
Each environment was characterized by 139 environmental covariates (EC) in Rincent et al., 2019. Seventy-two of these covariates were estimated using climatic data (temperature, radiation), and 67 were estimated with SiriusQuality CGM (Martre et al., 2006) as dry matter stress index. These 139 covariates were estimated in each environment, transformed to standard normal distributions, and compiled in a matrix of dimension (42 × 139).
was then used to estimate a covariance matrix between environments (W). We applied the approach of Jarquín et al. (2014) in which all the ECs were used to estimate the covariance matrix. Environments with similar stresses are assumed to have similar GxE patterns. To compute W, we proceeded in two steps: first we computed the Euclidean distance matrix between environments (D ) with the matrix of environmental covariates ( ), and then the covariance matrix W was computed as: . W does not take into account the targeted response variable, but the EC are supposed to be helpful to give an estimate of the true unknown environmental covariance matrix, as they reflect the conditions experienced by the plants. W was used in the downstream analysis as an estimate of the environmental covariance matrix.

Selection of a Subset of Representative Environments
Because of the redundancy in the dataset due to strong similarities between the environments , and to limit computational burden, a subset of sixteen environments was selected. The selection was applied in such a way that the different kinds of HD/GY relationship were explored. The 16 environments were composed of four environments with a low linear correlation between HD and GY (absolute value of the correlation <0.3), four environments with an intermediate HD/GY correlation (absolute value of the correlation between 0.3 and 0.6), four environments with a strong HD/GY correlation (absolute value of the correlation above 0.6), and four environments for which a quadratic relationship explained more than a linear relationship ( Table 1). These four groups of environments were named "low, " "medium, " "high, " and "quadratic" according to the kind of relationship between HD and GY. The 16 environments were randomly sampled among the 42 environments of Rincent et al. (2019) with the constraint that each of the four groups of environments were represented by four environments. The prediction scenarios presented below were applied to these sixteen environments.

Prediction Objectives and the Corresponding Cross-Validation Schemes
Three prediction objectives were considered: the prediction of observed varieties in observed environments (oGoE) or in new environments (oGnE), and the prediction of new varieties in observed environments (nGoE). oGoE consists of predicting missing values in a MET (sparse testing), which typically corresponds to the situation faced by breeders when some observations are missing in their trial networks. oGnE and nGoE are more ambitious because predictions are made in an environment or for a variety without any phenotypic information on it. a The sixteen environments were grouped in four classes according to the kind of relationship between HD and GY: "low," "medium," or "high" correlation, or "quadratic" relationship. b α estimates in the 16 environments using model EG_HD. c α and d β estimates in the 16 environments using model EG_HD 2 . α and β correspond to the coefficients for the linear and quadratic regressions of GY on HD in each environment. The estimates of α and β presented here were obtained with the full dataset.
To evaluate the performance of the prediction models in these three situations, three cross-validation schemes were defined: CVrandom, CVnewG, and CVnewE.
Observed varieties in observed environments was addressed by the CVrandom scheme, that is a 6-folds cross-validation, in which the folds were randomly sampled from the dataset. oGnE was addressed by the CVnewE scheme, which is a leaveone environment-out scheme, in which a new environment is predicted. nGoE was addressed by the CVnewG scheme, with a division in six folds consisting in five randomly selected groups of varieties.
Predictive abilities were computed for each fold and each environment as the correlation between adjusted means and predictions in the test set. For each cross-validation scheme except CVnewE (leave-one-out scheme), the total procedure was repeated 10 times to get robust estimates of predictive abilities.

Reference Prediction Models
Three reference models with various levels of complexity were used to predict GY in the MET. In the first kind of models, the kinship matrix allowed sharing information between varieties, but there was no sharing of information between environments: I N E is an identity matrix of size the number of environments and ⊗ is the Kronecker product. µ j is a fixed environmental effect, G i is a random polygenic effect for variety i, and GE ij is a random effect corresponding to the interaction between variety i and environment j. In all models, ∈ ij ∼ N(0, Iσ 2 e ), is the residual effect. In model EG the prediction of a variety is the same for any environment (main effect model). In model EG_GxE, GEI are taken into account, but there is no sharing of information between environments. Consequently, this second model is only applicable to scenarios CVrandom and CVnewG.
In the third model, information between environments is shared thanks to an environmental covariance matrix estimated with the environmental covariates as proposed by Jarquín et al. (2014): W was estimated as described above.

The Crop-Growth Modeling -Trait Assisted Prediction Model (CGM-TAP)
Crop-growth modeling -trait assisted prediction is a trait assisted prediction model in which the secondary trait is observed for the calibration set but predicted for the test set using crop growth modelling (CGM) or gene-based crop-growth modelling (GBM) (Figure 1). The idea is to combine CGM and TAP to overcome the limit of GBM (not applicable yet to complex traits such as yield) and TAP (requires the phenotyping of the test set for the secondary trait in each environment). First (step 1), the CGM or GBM approach is applied to predict the secondary trait (HD) for FIGURE 1 | Schematic representation of CGM-TAP. CGM-TAP is a trait-assisted prediction approach in which the secondary trait is predicted using crop-growth modeling instead of being phenotyped.
Frontiers in Plant Science | www.frontiersin.org the test set in each environment. Then (step 2), a trait-assisted prediction model is applied to the full dataset with the secondary trait as covariate (observed for the calibration set, predicted with CGM or GBM for the test set) with an environment specific effect to predict the target trait (GY) in each environment. In more details, CGM-TAP consists in the following two steps: - Step 1: Prediction of the secondary trait (HD) for the test set using CGM. In this first step, the two main genetic parameters of the CGM involved in phenology (sensitivity to photoperiod "SLDL" and phyllochron "P") are estimated (CVnewE and CVrandom) or predicted (CVnewG) for the test set. These variety specific genetic parameters are then used as input to the CGM together with the daily climatic data (temperature, daylength) to predict HD for the test varieties in each environment. In scenarios CVnewE and CVrandom, the test varieties are scored for HD in some environments, and so a Bayesian algorithm  was used to directly estimate the genetic parameters of the test varieties. Briefly, the MCMC algorithm is a hybrid Gibbs sampler that updates the coordinates of the different parameters through a Metropolis-Hastings step using as proposal a Gaussian distribution centered on the previous value of the chain. In scenario CVnewG, there is no phenotypic information on the test varieties and so the genetic parameters cannot be estimated. So, a GBM approach was applied in which the genetic parameters are estimated for the calibration varieties, and then a GBLUP model (with major phenology markers as fixed effects) was used to predict the genetic parameters of the test varieties. -Step 2: Different trait-assisted prediction models were then used to predict the target trait (GY) for the test set: with h ij the heading date of variety i in environment j. h ij was observed for the calibration varieties and predicted for the test set. α j is the effect of HD on GY in environment j. Three other models similar to models 4-6 but with a quadratic effect of HD were also used: with h 2 ij the squared heading date (centered and scaled) of variety i in environment j. β j is the quadratic effect of HD on GY in environment j. All models are summarized in Table 2.
For CVnewE, α j and β j cannot be estimated in models 4-9 because there is no observation of HD and GY in the test environment even for the calibration set, and so an alternative procedure was applied. First, α j and β j were estimated using a linear or a quadratic regression of GY on HD in each environment except the test environment. The estimates of the regression coefficients were then used to fit a multiple linear regression on the EC measured in these environments. A stepwise forward-backward procedure (function "step" of R package "stats" with a penalty k = 3.7 for α and k = 8 for β) was applied to select the most relevant EC and estimate their effects. This calibrated model was then used to predict α j and β j in the test environment. The estimatesα j andβ j were plugged into models 4-9, which could then be used to predict GY in the test environment. We suppose here that the target environments are characterized by EC, so they are not totally unknown.
For each model, predictions were obtained as the sum of the BLUEs of the fixed effects and the BLUPs of the random effects. For each cross-validation scenario, the predictive ability (correlation between the predictions and the adjusted means) and root mean square error (RMSE) of HD in step 1 was estimated, as well as the RMSE of α j and β j for CVnewE. The predictive ability of GY in step 2 was evaluated for each model as presented above.

Relationship Between HD and GY in the Sixteen Environments
The variance of HD was highly variable from one environment to another (Figure 2). HD was constraint to 15 days in Vra13LN whereas it covered around 30 days in other environments (Gre14WD, Coi12WW).
It is specified if the effect was included (x) or excluded (-) from the model. For mixed effects (G and GEI) the structure of the covariance matrix is indicated.
Frontiers in Plant Science | www.frontiersin.org The average adjusted R 2 of the linear relationship was of 0.04, 0.26, 0.46, and 0.13 in the low, medium, high and quadratic classes, respectively. For the quadratic environments, the adjusted R 2 of the quadratic regression was of 0.27 on average, and higher than the adjusted R 2 of the linear regression. In all four environments of this class, the quadratic relationship was concave with an optimal HD.

Predictive Ability and RMSE of the Secondary Trait (HD) in the Different Scenarios
Before looking at the predictive abilities of GY with the different models, we looked at the predictive ability of the secondary trait (HD) in the different cross-validation scenarios ( Table 3). The average predictive abilities of HD were above 0.78 in all scenarios, and it was higher in scenarios CVrandom and CVnewE (0.91) than in CVnewG (0.78). This was expected as the genetic parameters of the CGM had to be predicted in scenario CVnewG. RMSE was low in the three scenarios (average always below 3.33 days), and it was lower for CVrandom (1.86 days) than for CVnewG (2.88 days) and CVnewE (3.33 days). There was a strong variability of RMSE between folds with a maximum of 6.63 days in scenario CVnewE.

Prediction of the Regression Coefficients of HD on GY (α and β) in Scenario CVnewE
For models EG_HD, EG_GxE_HD, and EG_GxW_HD, it was necessary to first predict the value of α to be able to run the models in CVnewE, because no phenotypic observation is available on the predicted environment and thus, the relationship between HD and GY could not be directly estimated. For the same reasons, both α and β should be predicted before running models EG_HD 2 , EG_GxE_HD 2 , and EG_GxW_HD 2 .
The multiple linear regression on the environmental covariates was moderately efficient to predict the value of α for the linear regression of HD on GY in the new environment ( Figure 3A). The RMSE was equal to 0.40 for observed α values varying between −0.99 and 1.91. For all environments except one, the environmental covariates were able to accurately predict the sign of the regression.
The predictions of α and β were less accurate for the quadratic linear regression of GY on HD, with RMSE equal to 0.44 for α, CVrandom corresponds to sparse testing, CVnewG to the prediction of unobserved varieties, and CVnewE to the prediction of unobserved environments. CVrandom and CVnewG were 6-folds cross-validations repeated ten times, and CVnewE was a leave-one-out scenario. Predictive abilities were computed as the correlation between adjusted means and predictions. RMSE are the root mean squared errors. and 0.50 for β (Figures 3B,C). For β, the sign of the regression was inaccurately predicted for three environments.

Predictive Ability of the Target Trait (GY) in the Different Scenarios
The predictive abilities for the target trait (GY) were highly variable between scenarios and between models ( Table 4). On average over the four kinds of environments (low, medium, high and quadratic), the highest predictive abilities were obtained in scenario CVrandom (0.73 for EG_GxW_HDpred and EG_GxW_HDpred 2 ). It was slightly lower for CVnewE (0.70 for EG_GxW_HDpred) and CVnewG (0.57 for EG_GxW_HDpred). The inclusion of GEI effects in the model (GxE and GxW) always increased predictive abilities in comparison to the main effect model EG, in particular when EC were used to estimate covariance between environments (GxW). On average over the four kinds of environments, the introduction of term GxW increased predictive abilities by 39% for scenario CVrandom, 37% for CVnewG, and 25% for CVnewE in comparison to the main effect model. For all scenarios and kinds of environment, the best model was always a CGM-TAP model, except in CVnewE for "low" and "quadratic" environments. The inclusion of HD as fixed effect (CGM-TAP) in the three reference models (EG, EG_GxE, and EG_GxW) always increased the average (over the four kinds of environments) predictive abilities. This increase strongly depended on the kind of environment and on the reference model considered. For the environment with a low correlation between GY and HD, and the environments with a quadratic relationship, the inclusion of HD as fixed effect in the reference models resulted in similar predictive abilities than those of the reference models, or to lower predictive abilities in CVnewE. Inclusion of HD increased predictive abilities of the reference models for environments with a medium GY/HD relationship. This increase was of 31, 30, and 5% in CVrandom, CVnewG and CVnewE, respectively, for model EG. As expected, the interest of CGM-TAP was highest for environments with a strong GY/HD relationship ("high"). For these environments the increase brought by CGM-TAP model was of 141, 127, and 132% in comparison to model EG in CVrandom, CVnewG and CVnewE, respectively. This increase was high but less pronounced for more complex reference models: it was of 7.1, 5.3, and 29.0% in comparison to model EG_GxW in CVrandom, CVnewG and CVnewE, respectively. Interestingly, in many cases, the simple CGM-TAP models (EG_HDpred, and ED_GxE_HDpred) performed as well as complex model involving environmental covariates (EG_GxW).
The introduction of a quadratic relationship in the reference models only slightly increased predictive abilities in scenario CVrandom is a 6-fold cross-validation scheme in which the test set was randomly sampled from the full dataset (sparse testing). CVnewG is a 6-fold cross-validation scheme, for which the objective is to predict new varieties in known environments. CVnewE is a leave-one environment-out scheme, in which the objective is to predict the performance of known varieties in a new environment. "low," "medium," "high," and "quadratic" correspond to groups of four environments with different kinds of HD/GY relationships ( Table 1 and Figure 2). No results are indicated for models involving the GxE term in scenario CVnewE, because these models are unable to make predictions in environments without any phenotypic data.
CVrandom. In other scenarios, the introduction of a quadratic relationship resulted in similar or lower abilities. It considerably decreased predictive abilities in scenario CVnewE, in which the regression coefficients (α and β) have to be predicted first using the EC (there is no phenotypic information in the test environments, and so the regression coefficients of the GY/HD relationship cannot be directly estimated).

DISCUSSION
In plant breeding, the knowledge that we have on the candidate varieties is limited because of the necessity to control phenotyping costs and the duration of the breeding cycles. This lack of information considerably limits the chance of identifying top varieties candidate for the registration process. Genomic prediction is a promising tool to screen early generation material (Heffner et al., 2010), because it allows transferring the information collected on the previous breeding cycles to other combination of genes and environmental conditions even in presence of GEI. Different models were proposed in the past years to predict GEI, which involved genetic and environmental characteristics (markers and EC) (Heslot et al., 2014;Jarquín et al., 2014;Cooper et al., 2016;Ly et al., 2018). We developed a new approach, called CGM-TAP, in which a secondary trait is predicted for the test set in each environment thanks to CGM. This secondary trait, easy to predict and capturing GEI, is then used as an environment specific covariate in usual GEI GS models. To sum up, it is a trait assisted prediction model, in which the secondary trait is predicted thanks to CGM instead of being phenotyped. Here we used HD and GY as secondary and target traits, respectively, to evaluate this new approach.
The results presented here reveal that in almost all situations, the CGM-TAP models performed at least as good as the corresponding reference models. It performed even better than the most sophisticated reference model involving EC (EG_GxW) with an increase of predictive ability of 7.1, 5.3, and 29.0% for the environments with the highest correlation between HD and GY ("high") in scenarios CVrandom, CVnewG and CVnewE, respectively. The increase of predictive ability was generally stronger when the relationship between the secondary and the target trait was strong, as expected. But even for the "low" environments, CGM-TAP models (without the quadratic relationship) performed at least as good as the reference models. In the three cross-validation scenarios there were even situations in which the simple CGM-TAP model including only the main effects and the secondary trait (EG_HD) performed as good as the most sophisticated reference model EG_GxW (Table 4).
However, the introduction of a quadratic effect in the CGM-TAP models did not lead to a higher predictive ability and was even detrimental for the "quadratic" environments in scenario CVnewE. This is because we were unable to predict the regression coefficients α and β sufficiently well in new environments (Figure 3). The use of a larger set of environments would probably help to better predict these coefficients. Another option would be to use CGM to directly predict the HD/GY relationship in new environments, but we were unable to do it accurately. In other prediction scenarios (CVrandom and CVnewG), the introduction of the quadratic effect resulted in similar or lower predictive abilities in comparison to the CGM-TAP models without quadratic effect. This may be due to the necessity to reach more precise prediction of the secondary trait if a quadratic relationship is fitted. Another explanation is that the quadratic relationships were not very pronounced in this dataset (Figure 2).
These promising results show that the predictions of the secondary trait by the CGM were sufficiently accurate to be able to run a trait-assisted prediction model (without quadratic relationship) in all three prediction scenarios without phenotyping the test set. The RMSE of the prediction of HD with the CGM were indeed relatively low (Table 3). It was equal to 1.86 days in scenario CVrandom, which is close to usual field phenotyping error for this trait. The interest of the CGM-TAP models is dependent on our ability to predict the secondary trait. This means that CGM-TAP should work better if the calibration set was bigger, or was phenotyped in a MET optimized for the genetic parameters estimation . Note that the use of the observed HD instead of the predicted HD as covariate in the CGM-TAP models resulted in much higher predictive abilities (results no shown). This means that CGM-TAP could work much better if we were able to predict HD more accurately.
Most studies on multi-trait and trait-assisted predictions are based on mixed models involving a genetic correlation between secondary and target traits (Rutkoski et al., 2016;Sun et al., 2017;Michel et al., 2018;Lado et al., 2018;Schulthess et al., 2018). Here, we proposed to introduce the secondary trait as covariate as in Crain et al. (2018), which allows modeling non-linear relationship between traits. This is useful in some situations with complex dependencies as illustrated by the quadratic relationship between HD and GY in some environments. This could also be the case with other secondary traits, for which there is an optimum value depending on the environment such as aerial or root biomass.
The choice of HD as a secondary trait was justified by the fact that it is easily predicted by CGM (Bogard et al., 2014), and it is an adaptive trait, and so is expected to be related to GY, at least in some environments (Richards, 1991;Semenov et al., 2014;Flohr et al., 2017). The relationship between phenology and productivity is important for most species (Andrade, 1995;Ouk et al., 2007), and so we can expect our strategy to be valuable for other crops. This relationship was confirmed in our dataset in which environments with strong HD/GY correlations or quadratic relationship could be found (Figure 2). However, the range of HD in breeding programs would be narrower and the correlation probably less pronounced. HD was a good candidate to test this approach, but other traits even more related to the target trait could be used. The important development of highthroughput phenotyping and of crop-growth modeling should make it possible to predict other yield related traits in the future. A suitable secondary trait to run CGM-TAP models should be accurately predicted with CGM and strongly related to the target trait. Traits related to light interception (LAI) or even yield components are secondary traits of choice and may be accurately predicted thanks to CGM in the near future with the help of HTP.
Note that the predictive abilities obtained in scenario CVnewG are probably higher in this study as would be in practice, as the folds of the cross-validation were randomly sampled. This indeed implies that varieties in the calibration set and in the test set are sampled from a same population, which is not the case in a real breeding scheme. Considering the reference models, the results of this study confirm that the introduction of an interaction (GEI) term in the prediction model increased predictive abilities in comparison to a main effect model for all scenarios. The model involving the environmental covariance matrix W (EG_GxW), performed better than model EG_GxE in scenario CVrandom, but similarly in scenario CVnewG. A major advantage of EG_GxW is that it allowed predicting in new environments, and these predictions were much more accurate than the main effect model. This is in agreement with previous publications (Jarquín et al., 2014;Malosetti et al., 2016;Rincent et al., 2019).

CONCLUSION
In this paper, we propose a new way to predict genotype by environment interactions (CGM-TAP). This is a trait-assisted prediction approach in which the secondary trait is predicted by a crop growth model, instead of being phenotyped. The relationship between the secondary trait and the target trait is environment specific and thus allow predicting environment specific effects. This approach applied to yield (target trait) and heading date (secondary trait) increased predictive abilities of all reference GEI models, even those involving sophisticated environmental covariates, for various prediction scenarios (new varieties, new environments and sparse testing). This proof of concept could be applied to other traits in breeding programs in the near future, and is a new way of combining crop growth modeling and genomic prediction in the context of GEI.

DATA AVAILABILITY STATEMENT
The data supporting the conclusions of this article can be found at https://doi.org/10.15454/TKMGCQ.

AUTHOR CONTRIBUTIONS
RR defined the prediction strategy and prediction models, ran the analysis, and wrote the manuscript. PR ran the analysis and reviewed the manuscript. JL reviewed the manuscript, participated in the discussions, and contributed to the data production with the BreedWheat Consortium. All authors contributed to the article and approved the submitted version.

FUNDING
Management of the wheat multi-environment trials was financially supported by the French National Research National Agency under Investment for the Future (BreedWheat project ANR-10-BTBR-03) and by FranceAgriMer. The Phéno3C platform was financially funded by the French National Research National Agency under the Investment for the Future Phenome project (ANR-11-INBS-12) and by the European Regional Development Fund (AV0011535). We thank H2020 INVITE for financial support (grant agreement no. 817970).

ACKNOWLEDGMENTS
We thank three reviewers for their relevant comments which improved the manuscript.