The phenotypic data used in this study are composed of 136 sorghum accessions selected within the Generation Challenge reference set (Billot et al., 2013). Phenotyping was conducted on the PhenoArch high-throughput phenotyping platform in Montpellier, France (Cabrera Bosquet et al., 2015) where all genotypes were followed for 45 days from September to October 2017. The Phenoarch platform is based on a PhenoWare™ system (Lyon, France) composed of a conveyor belt structure of 28 lanes carrying 60 carts of one pot. Plants were grown in polyvinyl chloride (PVC) 9L pots (0.19m diameter and 0.4m high) filled with a 30:70 mixture of a clay and organic compost. Line spacing was of 0.4m and row spacing was of 0.2m. Five grains per pot were sown on August 30^th, which were subsequently thinned to a single plant two weeks after sowing. Two water treatment scenarios were tested during the late vegetative phase (starting from eight fully expanded leaves): well-watered (WW) plants were kept at 132% of soil humidity (Fraction of Transpirable Soil Water, FTSW of 0.6) and water-deficit (WD) plants were dried-down until 60% of soil humidity (± 5%, FTSW of 0.22).

Within each treatment, genotypes were replicated four times and local interactions were reduced by separating the genotypes into seven classes depending on plant height. To control the environmental heterogeneity of the greenhouse, the 28 rows were separated into four complete blocks defined according to the environmental gradient of the greenhouse (light and temperature, Cabrera-Bosquet et al., 2016). Each block was further divided into seven sub-blocks wherein plant height classes were randomly assigned following a Youden square (28 genotypes per sub-block). A d-optimal design generation software (SAS procedure OPTEX) was used to assign the genotypes to the sub blocks following an alpha-lattice.

All along the experiment, day temperature was kept at 30°C and night temperatures at 23°C. Environmental variables, air temperature (°C), radiation (PFFD, µmol m-2 s-1) and air relative humidity (%) were continuously monitored at eight positions in the greenhouse, above the canopy and recorded every fifteen minutes.

Thirteen pictures (twelve sides and one top) were taken daily of each plant. Seven genotypes (one for each plant height class, at three different growth stages) were used in a side experiment to calibrate a model used to estimate the following traits on all plants: aboveground fresh weight (Biomaerofw), plant leaf area (PLA), and plant height (PHT). Weekly measurements by hand of additional traits were conducted to assess: the number of appeared (App) and ligulated (Lig) leaves on the main stem and the number of tillers (Tillernb). After 45 days, plants were harvested and final measurements were conducted: shoot (Biomaerofw), mainstem (Mainstemfw) and mainstem blade (Bladefw) biomass fresh weight, base mainstem diameter (Stemdiam), and last ligulated leaf length and width (used to compute area of last ligulated leaf, Arealfel). Table 1 summarizes all measured traits.

The integrative trait considered for genomic prediction in this study is a Best Linear Unbiased Predictor (BLUP) of above-ground biomass, estimated by a mixed model (Equation 1, modeled using the R package Asreml v3, Butler et al., 2009) and considering all 8 plant replicates (4 replicates per water treatment, two water treatments).

With a i the fixed effect of the replicate, b the fixed competition effect of any neighbor, n the number of neighbors, C j the random effect of any row j (except the first and second row), d 1 and d 2 the fixed effect of the first and second row. 𝟙 j = 1 and 𝟙 j = 2 the indicators of the first and second rows: their value is 1 if the pot belongs to the first (respectively second) row, and 0 otherwise. G_v the random genotypic effect of variety v, w_t the fixed effect of the watering treatment t and (Gw) _vt their interaction. H_ik the between sub-blocks error, i.e. random effect of the sub-block of the i ^th replicate to which the k ^th height class was assigned. E_ivt the residual error for variety v in replicate i of treatment t.

This mixed model was the result of a more complete analysis of the data and assumes that it partially compensates for the heterogeneity of the greenhouse by modelling it as row and column effects. The mixed model was selected after exploring several types of variance decomposition. For the estimation of genotypic parameters of the crop model, raw data from the four WW replicate plants were used alongside local weather data (see Cabrera-Bosquet et al., 2016), as environmental effects are formalized inside the crop model (see Ecomeristem Model).

Genotypic information of the 136 accessions was obtained through Genotyping by Sequencing (GBS). Sequencing libraries were prepared according to the GBS protocol as per Elshire et al. (2011), with the ApeK1 enzyme. Single-end sequencing was performed on an Illumina HiSeq2000 (at the Genotool platform in Toulouse, France). SNP calling was performed using the GATK pipeline. After the SNP calling step, imputation was then realized using Beagle v.4 (Browning and Browning, 2013) with a 1000-SNP window and an overlap of 500 SNP after filtering on missing data per SNP (60% maximum of missing data and 5% maximum for minor allele frequency to keep a SNP). The imputed genotype matrix available for further analysis contained 31 713 SNP on the whole set of 136 accessions.

