Multi-omic integration using interpretable machine learning reveals genetic mechanisms of trait variation and phenotypic plasticity in switchgrass

  1. Paulo Izquierdo
  2. Xiaoyu Weng
  3. Thomas E Juenger
  4. Jason Bonnette
  5. Yuko Yoshinaga
  6. Chris Daum
  7. Anna Lipzen
  8. Kerrie Barry
  9. Matthew Blow
  10. Melissa D Lehti-Shiu
  11. David B Lowry
  12. Shin-Han Shiu

  • Curated by eLife

    eLife logo

    eLife Assessment

    This study provides invaluable insights into the field of genomic prediction of traits by demonstrating how machine learning algorithms can integrate genomic and transcriptomic data and handle nonlinearity, providing important comparative prediction results across environments. The rationale and analyses conducted to integrate the two types of ~omics datasets in two environments are solid. The work will be of interest to plant geneticists and breeders, especially those interested in the genetics of environmental adaptation.

Reviewed by

Read the full article

Listed in

Log in to save this article

Abstract

Uncovering the genetic architecture of quantitative traits is challenging because polygenic control yields small individual gene effects and because gene–gene and genotype-by-environment interactions add further complexity. To understand the genetic basis of polygenic traits and their plasticity across environments, we integrated genome-wide SNPs and RNA-seq transcript data with interpretable statistical and machine learning models in a switchgrass (Panicum virgatum) diversity panel grown at contrasting field sites in Michigan and Texas. Notably, in addition to single environments, our trait prediction models were able to predict phenotypic differences, across environments i.e., plasticity. By interpreting trait prediction models with explainable machine learning methods, we identified genes that are the most predictive of flowering time and annual biomass production across environments, based on their associated gene expression levels and nearby SNPs. This approach recovered canonical flowering regulators and revealed novel, environment-specific candidate flowering genes. Further, transcriptome models consistently recovered more switchgrass genes homologous to experimentally validated genes in Arabidopsis and rice than SNP-based models. Feature interaction scores from the models also allow the identification of trait- and environment-dependent gene–gene interactions, where flowering time showed stronger and more abundant interactions than biomass. While some of the interactions identified are consistent with the link between flowering time and yield, most are novel predictors that need to be further evaluated. These findings provide insight into the genetic basis of trait variation and plasticity across environments and help prioritize candidate genes for functional validation and cultivar improvement.