Genomic prediction was performed either by using the Least Absolute Shrinkage and Selection Operator (LASSO) implemented in the glmnet R package (Friedman et al., 2010), Bayes C implemented in the BGLR R package (Perez and de los Campos, 2014) or by a multilayer Convolutional Neural Networks (CNN) constructed with Tensorflow 2.0 (Singh and Manure, 2020) in Python (Van Rossum and Drake, 2009). For each method, validation was performed using a k-fold cross-validation method (with k = 5). The composition of each fold was equal across methods.

Two scenarios were investigated in this study. In the first instance, the direct prediction of aboveground fresh weight was done using the abovementioned three regression methods. The 31 713 SNP were used as the explanatory variables and the observed Biomaerofw as the explained variable. In the second scenario, the CGM Ecomeristem was used. First its input parameters were predicted using the same three regression methods, i.e. the SNP were used as the explanatory variables and the eight genotypic parameters as the explained variables. Then the predicted parameters were used to simulate growth and development of each genotype. The Ecomeristem output variable “Biomaerofw” was used as prediction for the aboveground fresh weight and compared to the observed Biomaerofw. The two scenarios and the three regression methods were compared based on the Normalized Mean Absolute Error (NMAE, Equation 3) between predicted and observed phenotypic trait, as well as the predictive ability (PA, Equation 4).

Where y ^ i are the predicted values of either the integrative phenotypic trait or the CGM parameters, y i the observed values of either the trait (phenotyped on the Phenoarch platform) or the CGM parameters (estimated trough the DE optimization algorithm, see Parameter estimation of the CGM), σ the standard deviation and C o v the covariance.

The genetic markers used as input data were encoded as {-1, 0, 1} corresponding respectively to homozygous for the reference allele, heterozygous and the homozygous for the alternative allele.

Narrow-sense heritability for all phenotypic traits, and CGM parameters, were computed using the genotypic and phenotypic data described above. In a first step, the Genomic Relationship Matrix (GRM) using the method proposed by VanRaden (2008) and implemented in the snpReady R package (Granato et al., 2018) was generated. Then, marker-based heritability was estimated using the heritability R package (Kruijer et al., 2015). The results are presented in the Supplementary Table S2 .

The first scenario compared the three regression methods on the direct prediction of sorghum aboveground fresh weight ( Figure 3 ). LASSO yielded the least accurate results with a NMAE of 0.41 and a PA of 0.39. Bayes C had slightly better results with a NMAE of 0.4 and a PA of 0.47. Finally, genomic prediction using a CNN showed the best results with a NMAE of 0.22 and a PA of 0.53.

The genomic prediction of CGM parameters yielded variable results across predicted parameters but not across prediction methods ( Table 3 ). NMAE varied from 0.87 for the least accurately predicted parameter, SLAp, to 0.05 for the most accurately predicted parameter, Ligulo_init. The average NMAE for each method was of 0.26 (SD = 0.27) for LASSO, 0.24 (SD = 0.24) for Bayes C and 0.24 (SD = 0.24) for CNN.

The prediction of Biomaerofw using Ecomeristem ( Figure 4 ) was slightly more precise when using parameters predicted with a CNN (NMAE 0.19, PA 0.63). LASSO and Bayes C still yielded parameters that resulted in good predictions (NMAE 0.20 and 0.19 respectively and PA 0.52 and 0.61).

Biomass components (secondary traits) were also simulated using Ecomeristem, based on the prediction of genotypic parameters through a convolutional neural network and are presented in Supplementary Figure 1 . Mainstem biomass (Mainstemfw, NMAE = 0.18, r = 0.72) as well as plant height (Pht, NMAE = 0.18, r = 0.6) were accurately predicted. The predictive ability for both individual leaf size (Arealfel, NMAE = 0.13, r = 0.53) and leaf biomass (Bladefw, NMAE = 0.22, r = 0.47) were lower.

Integrative traits are under the influence of (i) many genes and genetic interactions, especially epistasis (Zuk et al., 2012), and (ii) Genotype x Environment interactions (Rincent et al., 2017). These interactions are, by nature, non-linear. Therefore, to assess the capacity of linear regression models in producing accurate predictions of these traits, LASSO and Bayes C were compared to a non-linear regression method, specifically a convolutional neural network (CNN). The results revealed that LASSO exhibited the least accurate predictions for the studied integrative trait. While Bayes C demonstrated higher prediction accuracy, it still faced challenges in predicting aboveground fresh weight compared to convolutional neural networks (almost twice as high NMAE, Figure 3 ). These results are in line with previous studies comparing LASSO and Bayesian methods (e.g. Howard et al., 2014).

The use of non-linear regression methods, capable of efficiently considering interactions between polymorphism sites, could potentially address the challenge of epistasis in predicting integrative traits. Results showed that CNN, as a non-linear regression method, surpassed both LASSO and Bayes C in predicting total biomass fresh weight. The complexity of sorghum biomass growth, influenced by numerous interdependent biological processes and characterized as polygenic (Habyarimana et al., 2020), underscores the relevance of considering epistatic interactions. Previous studies have shown the influence of epistasis on complex trait variation, such as growth rate or crop yield (Kroymann and Mitchell-Olds, 2005; Melchinger et al., 2007). Given that integrative traits are influenced by myriad of local and distal SNP-SNP interactions (Qian et al., 2017), convolutional neural networks, by design, emerge as a suitable tool for accommodating these epistatic interactions in the prediction of these traits. Recent studies proposed the application of neural networks (NN), and in particular CNNs in genomic prediction of complex human and animal traits (e.g. Pook et al., 2020). These studies show that NNs can exhibit either slightly better, similar or inferior prediction accuracies compared to classical linear regression. Notably, the outcome is highly dependent on the set of SNPs, the genetic architecture of the trait (e.g. Bellot et al., 2018), and the architecture of the NN itself (Zhu et al., 2021).

The use of NNs in genomic prediction is obviously case-dependent, and it is apparent that more research is needed to adapt these non-linear methods to the Genomic Selection process. Specifically, this involves understanding how genetic data is processed by the NN, and assessing their impact on prediction outcomes (e.g. Verplaetse et al., 2023). In contrast to classical linear regression, a universal neural network architecture that suits all cases does not exist. Apart from estimating regression parameters, the entire structure (including the number, type, and size of layers, activation functions, etc.) must be adapted to the specific prediction problem at hand.

The Ecomeristem crop growth model was used to provide a more comprehensive consideration of causal non-linear environmental effects, as well as interactions between biomass components on biomass growth, potentially leading to enhanced prediction accuracies. Biomass is a highly polygenic trait, Habyarimana et al. (2020) detected significant marker-trait associations across eight of the then sorghum chromosomes. It is also under the influence of epistatic interactions: Brown et al. (2008) highlighted epistatic interactions between two major dwarfing QTL, Ishimori et al. (2020) showed the important role of epistasis for total biomass as well as stem length. Moreover, biomass components show differences in the dynamics of the effects of different genes depending on environmental conditions as well as developmental stage (Mu et al., 2022).

In this context, we proceeded to predict genotypic parameters that govern the equations formalizing biological processes where GxE, as well as interactions between resulting intermediate traits take place. Not all CGM parameters where predicted with the same accuracy. The least accurate predicted parameter was associated with the control of specific leaf area decrease between successive leaves (SLAp), with an average NMAE of 0.81. This can be attributed to the limited variability observed on this parameter across genotypes, with 40% of them having identical values, although SLAp exhibited high heritability (0.616). The challenge in accurately predicting SLAp emphasizes the influence of parameter characteristics and parameter estimation quality on the efficiency of genomic predictions of crop model parameters. Variations in prediction accuracy across other parameters can be explained by their heritability (see Supplementary Table S2 ). Parameters such as Plasto_init, Phyllo_init and Ligulo_init, characterized by a lower heritability, resulted in less accurate predictions. Conversely, parameters with higher heritability, including MGR_init, LL_to_INL, and ICt, exhibited predictions that are more accurate.