Article activity feed

  1. eLife

    eLife Assessment

    This study provides invaluable insights into the field of genomic prediction of traits by demonstrating how machine learning algorithms can integrate genomic and transcriptomic data and handle nonlinearity, providing important comparative prediction results across environments. The rationale and analyses conducted to integrate the two types of ~omics datasets in two environments are solid. The work will be of interest to plant geneticists and breeders, especially those interested in the genetics of environmental adaptation.

  2. eLife

    Reviewer #2 (Public review):

    The authors examine the relative contributions of SNP and transcriptomic information to predicting trait variation and phenotypic plasticity in switchgrass. Using a diversity panel of 426 genotypes evaluated in Texas and Michigan, USA., this study demonstrates that transcriptomic predictors could outperform genomic markers in predicting environmental differences for traits like biomass yield and tiller count, highlighting the trait-dependent advantage of transcriptomic information over SNPs. However, combining SNPs and transcriptomic predictors generally yielded limited improvement in prediction accuracy.

    A principal limitation of this study concerns the experimental design and the robustness of transcriptomic predictors. With only one transcriptomic sample per genotype in each environment, within-genotype variation in transcript abundance cannot be estimated, making it difficult to distinguish genotype-specific expression responses from individual-plant variation. Unlike SNP genotypic information, which provides a stable representation of a genotype's inherited genetic variation underlying phenotypic differences, transcript abundance is highly sensitive to environmental conditions, developmental stages, and physiological state at the time of sampling. As a result, the stability of the transcriptomic predictors and their predictive performance across independent environments, years, or sampling conditions remain untested. Furthermore, the stability of feature selection across different training populations has not been established, raising additional questions about the reproducibility of the identified candidate genes. Together, these limitations constrain their generalizability and the biological interpretation of transcriptomic associations with phenotypic plasticity.

    Overall, this study provides valuable comparative evidence regarding the predictive performance of genomic and transcriptomic information. Its principal contribution lies in demonstrating that the relative predictive advantages of these information sources are trait-dependent and their integration does not necessarily improve prediction accuracy. While the predictive findings are supported, the biological interpretations extend beyond the evidence presented, as the identified transcriptomic predictors and their interactions do not establish causal genetic effects, epistasis, or direct regulatory relationships. Therefore, the proposed genetic mechanisms and environment-specific interactions should be regarded as hypotheses requiring independent biological validation, rather than established mechanistic findings.

  3. eLife

    Author response:

    The following is the authors’ response to the original reviews.

    eLife Assessment:

    The study by Izquierdo and colleagues provides important insights into the field of genomic and transcriptomic prediction of traits across multiple environments. The rationale and analyses conducted to integrate the two types of ~omics datasets across two environments are solid. However, some clarification would be appreciated in the presentation of the results, and adding some statistical control to clarify how the predictors were selected, or assessing their importance using the SHAP framework, would further consolidate the findings.

    We appreciate the positive feedback and the constructive evaluation of our work. In response, we have added further clarification to the Results section and provided additional detail on the feature selection procedure to facilitate understanding of how predictors were identified within the modeling framework.

    Public Reviews:

    Reviewer #1 (Public review):

    Summary:

    P. Izquierdo et al. investigated the genetic determinism of various traits of interest in switchgrass using large-scale genomic and transcriptomic data. More specifically, they worked on a diversity panel comprising 426 genotypes evaluated in common-garden experiments at two locations (Michigan and Texas). The phenotypic and genomic data were already published. In this work, they produced transcriptomic data for each of the 426 genotypes at each site, and they carried out phenotype predictions using genomic and transcriptomic data separately or together. While they were moderately correlated at each location, both omic information appeared to be complementary for the prediction of phenotype. To further exploit the fact that they have data across two locations, they computed differences for phenotypes and transcripts between locations as indicators of trait and transcript plasticity, respectively. They built predictive models of trait plasticity using genomic information and transcript plasticity, which proved to be quite accurate for traits affected by GxE. Finally, they made use of SHAP values from predictive models of flowering time and biomass at each location, as well as for their plasticity, to gain insight into their genetic determinism. These SHAP values provide the importance of the predictive features (SNP and/or transcripts) for trait prediction. This allowed them to confirm some candidate genes and to propose new candidates for both traits.

    We appreciate the reviewer's thorough and constructive summary of our work. We would like to clarify that while the genomic data were previously published, the transcriptomic and phenotypic data were generated as part of the present study. We have revised the text to make this distinction explicit.

    Strengths:

    I found this study interesting and rich. I think the sample size (426 genotypes) is large enough to support the findings. The use of a modern machine-learning approach (XGBoost) together with SHAP indices to find interesting features and get insights into the biological mechanisms underlying flowering time and biomass production is quite original. The methodology employed is globally sound. I also like the fact that the authors accounted implicitly for the population structure by providing a baseline prediction using the first 5 PCs.

    Weaknesses:

    While the methodology is globally sound, I sometimes had difficulties following exactly what was done. This is partly due to the fact that the authors used 2 omics (SNPs and transcripts) to predict phenotypes, and sometimes, in the results, it is not clear which of the 2 is the focus. This was especially the case for the importance of the features and the interpretability of the models, where I found it sometimes hard to tell whether the analysis was done on SNPs or transcripts.

    We have revised the text throughout to clarify when models using SNP, transcript (T), or both data types (S+T) were the focus of the analysis.

    Also, regarding the methodology, I did not understand why the authors needed to perform a feature selection approach. Maybe it was required to perform the interaction analysis, which could not be deployed on all the features? But regarding the importance of the features, I do not get the added value of the selection over the direct use of SHAP indices when using all features. Maybe this is because I am not a specialist in this kind of approach, but maybe the authors could add more details to explain the rationale behind the feature selection.

    Feature selection was performed to identify the subset of features that optimized predictive accuracy while reducing model complexity. By retaining only the most informative features, the resulting models achieve higher predictive accuracy than models trained on all features, and model interpretation is simplified. The rationale and procedure are described.

    Regarding SHAP interactions specifically, the number of features used to estimate SHAP interaction scores was limited to 5,000 due to the high computational demands of this analysis. Interaction scores are estimated for each individual in the training population, producing a 5,000 × 5,000 interaction matrix per individual. With 306 individuals in the training population and 10 independent model runs, this corresponds to 3,060 matrices of this size.

    Reviewer #2 (Public review):

    Summary:

    The authors aimed to evaluate whether integrating genomic (SNP) and transcriptomic information with machine learning can improve phenotypic prediction of polygenic traits across environments. The manuscript explored not only the predictability across models and predictor feature sets, but also attempted to identify meaningful genes and interactions underlying trait variation.

    Strengths:

    The main strength of the manuscript is its integration of SNP, transcriptomic, and phenotype datasets for 426 sorghum genotypes between Texas and Michigan. It provides a systematic comparison of predictor types (SNP versus transcriptomic abundance) and model strategies to integrate them.

    Weaknesses:

    (1) Experimental Design

    The experimental design raises several concerns that should be clarified before strong biological conclusions are drawn from the transcriptomic analyses.

    First, the transcriptomic sampling is not well aligned with the developmental stages most relevant to the phenotypes being modeled. Leaf tissue was collected at a single time point in each environment, whereas traits such as flowering time, biomass, tiller count, and panicle height arise from developmental processes occurring over extended and potentially distinct temporal windows. Consequently, the measured expression profiles are likely to reflect physiological states specific to the sampling dates (May 5-6 in Texas and June 22-24 in Michigan) rather than the regulatory processes underlying the target phenotypes.

    We selected the pre-flowering stage because it represents a critical developmental transition during which plants remain largely vegetative while initiating reproductive development. Gene expression at this stage is therefore expected to capture regulatory processes associated with flowering initiation while also reflecting pathways involved in vegetative growth and biomass accumulation. In our previous work in maize (Azodi et al., 2020), we demonstrated that transcriptomic profiles obtained at a single developmental stage can contain information predictive of complex traits. Although the expression profiles analyzed here do not capture the full temporal dynamics underlying traits such as flowering time and biomass, they provide a biologically meaningful snapshot of molecular variation relevant to these traits. We have revised the manuscript to acknowledge this limitation and to clarify that our findings should be interpreted as associations between gene expression measured at the pre-flowering stage and the phenotypes of interest, rather than as evidence that all identified genes directly regulate these traits throughout development.

    Second, the phrase "haphazardly randomized" is questionable for a field experiment. It is unclear whether the design included formal randomization, blocking, row/column structure, or spatial correction. Without explicit accounting for spatial field heterogeneity, environmental variation within sites may confound genotype and transcriptomic effects.

    We thank the reviewer for this comment and agree that the original wording was unclear. The text has been revised to more precisely describe the experimental design. The establishment and management of the common gardens are described in detail in Lovell et al. (2021). Briefly, the diversity panel was established as a single-block planting at each location, with genotypes distributed throughout the field and arranged in a honeycomb design in which each individual had four nearest neighbors located 1.5 m apart. We have removed the phrase "haphazardly randomized," which was imprecise.

    Because only a single replicate of each genotype was established at each site, formal spatial correction was not performed, as genotype and within-site spatial effects cannot be estimated independently in this design. We acknowledge that additional replication would have improved our ability to account for fine-scale environmental heterogeneity.

    However, several aspects of the experimental design reduce the likelihood that spatial variation substantially confounded the results. The diversity panel is part of a coordinated common garden network spanning ten field sites across the United States, designed primarily to evaluate genotype-by-environment interactions across contrasting environments rather than within-site spatial variation. In addition, the AP13 reference genotype was distributed throughout the experimental plots and did not show evidence of substantial spatial heterogeneity at the Michigan site, suggesting that within-site environmental variation was limited. We have added text acknowledging this limitation and its implications for interpreting genotype and genotype-by-environment effects.

    Third, the Methods do not clearly describe biological replication for RNA-seq. If each genotype-by-environment combination were represented by a single transcriptomic sample, then within-genotype expression variance cannot be estimated. This is important because transcript abundance is highly sensitive to microenvironment, sampling time, tissue status, developmental stage, and technical variation. The absence of replication significantly weakens confidence in gene-level feature importance and gene-gene interaction claims.

    We appreciate this thoughtful comment and fully acknowledge that the absence of biological replication for RNA-seq within each genotype-by-environment combination is an important limitation of this study. Because the diversity panel consisted of a single plant per genotype at each field site, transcriptomic data were generated from one biological sample per genotype-by-environment combination, and within-genotype expression variance could not be directly estimated. We recognize that transcript abundance can be influenced by microenvironmental conditions, developmental stage, tissue status, sampling time, and technical variation, all of which may affect confidence in gene-level feature importance and inferred interaction effects.

    The scale of the experiment, involving 426 genotypes evaluated across contrasting common gardens, made biological replication for RNA-seq financially and logistically impractical. This experimental design was established to maximize representation of genetic diversity and environmental variation across a large multi-environment common garden network rather than to obtain replicated transcriptomic measurements within individual sites. To minimize uncontrolled variation, tissue was collected at both locations at approximately the same pre-flowering developmental stage, and standardized protocols for tissue collection, library preparation, sequencing, and data processing were applied across all samples.

    Although individual genotype-by-environment combinations were not biologically replicated, the analyses leverage substantial replication at the level of alleles and population structure, as individual alleles and transcripts are represented across many genetically distinct accessions spanning multiple ecotypes. This is analogous to many human population-scale transcriptomic studies, where each genotype is typically represented by a single biological sample and inference is made across genetically diverse individuals rather than replicated measurements of the same genotype.

    Importantly, the goal of this study was to identify transcriptomic patterns that are predictive of phenotypic variation across genotypes and environments, rather than to establish causal regulatory relationships. Accordingly, the genes and interaction patterns identified through model interpretation should be viewed as candidate predictive associations and hypothesis-generating signals that warrant validation in future studies with replicated sampling and controlled experiments. We have revised the manuscript to more clearly acknowledge these limitations and to align the interpretation of gene importance and interaction analyses with the predictive scope of the study.

    Four, the analysis of expression differences across environments is based on a simple subtraction (TX - MI) followed by correlation with genetic similarity. This approach is not standard in transcriptomic analysis and does not account for variability, replication, or statistical uncertainty. Conventional methods for assessing differential expression and genotype-by-environment interactions rely on model-based frameworks that explicitly estimate variance components and test for interaction effects. Without such modeling, the observed expression differences may reflect noise or confounding factors rather than genotype-driven responses.

    We appreciate this comment. Our approach to estimating expression differences across environments was intended to quantify transcriptomic plasticity and evaluate whether variation in gene expression responses across environments could help explain and predict phenotypic plasticity. While this is not a conventional differential expression framework, our objective was not to identify differentially expressed genes or formally test genotype-by-environment interaction effects at the transcript level.

    Expression differences were calculated as TX − MI using log-transformed expression values. On the logarithmic scale, this difference is equivalent to the log expression ratio between environments and therefore provides a straightforward measure of the magnitude and direction of transcriptomic responses across locations. We used this metric as a descriptive measure of transcriptomic plasticity that could be related to genetic similarity and incorporated into predictive models of phenotypic plasticity. Consistent with this rationale, models based on expression differences across environments were able to predict differences in phenotypic performance between locations.

    (2) SHAP contribution values

    Although SHAP is a well-established framework for decomposing model predictions into feature-level contributions, its use in this manuscript raises several concerns regarding interpretation, statistical validity, and biological inference.

    First, SHAP values quantify the contribution of features within the fitted model, conditional on the joint distribution of inputs and the model structure. They do not represent causal effects or direct biological importance.

    We agree that SHAP values do not represent causal effects or direct measures of biological importance. Rather, SHAP quantifies the contribution of features to model predictions within the context of a fitted model and the observed input distribution. Our rationale for using SHAP was not to infer causality, but to identify features that are consistently informative for predicting phenotypic variation. The underlying hypothesis is that features contributing strongly to accurate predictions may be associated with biological processes influencing the trait and can therefore serve as candidate genes or genomic regions for further investigation. This strategy has been widely adopted in genomics and plant biology, where interpretable machine learning methods are used to extract biological insights from predictive models. Accordingly, the genes and interactions highlighted by SHAP should be interpreted as predictive associations and candidate mechanisms requiring validation through independent experiments and functional studies. We have revised the manuscript to more clearly distinguish predictive importance from causal biological effects.

    There is a difference where SHAP values are often in log-odds and the regression model uses absolute units.

    We thank the reviewer for this comment. We believe there may be a misunderstanding regarding the scale of the SHAP values reported in this study. The models analyzed here are regression models predicting continuous traits, not classification models. Consequently, SHAP values are expressed in the same units as the predicted phenotype and represent the additive contribution of each feature to individual predictions. The interpretation of SHAP values in log-odds units is specific to certain classification models where predictions are generated on the log-odds scale. We have revised the manuscript to clarify this point and avoid potential confusion.

    Without a fair evaluation of model fit, the interpretation of SHAP values needs to take a cautious step because a model could fit poorly when a feature shows very high SHAP values.

    We agree that SHAP values should be interpreted in the context of model performance. For this reason, model performance was evaluated prior to feature interpretation. Figure 1 presents the predictive performance of all models, and we focused our biological interpretation on flowering time and biomass because these traits are agronomically important and showed relatively stronger predictive performance compared to other traits, providing a more reliable basis for interpreting feature importance. We have revised the manuscript to emphasize that model interpretation should be considered in conjunction with model accuracy.

    In genomic data, where features are highly correlated due to linkage disequilibrium and co-expression, SHAP values can distribute contribution values across correlated variables in ways that are not uniquely identifiable. As a result, features highlighted as "important" may reflect correlation structure rather than true functional relevance.

    We agree that multicollinearity can complicate the interpretation of feature importance, as correlated features may share predictive information, resulting in importance being distributed among them or assigned preferentially to one feature over another. This issue is particularly relevant in genomic datasets due to linkage disequilibrium among SNPs and co-expression among genes. To reduce redundancy among genetic markers, SNPs were filtered using a linkage disequilibrium threshold of r2 < 0.1 prior to model training, substantially reducing the number of highly correlated markers. We acknowledge, however, that correlation among features, particularly among co-expressed genes, remains present and may influence SHAP-based feature attribution.

    To improve the robustness of feature interpretation, models were trained across 10 independent runs and SHAP values were averaged across runs, reducing sensitivity to individual model fits and providing more stable estimates of feature contributions. Nevertheless, SHAP values should not be interpreted as uniquely identifying the most functionally relevant feature within a correlated group of predictors. The genes and genetic variants highlighted by SHAP should therefore be viewed as candidate predictive features associated with trait variation rather than definitive evidence of causal biological relevance.

    This correlative structure can be exacerbated in this manuscript because of the use of TPM-normalized transcript abundances as predictor variables without biological replicates. Assume the estimates of transcript abundances are robust, TPM values are compositional, with a constant-sum constraint that creates dependencies among all genes that induce negative correlations. This issue is particularly relevant for the interpretation of gene importance and interaction effects, where correlated predictors can lead to unstable and non-unique attributions. This biological interpretation of transcript-based features remains uncertain.

    We agree that TPM-normalized expression values are compositional and may introduce dependencies among genes due to the constant-sum constraint. TPM normalization was selected because it facilitates comparisons of transcript abundance across samples by accounting for differences in sequencing depth and transcript length. Expression values were further transformed using log10(TPM + 1), which reduces skewness and stabilizes variance, making the data more suitable for machine learning analyses. We acknowledge that compositional effects and gene co-expression may influence feature attribution and the interpretation of gene importance. However, the primary objective of this study was predictive rather than inferential. In this context, log-transformed TPM values provide a standardized and comparable representation of transcript abundance across genotypes and environments, serving as consistent input for predictive models. Consequently, genes identified through SHAP analyses should be interpreted as candidate predictive features associated with trait variation, rather than as definitive evidence of direct causal biological importance.

    (3) Result interpretation

    For example, in page 11, "plasticity SNP- and transcriptomic-based models generally outperformed single-environment models for traits with low cross-environment correlation, such as green-up (Fig. 2c, r = -0.13, p < 8.3 × 10-3) and tiller count (Fig. 2f, r = -0.08, p = 0.1) (Supplementary Fig. S1)." is too broad. For green-up, the Diff model appears much better than MI, but not clearly better than TX.

    We thank the reviewer for this observation and agree that the original wording was too broad. We have revised the text to more accurately reflect the observed pattern, noting that difference-based models showed relatively stronger predictive performance for traits with low cross-environment correlation, while generally performing worse than single-environment models for traits with high cross-environment correlation. We believe this pattern is biologically meaningful, as transcriptomic differences across environments are expected to be most informative when phenotypic responses differ substantially between locations. We have revised to more accurately characterize these results.

    And, same page 11, "...Diffexp was more predictive than SNPs for trait plasticity in biomass, flowering time, and tiller count..." only holds true for biomass, not flowering time, or tiller count.

    We thank the reviewer for this observation and agree that the original statement was too strong. Diffexp models showed higher predictive performance for biomass and tiller count, whereas performance for flowering time was comparable between Diffexp and SNP models. We have revised to more accurately characterize these results.

    The aspect of "complementary information" between SNP and transcriptomic models in page 12 is stronger than what is supported by Figure 2. Figure 2 shows different predictive performance, but it does not by itself demonstrate complementarity. Establishing complementarity requires evidence that combining SNP+T improves prediction consistently or captures distinct, non-overlapping signals. Yet the preceding section says SNP+T outperformed either single data type in only 15% of cases, with modest gains. This is confusing. Also, there was not G+T in Figure 2; it is SNP+T.

    We agree that the term "complementary" was stronger than what the results in Figure 2 directly support. We have revised the text to describe these findings more cautiously, emphasizing that SNP and transcriptomic models may capture partially distinct aspects of biological variation. In the following section, we further address this point by showing that the two model types identify largely distinct sets of predictive features, suggesting that they reflect different sources of information relevant to trait prediction. We also thank the reviewer for noting the SNP+T vs. G+T inconsistency, which has been corrected throughout the manuscript.

    Recommendations for the authors:

    Reviewer #1 (Recommendations for the authors):

    (1) I found that the title does not properly reflect the content of the paper. The term "artificial intelligence" is a bit misleading: I was expecting deep learning-based predictions. GxE and plasticity are clear strengths of the paper that are not present at all in the title. The same holds for multi-omic integration. Maybe you could consider revisiting the title.

    We thank the reviewer for this suggestion and have revised the title to: "Multiomic integration using interpretable machine learning reveals genetic mechanisms of trait variation and phenotypic plasticity in switchgrass." This revision highlights the integration of genomic and transcriptomic data, explicitly mentions phenotypic plasticity, and replaces "artificial intelligence" with "interpretable machine learning," which more accurately describes the modeling framework employed.

    (2) In general, the quality of the figures is very poor, and it did not allow me to perform a deep assessment of their content, which is sometimes hardly visible.

    We agree that figure quality in the submitted PDF was suboptimal. The reduction in image resolution occurred during conversion from the original Word document to PDF. We have uploaded a revised manuscript with higher-resolution figures to ensure all elements are clearly visible.

    (3) To estimate plasticity in phenotypes and transcripts, you use differences between locations. Other options have been proposed in the literature, like ratios, for example. I think the choice of the difference is probably a good one, but maybe you could justify it in the introduction.

    We chose to estimate phenotypic and transcriptomic plasticity using differences between locations because this metric preserves both the magnitude and direction of environmental responses while avoiding the instability associated with ratio-based measures, particularly when values are small. Furthermore, because transcript abundance values were modeled using log10(TPM + 1)-transformed data, differences between environments are conceptually equivalent to log expression ratios, as log(A) − log(B) = log(A/B). This justification has been added to the manuscript.

    (4) Figure 1. I wonder why you computed correlation matrices between genotypes for transcriptomic data instead of simply covariance matrices, which would be closer to the genomic relationship matrix.

    Our intention in Figure 1 was to visualize overall similarity among samples based on transcriptomic expression profiles across genes. We used Pearson correlation because it is widely applied in transcriptomic analyses to capture relative expression patterns while remaining independent of absolute expression magnitude and scale differences among samples. Covariance matrices, by contrast, are strongly influenced by feature variance and scale, which can complicate interpretation of overall transcriptomic similarity in high-dimensional expression datasets.

    (5) Figure 1. You calculate correlations between phenotypic distance and genomic distances for the different traits, but I think this is somehow related to the genomic heritability estimation. Could you please elaborate on this or provide the rationale under your analysis?

    We agree that the correlation between genomic relationships and phenotypic similarity is conceptually related to genomic heritability, as both analyses evaluate the extent to which genetic similarity is associated with phenotypic variation. However, their objectives differ. The correlation analysis between SNP, eCor, and phenotypic distance was intended as an exploratory assessment of whether genetically or transcriptomically similar genotypes tend to exhibit similar phenotypic responses across traits and environments. Genomic heritability estimation, in contrast, provides a formal variance component framework to quantify the proportion of phenotypic variance explained by genome-wide marker information. While related conceptually, these analyses address complementary but distinct aspects of the relationship between genomic, transcriptomic, and phenotypic variation.

    (6) L163-164: "Because genotypes with similar genetic information are expected to exhibit similar trait value": this is the definition of the heritability, so I would say that this is true to the extent of the heritability of the trait... Maybe you could consider modifying this sentence.

    Thank you for this suggestion. Our intention was not to invoke heritability directly, but rather to assess the associations among kinship, eCor, and phenotypic similarity within each environment and to compare the patterns captured by these different data types. We agree that the original wording could be interpreted as implying a deterministic relationship between genetic similarity and phenotypic resemblance. We have revised the text to focus on the association between genetic relatedness and similarity in expression and phenotypic profiles.

    (7) Figure 2. For the prediction with transcriptomic data, I understood that you matched data by site for the prediction, which is OK. But I wonder what would happen if you had used transcriptomic data from the other site and vice versa for the prediction, how much predictability you would have lost? I guess that this would somehow be related to the extent of GxE for each trait. Maybe it would be something interesting to look at?

    We thank the reviewer for this thoughtful suggestion. While cross-environment prediction is an interesting avenue for future research, our primary objective was not to maximize prediction across environments, as predictive ability would be expected to decrease substantially under strong G×E effects. Instead, we focused on single-environment models to: (1) identify transcriptomic features associated with trait prediction within each environment, and (2) evaluate which features were consistently important across environments. Furthermore, because G×E effects were observed, we developed TX − MI difference models to identify transcriptomic features associated with phenotypic plasticity. We therefore believe that the within-environment and plasticity-focused analyses were more closely aligned with the objectives of this study than cross-environment prediction would have been.

    (8) The calculation of local SHAP values for each genotype mentioned L331 and the clustering of genotypes based on such values, as presented in Figure 4, are not presented in the method section. Could you please provide more details?

    The estimation of SHAP values is described in the Methods section. Briefly, local SHAP values were calculated for each genotype, and genotypes were subsequently clustered based on these values using hierarchical clustering with Euclidean distance.

    (9) Gene-by-gene interactions paragraph: I am a bit confused here regarding which omic was used for looking at the interaction. Naively, I would have gone with the SNPs, but maybe you would rather use the transcripts instead? This would explain why in Figure 5 c-h, the x axes seemed to correspond to transcripts rather than SNPs. Could you please clarify and explain why you focused here on transcripts rather than SNPs?

    We thank the reviewer for this comment. To clarify, gene–gene interactions were assessed using transcriptomic models. This choice was motivated by three observations: (1) transcriptomic features were enriched among the most important predictors in the multi-omic models; (2) transcriptomic models identified a greater number of benchmark genes than SNP-based models; and (3) gene expression data provide a more biologically interpretable framework for assessing gene–gene interactions, as transcripts more directly reflect gene activity and regulatory processes than individual genetic variants. We have revised the manuscript to make this rationale explicit.

    (10) Materials and Methods:

    - Were there replicates of the genotypes within each site?

    This diversity panel is part of a multi-institutional project maintaining 10 common gardens across the United States. Given the scale of the experiment and the large number of genotypes evaluated, the design included a single plant per genotype at each site rather than within-site replication. We acknowledge this as a limitation and have discussed it explicitly in the revised manuscript.

    - Have the data (phenotypic/transcriptomic) been corrected for spatial effects?

    Because only a single replicate per genotype was evaluated at each site, genotype and local spatial effects could not be estimated independently. Applying a spatial correction under these conditions could inadvertently remove true biological variation. The randomized honeycomb design was intended to minimize systematic spatial confounding; however, we acknowledge that residual environmental heterogeneity may still contribute to phenotypic and transcriptomic variation.

    - L502-503: Could you please write the model formally instead of R code?

    The model has been written formally.

    - Genomic heritability estimates: I wonder how useful they really are? I suppose you could have used the model L502-503 to estimate genomic heritability in general across the 2 sites.

    We agree that the mixed model described in the revised manuscript could be used to estimate genomic heritability jointly across both sites. In this study, however, genomic heritability was estimated within individual environments to evaluate trait predictability within each location and to assess the heritable component of phenotypic plasticity. Genomic heritability also provides an important benchmark for interpreting prediction results, as the theoretical upper limit of prediction accuracy is approximately proportional to the square root of heritability under quantitative genetic theory.

    - L523: I am not sure what LD> 0.1 means. Could you please give more details, e.g., which measure of LD, which software, and which parameters?

    We have clarified this point in the Methods. LD pruning was performed using the +prune plugin in BCFtools, with pairwise LD measured as the squared Pearson correlation coefficient (r2). SNPs with r2 > 0.1 within a 5 Mb sliding window were removed to retain a set of more independent markers.

    - L526: How many leaves were sampled per plant? Or maybe a particular leaf was chosen?

    For each plant, we collected the fully expanded leaves at the highest point on each tiller. The proximal ~4 cm (length of microcentrifuge tube) of that harvested leaf was then flash-frozen in liquid nitrogen for subsequent RNA extractions.

    - L592-597: I understand that training and test sets were stratified according to biomass (log-transformed) in TX. Could you please explain why this trait/location was chosen for the stratification?

    Training/test stratification: We chose log-transformed biomass in TX for stratification because biomass is the primary agronomic trait of interest in switchgrass and exhibits greater phenotypic variation in TX compared to MI, allowing stratification to better preserve the full range of values in both the training and testing sets. The same training/testing split was subsequently applied to all analyses to ensure consistent and fair comparisons across traits, environments, and data types.

    (11) Data availability: Maybe the read counts could also be provided on top of the raw RNA sequences. This would ensure, to some extent, reproducibility of the analyses, especially because the code for processing RNA-seq reads does not seem to be released in the GitHub repository.

    The RNA-seq pipeline used in this study is described in the revised manuscript. The pipeline code is maintained in an internal GitLab repository at JGI and is therefore only accessible within JGI systems. To facilitate reproducibility, we have deposited the TPM count matrix used in our analyses in the following Zenodo repository: https://zenodo.org/records/21781035

    (12) Code availability: The code to run Bayesian regression models does not seem to be released.

    The code for the Bayesian regression models has been added to the manuscript repository.

    Reviewer #2 (Recommendations for the authors):

    The manuscript will benefit the most from additional biological replicates for RNAseq, as well as the phenology alignment with data collection. Further recommendations for improving the clarity in experimental design and alignment between claims and evidence.

    We agree that additional biological replication for RNA-seq and transcriptomic sampling across multiple developmental stages would further strengthen the study. However, given the scale of the experiment, 426 genotypes evaluated across contrasting field environments, generating replicated transcriptomic datasets at multiple developmental stages was not logistically or financially feasible. Despite these constraints, the experimental design enabled characterization of transcriptomic and phenotypic variation across a large natural diversity panel and the identification of predictive patterns associated with trait variation across environments. We have revised the manuscript to more clearly describe the experimental design, sampling strategy, and study limitations, and to ensure that conclusions are appropriately aligned with the evidence presented.

    Minor corrections:

    (1) MAF threshold in Line 522, page 27. It should be > 0.05?

    Corrected to > 0.05.

    (2) SNP+T versus G+T.

    Revised to SNP+T consistently throughout the manuscript.

    (3) Include FDR-adjusted p-value thresholds.

    FDR-adjusted p-value thresholds: For the Pearson correlations reported in this study, each represents a single pre-specified comparison rather than a family of simultaneous tests, and raw p-values are therefore appropriate. For the benchmark gene enrichment and Wilcoxon rank-sum tests, where multiple comparisons were performed across models and traits, FDR correction has been applied using the Benjamini– Hochberg procedure and the reported p-values have been updated accordingly.

    (4) Avoid ambiguous phrasing like "generally consistent" and "generally outperformed".

    Revised throughout the manuscript.

    (5) The mention of artificial intelligence in the Abstract is largely conceptual and reads more like discussion material rather than the main finding of the manuscript, per eLife's Abstract requirements.

    The mention of artificial intelligence has been revised in the Abstract in accordance with eLife's guidelines.

  4. Version published to 10.7554/elife.111208.2 on eLife
  5. Version published to 10.7554/elife.111208 on eLife
  6. Version published to 10.64898/2026.03.06.710154v3 on bioRxiv
  7. Version published to 10.64898/2026.03.06.710154v2 on bioRxiv
  8. eLife

    eLife Assessment

    The study by Izquierdo and colleagues provides important insights into the field of genomic and transcriptomic prediction of traits across multiple environments. The rationale and analyses conducted to integrate the two types of ~omics datasets across two environments are solid. However, some clarification would be appreciated in the presentation of the results, and adding some statistical control to clarify how the predictors were selected, or assessing their importance using the SHAP framework, would further consolidate the findings.

  9. eLife

    Reviewer #1 (Public review):

    Summary:

    P. Izquierdo et al. investigated the genetic determinism of various traits of interest in switchgrass using large-scale genomic and transcriptomic data. More specifically, they worked on a diversity panel comprising 426 genotypes evaluated in common-garden experiments at two locations (Michigan and Texas). The phenotypic and genomic data were already published. In this work, they produced transcriptomic data for each of the 426 genotypes at each site, and they carried out phenotype predictions using genomic and transcriptomic data separately or together. While they were moderately correlated at each location, both omic information appeared to be complementary for the prediction of phenotype. To further exploit the fact that they have data across two locations, they computed differences for phenotypes and transcripts between locations as indicators of trait and transcript plasticity, respectively. They built predictive models of trait plasticity using genomic information and transcript plasticity, which proved to be quite accurate for traits affected by GxE. Finally, they made use of SHAP values from predictive models of flowering time and biomass at each location, as well as for their plasticity, to gain insight into their genetic determinism. These SHAP values provide the importance of the predictive features (SNP and/or transcripts) for trait prediction. This allowed them to confirm some candidate genes and to propose new candidates for both traits.

    Strengths:

    I found this study interesting and rich. I think the sample size (426 genotypes) is large enough to support the findings. The use of a modern machine-learning approach (XGBoost) together with SHAP indices to find interesting features and get insights into the biological mechanisms underlying flowering time and biomass production is quite original. The methodology employed is globally sound. I also like the fact that the authors accounted implicitly for the population structure by providing a baseline prediction using the first 5 PCs.

    Weaknesses:

    While the methodology is globally sound, I sometimes had difficulties following exactly what was done. This is partly due to the fact that the authors used 2 omics (SNPs and transcripts) to predict phenotypes, and sometimes, in the results, it is not clear which of the 2 is the focus. This was especially the case for the importance of the features and the interpretability of the models, where I found it sometimes hard to tell whether the analysis was done on SNPs or transcripts.

    Also, regarding the methodology, I did not understand why the authors needed to perform a feature selection approach. Maybe it was required to perform the interaction analysis, which could not be deployed on all the features? But regarding the importance of the features, I do not get the added value of the selection over the direct use of SHAP indices when using all features. Maybe this is because I am not a specialist in this kind of approach, but maybe the authors could add more details to explain the rationale behind the feature selection.

  10. eLife

    Reviewer #2 (Public review):

    Summary:

    The authors aimed to evaluate whether integrating genomic (SNP) and transcriptomic information with machine learning can improve phenotypic prediction of polygenic traits across environments. The manuscript explored not only the predictability across models and predictor feature sets, but also attempted to identify meaningful genes and interactions underlying trait variation.

    Strengths:

    The main strength of the manuscript is its integration of SNP, transcriptomic, and phenotype datasets for 426 sorghum genotypes between Texas and Michigan. It provides a systematic comparison of predictor types (SNP versus transcriptomic abundance) and model strategies to integrate them.

    Weaknesses:

    (1) Experimental Design

    The experimental design raises several concerns that should be clarified before strong biological conclusions are drawn from the transcriptomic analyses.

    First, the transcriptomic sampling is not well aligned with the developmental stages most relevant to the phenotypes being modeled. Leaf tissue was collected at a single time point in each environment, whereas traits such as flowering time, biomass, tiller count, and panicle height arise from developmental processes occurring over extended and potentially distinct temporal windows. Consequently, the measured expression profiles are likely to reflect physiological states specific to the sampling dates (May 5-6 in Texas and June 22-24 in Michigan) rather than the regulatory processes underlying the target phenotypes.

    Second, the phrase "haphazardly randomized" is questionable for a field experiment. It is unclear whether the design included formal randomization, blocking, row/column structure, or spatial correction. Without explicit accounting for spatial field heterogeneity, environmental variation within sites may confound genotype and transcriptomic effects.

    Third, the Methods do not clearly describe biological replication for RNA-seq. If each genotype-by-environment combination were represented by a single transcriptomic sample, then within-genotype expression variance cannot be estimated. This is important because transcript abundance is highly sensitive to microenvironment, sampling time, tissue status, developmental stage, and technical variation. The absence of replication significantly weakens confidence in gene-level feature importance and gene-gene interaction claims.

    Four, the analysis of expression differences across environments is based on a simple subtraction (TX - MI) followed by correlation with genetic similarity. This approach is not standard in transcriptomic analysis and does not account for variability, replication, or statistical uncertainty. Conventional methods for assessing differential expression and genotype-by-environment interactions rely on model-based frameworks that explicitly estimate variance components and test for interaction effects. Without such modeling, the observed expression differences may reflect noise or confounding factors rather than genotype-driven responses.

    (2) SHAP contribution values

    Although SHAP is a well-established framework for decomposing model predictions into feature-level contributions, its use in this manuscript raises several concerns regarding interpretation, statistical validity, and biological inference.

    First, SHAP values quantify the contribution of features within the fitted model, conditional on the joint distribution of inputs and the model structure. They do not represent causal effects or direct biological importance. There is a difference where SHAP values are often in log-odds and the regression model uses absolute units. Without a fair evaluation of model fit, the interpretation of SHAP values needs to take a cautious step because a model could fit poorly when a feature shows very high SHAP values.

    In genomic data, where features are highly correlated due to linkage disequilibrium and co-expression, SHAP values can distribute contribution values across correlated variables in ways that are not uniquely identifiable. As a result, features highlighted as "important" may reflect correlation structure rather than true functional relevance.

    This correlative structure can be exacerbated in this manuscript because of the use of TPM-normalized transcript abundances as predictor variables without biological replicates. Assume the estimates of transcript abundances are robust, TPM values are compositional, with a constant-sum constraint that creates dependencies among all genes that induce negative correlations. This issue is particularly relevant for the interpretation of gene importance and interaction effects, where correlated predictors can lead to unstable and non-unique attributions. This biological interpretation of transcript-based features remains uncertain.

    (3) Result interpretation

    For example, in page 11, "plasticity SNP- and transcriptomic-based models generally outperformed single-environment models for traits with low cross-environment correlation, such as green-up (Fig. 2c, r = -0.13, p < 8.3 × 10⁻³) and tiller count (Fig. 2f, r = -0.08, p = 0.1) (Supplementary Fig. S1).", is too broad. For green-up, the Diff model appears much better than MI, but not clearly better than TX.

    And, same page 11, "...Diffexp was more predictive than SNPs for trait plasticity in biomass, flowering time, and tiller count..." only holds true for biomass, not flowering time, or tiller count.

    The aspect of "complementary information" between SNP and transcriptomic models in page 12 is stronger than what is supported by Figure 2. Figure 2 shows different predictive performance, but it does not by itself demonstrate complementarity. Establishing complementarity requires evidence that combining SNP+T improves prediction consistently or captures distinct, non-overlapping signals. Yet the preceding section says SNP+T outperformed either single data type in only 15% of cases, with modest gains. This is confusing. Also, there was not G+T in Figure 2; it is SNP+T.

  11. Version published to 10.7554/elife.111208.1 on eLife
  12. Mentioned by @geneticsmlpapers.bsky.social on Bluesky

    Uncovering genetic mechanisms underlying trait variation in switchgrass using explainable artificial intelligence https://www.biorxiv.org/content/10.64898/2026.03.06.710154v1

  13. Mentioned by @biorxiv-genetic.bsky.social on Bluesky

    Uncovering genetic mechanisms underlying trait variation in switchgrass using explainable artificial intelligence https://www.biorxiv.org/content/10.64898/2026.03.06.710154v1

  14. Mentioned by @biorxivpreprint.bsky.social on Bluesky

    Uncovering genetic mechanisms underlying trait variation in switchgrass using explainable artificial intelligence https://www.biorxiv.org/content/10.64898/2026.03.06.710154v1

  15. Version published to 10.64898/2026.03.06.710154v1 on bioRxiv