While errors are not directly comparable between parameters, intermediate traits, and the integrative trait, the simulation of Biomaerofw with Ecomeristem showed higher accuracies than the direct prediction of the trait, as illustrated in Figures 4C and D . This observation aligns with recent studies (e.g. Heslot et al., 2014; Tolhurst et al., 2022; Filho et al., 2023) which highlighted the enhanced prediction accuracies achieved through the incorporation of environmental effects. Interestingly some studies, like the one of Widener et al. (2021), suggest that the inclusion of environmental covariates may have limited impact, if any, on improving predictions in extreme environments. Moreover, they suggest that the composition of the calibration set of environments may hold greater significance as they found that only a subset of the available environments was needed to accurately predict GEBV. In addition, Rogers and Holland (2022) showed that the environmental similarity between training and test sets had a great impact on prediction accuracies. Phenomic selection (PS, Rincent et al., 2018) is one way of considering the impact of GxE. It is suggested that environmental variation is captured by the spectra and could thus enhance prediction accuracies compared to genomic prediction (e.g. Lane et al., 2020; Robert et al., 2022). On the other hand, the integration of crop models, as proposed in our study, introduces promising avenue to include causal environmental effects into GP. By providing causal relationships between environmental variation and individual traits, crop models have the potential to offer insights of the environmental effect on phenotypic plasticity. For example, the Ecomeristem model has already demonstrated its ability to predict plant behavior in new environments, i.e. not used during the estimation of genotypic parameters, or for alternate crop management (Larue et al., 2019). In addition, CGM simulate integrative traits as a result of interactions between intermediate traits and as responses to plant state. Therefore, the integrated CGM-GP approach could also answer the varying effects of genes depending on plant developmental stage. However, further research is imperative to thoroughly assess and quantify the impact of the integration of crop models with genomic prediction on prediction accuracy of the genetic values of candidates to selection, particularly in the context of multi-environmental trials. Our dataset was composed of a single environment. While GxE was still introduced by considering individual plant micro-environment inside the greenhouse (see Phenotypic data section), it is still necessary to validate this approach in a multi-environment setting. Several studies have explored this question (e.g. Reymond et al., 2003; Technow et al., 2015; Jighly et al., 2023). Nevertheless, it is important to note that their methodologies differ from the approach proposed in our study, particularly in terms of how crop model parameters are predicted through the use of genomic prediction: in their approach, genotypic parameters are sampled in each iteration of a Bayesian model and used for predicting observed traits through the CGM that is then used to update the marker effects. They are thus estimating parameters at the same time as the markers effects. In our approach, the parameter values are considered to be known. More specifically, the parameters are estimated in a preliminary step by model inversion using an optimization algorithm (here, a metaheuristic called Differential Evolution). These “observed” parameter values are then used to train a prediction model using marker information with no feedback of the CGM performance, relying solely on the difference between predicted and “observed” parameter values. Once trained, the prediction model is then used to predict the value of these genotypic parameters for the validation set and are then fed to the CGM to predict the phenotypic traits. Furthermore, these studies often focus on either a limited number of markers or crop model parameters; or they consider a reduced set of physiological processes within the crop model. The differences in methods and areas of interest highlights the need for a comprehensive study into the potential benefits of employing crop models for predicting integrative traits.

As mentioned previously, the prediction accuracy of the CGM-GP approach is highly dependent on the calibration set, the quality of parameter estimation and also the relevance of the crop model (Rincent et al., 2017). Indeed, the main limiting factor of the CGM-GP approach is how well the crop model is able to simulate the integrative trait of interest. In our study, the “optimal” parameter set yielded an average NMAE of 0.065 on Biomaerofw ( Supplementary Table S1 ). A perfect prediction of the CGM parameters could not yield a more precise result than the “optimal” parameter set. In any case, our study shows that the CGM-GP approach improves the predictive accuracy of an integrative trait compared to its direct prediction.

In this study, we showed that using non-linear regression methods (convolutional neural networks) in the genomic prediction process yielded higher prediction accuracies than classical linear methods, such as LASSO or Bayes C. However, CNN superiority is less significant when the prediction target is CGM parameters rather than the integrative trait ( Figure 3 ). Table 3 shows that, across all CGM parameters, the NMAE is similar for all regression methods. The same trend in parameter prediction accuracy as highlighted before is true for the three regression methods. As mentioned previously, it is suggested that the crop models decompose integrative traits into potentially more heritable intermediate traits (organ size, number, etc.) that are simulated through response curves to the environment, defined by genotypic parameters and can thus be used to predict the genetic variability of these integrative traits (Reymond et al., 2003; Parent and Tardieu, 2014). Indeed, the heritability of the integrative trait Biomaerofw was of 0.714 while the elementary traits composing aboveground biomass systematically harbored higher heritability with the exception of the number of tillers (see Supplementary Table S2 ). Most of these traits were indeed simulated with higher accuracy than the direct prediction of aboveground biomass (see Supplementary Figure S1 ).

It is interesting to note that the prediction of the trait through a CGM is also more accurate for the LASSO and Bayes C methods compared to the direct prediction of the trait. These results may show that genotypic parameters are under less complex genetic control than the integrative traits of interest. It could also indicate the robustness of the crop model for parameter sets that deviate slightly from the observed values. Lastly, considering causal effects of the environment on intermediate traits, as formalized in crop models, could enhance prediction accuracies of the integrative trait.