Reciprocal BLUP: A Predictability-Guided Multi-Omics Framework for Plant Phenotype Prediction

Reciprocal BLUP: A Predictability-Guided Multi-Omics Framework for Plant Phenotype Prediction | bioRxiv /* */ /* */ <!-- <!-- /*! * yepnope1.5.4 * (c) WTFPL, GPLv2 */ (function(a,b,c){function d(a){return"[object Function]"==o.call(a)}function e(a){return"string"==typeof a}function f(){}function g(a){return!a||"loaded"==a||"complete"==a||"uninitialized"==a}function h(){var a=p.shift();q=1,a?a.t?m(function(){("c"==a.t?B.injectCss:B.injectJs)(a.s,0,a.a,a.x,a.e,1)},0):(a(),h()):q=0}function i(a,c,d,e,f,i,j){function k(b){if(!o&&g(l.readyState)&&(u.r=o=1,!q&&h(),l.onload=l.onreadystatechange=null,b)){"img"!=a&&m(function(){t.removeChild(l)},50);for(var d in y[c])y[c].hasOwnProperty(d)&&y[c][d].onload()}}var j=j||B.errorTimeout,l=b.createElement(a),o=0,r=0,u={t:d,s:c,e:f,a:i,x:j};1===y[c]&&(r=1,y[c]=[]),"object"==a?l.data=c:(l.src=c,l.type=a),l.width=l.height="0",l.onerror=l.onload=l.onreadystatechange=function(){k.call(this,r)},p.splice(e,0,u),"img"!=a&&(r||2===y[c]?(t.insertBefore(l,s?null:n),m(k,j)):y[c].push(l))}function j(a,b,c,d,f){return q=0,b=b||"j",e(a)?i("c"==b?v:u,a,b,this.i++,c,d,f):(p.splice(this.i++,0,a),1==p.length&&h()),this}function k(){var a=B;return a.loader={load:j,i:0},a}var l=b.documentElement,m=a.setTimeout,n=b.getElementsByTagName("script")[0],o={}.toString,p=[],q=0,r="MozAppearance"in l.style,s=r&&!!b.createRange().compareNode,t=s?l:n.parentNode,l=a.opera&&"[object Opera]"==o.call(a.opera),l=!!b.attachEvent&&!l,u=r?"object":l?"script":"img",v=l?"script":u,w=Array.isArray||function(a){return"[object Array]"==o.call(a)},x=[],y={},z={timeout:function(a,b){return b.length&&(a.timeout=b[0]),a}},A,B;B=function(a){function b(a){var a=a.split("!"),b=x.length,c=a.pop(),d=a.length,c={url:c,origUrl:c,prefixes:a},e,f,g;for(f=0;f<d;f++)g=a[f].split("="),(e=z[g.shift()])&&(c=e(c,g));for(f=0;f<b;f++)c=x[f](c);return c}function g(a,e,f,g,h){var i=b(a),j=i.autoCallback;i.url.split(".").pop().split("?").shift(),i.bypass||(e&&(e=d(e)?e:e[a]||e[g]||e[a.split("/").pop().split("?")[0]]),i.instead?i.instead(a,e,f,g,h):(y[i.url]?i.noexec=!0:y[i.url]=1,f.load(i.url,i.forceCSS||!i.forceJS&&"css"==i.url.split(".").pop().split("?").shift()?"c":c,i.noexec,i.attrs,i.timeout),(d(e)||d(j))&&f.load(function(){k(),e&&e(i.origUrl,h,g),j&&j(i.origUrl,h,g),y[i.url]=2})))}function h(a,b){function c(a,c){if(a){if(e(a))c||(j=function(){var a=[].slice.call(arguments);k.apply(this,a),l()}),g(a,j,b,0,h);else if(Object(a)===a)for(n in m=function(){var b=0,c;for(c in a)a.hasOwnProperty(c)&&b++;return b}(),a)a.hasOwnProperty(n)&&(!c&&!--m&&(d(j)?j=function(){var a=[].slice.call(arguments);k.apply(this,a),l()}:j[n]=function(a){return function(){var b=[].slice.call(arguments);a&&a.apply(this,b),l()}}(k[n])),g(a[n],j,b,n,h))}else!c&&l()}var h=!!a.test,i=a.load||a.both,j=a.callback||f,k=j,l=a.complete||f,m,n;c(h?a.yep:a.nope,!!i),i&&c(i)}var i,j,l=this.yepnope.loader;if(e(a))g(a,0,l,0);else if(w(a))for(i=0;i (function(w,d,s,l,i){w[l]=w[l]||[];w[l].push({'gtm.start':new Date().getTime(),event:'gtm.js'});var f=d.getElementsByTagName(s)[0];var j=d.createElement(s);var dl=l!='dataLayer'?'&l='+l:'';j.src='//www.googletagmanager.com/gtm.js?id='+i+dl;j.type='text/javascript';j.async=true;f.parentNode.insertBefore(j,f);})(window,document,'script','dataLayer','GTM-M677548'); Skip to main content Home About Submit ALERTS / RSS Search for this keyword Advanced Search New Results Reciprocal BLUP: A Predictability-Guided Multi-Omics Framework for Plant Phenotype Prediction View ORCID Profile Hayato Yoshioka , View ORCID Profile Gota Morota , View ORCID Profile Hiroyoshi Iwata doi: https://doi.org/10.1101/2025.11.27.690927 Hayato Yoshioka 1 Graduate School of Agricultural and Life Sciences, Univ. of Tokyo , Tokyo, Japan Find this author on Google Scholar Find this author on PubMed Search for this author on this site ORCID record for Hayato Yoshioka Gota Morota 1 Graduate School of Agricultural and Life Sciences, Univ. of Tokyo , Tokyo, Japan Find this author on Google Scholar Find this author on PubMed Search for this author on this site ORCID record for Gota Morota Hiroyoshi Iwata 1 Graduate School of Agricultural and Life Sciences, Univ. of Tokyo , Tokyo, Japan Find this author on Google Scholar Find this author on PubMed Search for this author on this site ORCID record for Hiroyoshi Iwata For correspondence: hiroiwata{at}g.ecc.u-tokyo.ac.jp Abstract Full Text Info/History Metrics Data/Code Preview PDF Abstract Sustainable improvement of crop performance requires integrative approaches that link genomic variation to phenotypic expression through intermediate molecular layers. Here, we present Reciprocal Best Linear Unbiased Prediction (Reciprocal BLUP), a predictability-guided multi-omics framework that quantifies cross-layer relationships among the genome, metabolome, and microbiome to enhance phenotype prediction. Using a panel of 198 soybean accessions grown under well-watered and drought conditions, we first evaluated four direction-specific prediction models (genome→microbiome, genome→metabolome, metabolome→microbiome, and microbiome→metabolome) to estimate the predictability of individual omics features. Then, we evaluated whether subsets of features with high cross-omics predictability improve phenotype prediction. These cross-layer models identify features that play physiologically meaningful roles within multi-omics systems, enabling us to prioritize variables that capture coherent biological signals enriched for phenotype-relevant information. As a result, metabolome features were highly predictable from microbiome data, whereas microbiome predictability from metabolomic data was weaker and more environment-dependent, revealing an asymmetric relationship between these layers. In the subsequent phenotype prediction analysis, the model incorporating predictability-based feature selection substantially outperformed models using randomly selected features and also achieved prediction accuracies comparable to those of the full-feature model. Under drought, phenotype prediction models based on metabolomic or microbiomic kernels (MetBLUP or MicroBLUP) outperformed the genomic baseline (GBLUP) for several biomass-related traits, indicating that environment-responsive omics layers captured phenotypic variation not explained by additive genetic effects. Our results highlight hierarchical interactions among genomic, metabolic, and microbial systems, with the metabolome functioning as an integrative mediator linking genotype, environment, and microbiome composition. The Reciprocal BLUP framework provides a biologically interpretable and practical approach for integrating multi-omics data, improving phenotype prediction, and guiding omics-based feature selection in plant breeding. 1. Introduction Sustainable increases in agricultural productivity are essential to meet the demands of a growing global population. To address these challenges, modern plant breeding increasingly leverages genomic prediction approaches such as Genomic Best Linear Unbiased Prediction (GBLUP), which enable accurate estimation of breeding values using genome-wide marker information [ 1 ]. While genomic information alone can capture substantial genetic variance, complex traits are also shaped by intermediate molecular and ecological processes that are not fully explained by additive genetic effects [ 2 , 3 ]. Recent advances in multi-omics technologies have made it possible to profile diverse biological layers—genome [ 4 ], metabolome [ 5 ], and microbiome [ 6 ]—within the same individuals. These intermediate molecular layers provide unprecedented opportunities to unravel the complex causal relationships among biological systems, as they can act as mediators linking genotypes to phenotypes [ 7 , 8 ]. Among them, rhizosphere microorganisms play a pivotal role in plant growth and resilience by enhancing nutrient uptake, suppressing pathogens, and improving tolerance to abiotic stresses such as drought [ 9 – 11 ]. Understanding the structure and function of the rhizosphere microbiome is therefore essential for deciphering the mechanisms of plant–microbe symbiosis and their contributions to stress adaptation [ 12 ]. In addition, the metabolome provides a complementary dimension by capturing the dynamic physiological state of plants under varying environmental conditions. Because metabolite levels integrate genetic, environmental, and microbial influences, they can serve as sensitive biomarkers of stress responses and genotype-by-environment interactions [ 13 , 14 ]. Despite these advances, most existing frameworks treat all omics features equally, potentially diluting biologically meaningful signals. While Christensen et al. [ 15 ] and Zhao et al. [ 16 ] focused on how to handle missing intermediate omics information within genetic evaluation frameworks, the distinct question of feature selection—identifying which omics variables are biologically relevant and predictive—has been less explored. The need for feature selection arises in part from the characteristics of intermediate omics layers. These layers are high-dimensional and often contain substantial redundancy and noise, making it challenging to determine which features truly contribute to transmitting biologically meaningful signals toward the phenotype. This motivates the need for biologically informed strategies to prioritize omics features. In this study, we aim to develop a reciprocal multi-omics framework that leverages predictability among genomic, metabolomic, and microbiomic layers to guide feature selection. We quantify interactions and variance components across these interconnected layers in a plant–microbe system, and evaluate whether cross-omics–informed feature subsets enhance phenotype prediction under contrasting environmental conditions. To support this framework, we introduce the concept of cross-omics predictability—the extent to which one omics layer can be predicted from another—as a biologically grounded criterion for feature prioritization. This metric highlights features that occupy functionally meaningful positions within the information flow from genotype to phenotype. Highly predictable features tend to represent coherent biological variation, including upstream–downstream regulatory relationships, stable genetic or ecological influences, and consistent responses to environmental perturbations. Therefore, prioritizing such features provides a principled strategy for dimensionality reduction that enriches for phenotype-relevant signals while mitigating noise, potentially improving phenotype prediction accuracy. In the following sections, we first describe the multi-omics datasets and experimental design employed in this study ( Section 2.1 ). We then introduce our predictability-based integration framework, outlining the implementation of inter-omics prediction and the identification of predictable features ( Section 2.2 ), and describe how these selected omics features were used for phenotype prediction relative to genomic and random baselines ( Section 2.4 ). The results of these analyses are presented in Section 3 , followed by a discussion of their biological implications for understanding plant–microbe–environment interactions and improving genomic prediction frameworks ( Section 4 ). 2. Materials and Methods 2.1. Soybean Multi-Omics Data We analyzed multi-omics datasets collected from a common panel of N = 198 soybean accessions. The panel was drawn from the Global Soybean Minicore Collection [ 4 ]. The datasets included: Whole-genome genotypes (genome) : a 198 × 425,858 matrix, where each column corresponds to a single nucleotide polymorphism marker in the soybean genome [ 4 ]. Rhizosphere metabolome : a 198 × 265 matrix of metabolome features obtained from rhizosphere samples. Each column represents the normalized peak area of a distinct metabolome feature detected by tandem mass spectrometry [ 17 ]. Rhizosphere microbiota profiles (microbiome) : a 198 × 16,457 matrix generated by 16S rRNA gene amplicon sequencing of DNA extracted from rhizosphere samples. Each column represents the relative abundance of an amplicon sequence variant (ASV) [ 17 ]. Plant phenotypes : a 198 × 9 matrix. These traits include biomass-related features such as shoot and leaf fresh/dry weights, growth stage, plant height, stem length, number of nodes, and number of branches. All measurements were obtained from field trials conducted in 2019 under two watering conditions: a well-watered regime (control) and a water-limited regime (drought), enabling assessment of environmental effects. Microbiome and metabolome datasets are identical to those described in Dang et al. [ 17 ], and were collected from the same field plots as the phenotypic measurements to ensure sample consistency across omics layers. Microbiome data were filtered to remove chloroplast and mitochondrial sequences, and all datasets were subsequently scaled prior to downstream analyses. The details of the data collection are described in Appendix A. 2.2. Predictability analysis among omics layers We implemented a Best Linear Unbiased Predictor (BLUP)-based framework to integrate multiple omics layers through predictability-guided feature selection and phenotype prediction. Although BLUP was originally formulated to predict random effects based on expected pedigree relationships, its linear mixed-model framework was later extended to incorporate genome-wide markers (GBLUP). This extension enables efficient analysis of high-dimensional genomic data [ 1 ]. This framework quantifies how omics features represent genetic and environmental signals, identifies stable high-predictability features across environments, and evaluates their contribution to the prediction of phenotypic values. The analysis focuses on interactions among genome, microbiome, and metabolome layers under both control and drought conditions. To characterize how different omics layers capture genetic and environmental information, we performed four direction-specific prediction analyses. We employed BLUP models, each fitted with a single omics-derived relationship matrix ( Fig. 1 ): Download figure Open in new tab Figure 1. Conceptual overview of the modelling framework. We investigated how genomic (G), microbiome (Micro), and metabolome (Met) layers contribute to phenotype prediction across environments. Four direction-specific BLUP models were considered: (1) genome → microbiome, (2) genome → metabolome, (3) metabolome → microbiome, and (4) microbiome → metabolome. The microbiome and metabolome were sampled from the rhizosphere. The omics features were ranked by predictability in these directions, and the selected top features were evaluated in phenotype prediction. Model 1: genome → microbiome Model 2: genome → metabolome Model 3: metabolome → microbiome Model 4: microbiome → metabolome Let M = ( m 1 , …, m q ) denote the omics matrix to be predicted (e.g., metabolome or microbiome), and X the explanatory omics used as predictors. Here, when M represents the metabolome, m j denotes the vector of observed values for the j -th metabolite; whereas when M represents the microbiome, m j denotes the vector of observed abundances for the j -th ASV. For each response vector m j , we fit a BLUP model with a single omics-derived relationship matrix: where α j is the fixed-effect intercept for feature m j , Z is the incidence matrix linking samples to the random effects, a j represents the vector of omics-based random effects, and η j is the residual error term. The variance component denotes the variance of the omics-based random effects, whereas represents the residual variance, both specific to feature m j . The matrix K X represents the covariance structure of the omics-derived random effects. It is constructed from the explanatory omics data X as a linear relationship matrix: where q X denotes the number of explanatory features. Predictive performance was evaluated by 5-fold cross-validation, and quantified by Pearson’s correlation (COR) between observed and predicted values: 2.3. Feature selection based on predictability Features were ranked according to their predictability score r j , which quantifies how well each omics feature can be predicted from other data sources. Instead of applying a fixed correlation threshold, we selected features based on their relative predictability ranking. Specifically, we considered the top 100%, 50%, 25%, 12.5%, and 6.25% of features, corresponding to progressively stricter levels of selection. For each environment e ∈ {control, drought} , the top predictable features were defined as the top p % of features ranked by : Here, p denotes the proportion of features retained. Because a given omics layer may be predictable from multiple information sources, we combined results across relevant models. For the metabolome, three feature selection strategies were applied: top:G — metabolites highly predictable from the genome; top:Micro — metabolites highly predictable from the microbiome; top:G+Micro — metabolites highly predictable from genome and microbiome. Similarly, for the microbiome, features were selected under three analogous rules: top:G — microbes predictable from the genome; top:Met — microbes predictable from the metabolome; top:G+Met — microbes predictable from genome and metabolome. These top-feature sets were used for downstream phenotype prediction analysis. 2.4. Phenotype prediction using selected omics features To evaluate the contribution of predictable omics features to phenotype prediction, we constructed linear kernels—used here as sample-to-sample covariance structure (relationship matrices)—from the features in each omics layer as follows: where M ∗ represents the feature matrix used to compute the kernel and q ∗ is the number of features included. We compared three feature sets: Full omics : all available features in the omics layer were used to construct the kernel, representing the baseline performance without feature filtering. Top omics : only the top p % of features ranked by predictability in environment e were used. These represent the biologically and statistically predictable subset identified through the predictability analysis. Random omics : a random subset of features with the same size as , generated for each iteration index b ( b = 1, …, 10) to benchmark random selection against predictability-based selection. These three settings enabled us to disentangle the effects of feature selection and model complexity on phenotype prediction accuracy. Phenotype prediction was then performed using a BLUP model: Here, y denotes the vector of phenotypic values, μ is the fixed-effect intercept, and Z is the design matrix from samples to random effects. The vector u represents random effects with variance component , whereas ε denotes residual errors with variance . Models based on metabolome and microbiome kernels are denoted as MetBLUP and MicroBLUP , respectively. The GBLUP model using the additive genomic relationship matrix K G , calculated according to Eq. (2) , was included as a baseline. To account for variability in random feature selection, the entire 5-fold cross-validation procedure (using five different seeds) was repeated 10 times, with a new random subset generated at each iteration ( b = 1, …, 10). All analyses were implemented in R [ 18 ] using the RAINBOWR package [ 19 ], and visualizations were produced with ggplot2 [ 20 ]. 3. Results 3.1. Predictability patterns across omics layers and environments Pairwise prediction performance across omics layers is summarized in Fig. 2 . Microbiome prediction was more accurate when metabolome features were used as predictors (Model 3: Met → Micro) than when genome features were used as predictors (Model 1: G → Micro), under both control and drought conditions. Conversely, the prediction of the metabolome was more accurate when microbiome features served as predictors, particularly under drought conditions (Model 4: Micro → Met). Download figure Open in new tab Figure 2. Prediction accuracy for each omics by method and environment. Left: (a) Metabolome; Right: (b) Microbiome. For the microbiome, Model 1 (G → Micro) and Model 3 (Met → Micro); for the metabolome, Model 2 (G → Met) and Model 4 (Micro → Met) are compared under Drought (D) and Control (C) conditions. Diagonal panels show the distributions of prediction accuracy of each omics feature for each setting; lower triangles display pairwise scatter plots of prediction accuracy for each omic feature (the dashed line indicates y = x ); upper triangles report correlation coefficients between settings. The similar pattern was observed in the number of features selected based on predictability thresholds ( Table 1 ). Model 1 (G → Micro) selected far fewer features with r > 0.1, only 10/1193 (0.8%) under control and 35/893 (3.9%) under drought conditions. In contrast, Model 3 (Met → Micro) selected a larger number of metabolome features with high predictability ( r > 0.1), identifying 132/1193 (11.1%) under control and 127/893 (14.2%) under drought conditions. View this table: View inline View popup Download powerpoint Table 1. Number of selected features. The correlation (r) threshold for each direction-specific prediction model. The Full column shows the total number of available features without correlation filtering. Similarly, Model 4 (Micro → Met) under drought conditions identified more predictable metabolome features than under control conditions, suggesting enhanced metabolic signaling to microbial structure in water-limited environments. 3.2. Robust metabolome features across conditions and models A Venn diagram summarizing features exceeding r > 0.1 across the four models and two environments ( Fig. 3 ) revealed a set of robust metabolome features consistently predictable regardless of model direction or environment. Three metabolites—glutamic acid, daidzein, and genistin—were consistently selected across all conditions. These compounds have been implicated in root exudation, plant–microbe interactions, and stress responses in the rhizosphere, underscoring their biological relevance [ 21 – 23 ]. Download figure Open in new tab Figure 3. The selected microbiome and metabolome features under each condition. Drought: yellow; Control: blue. Identified using Pearson correlation ( r > 0.1). (A) For the prediction of the metabolome, genomic or microbiome features were used as inputs. (B) For the prediction of the microbiome, genomic or metabolome features were used as inputs. In contrast, no single bacterial taxon was consistently predictable across environments or models. Microbiome feature selection exhibited pronounced environment-and genome-dependent variation. 3.3. Phenotype prediction using predictability-selected omics features 3.3.1. MetBLUP: Phenotype prediction by selected metabolome We next assessed whether predictability-based feature selection enhanced phenotype prediction. For metabolome-based models (MetBLUP; Fig. 4 ), models based on the most predictable features consistently outperformed those built from randomly selected subsets of equal size, and their predictive accuracy was often comparable to or greater than that of full-feature models. Download figure Open in new tab Figure 4. Phenotype prediction via MetBLUP across selection thresholds. The x-axis shows selection percentages, with the number of selected features at each threshold in parentheses. Solid lines represent models based on metabolome features selected by different methods: “top:G” (features selected using genomic information), “top:Micro” (features selected using microbiome information), and “top:G+Micro” (features jointly selected using both genomic and microbiome data). Dashed gray lines correspond to random feature selections matched by feature count. The dotted red horizontal lines indicate the prediction accuracy obtained by GBLUP, used here as a genomic baseline. For each point, the marker indicates the mean prediction correlation from 5-fold CV with five different seeds. For the random model, this procedure was repeated 10 times; the error bars represent the standard deviation across iterations. For the other models (top:G, top:Micro, top:G+Micro) , the kernel is identical across iterations; accordingly, no error bars are shown. Overall, reducing the feature set from 265 to 34 features (25%) resulted in only a negligible reduction in prediction accuracy compared with using the full 265-feature set. At each feature-count level, the predictability-based subsets consistently outperformed random subsets matched by feature count. Notably, the performance difference between the selected subsets and the random subsets became larger as the number of features decreased. For most phenotypes, MetBLUP outperformed the genomic baseline (GBLUP) under drought conditions, and achieved comparable accuracy under control conditions. When comparing the three feature-selection rules (top:G, top:Micro, and top:G+Micro), distinct trends emerged depending on the trait. The top:Micro performed best for traits with lower genomic predictability (e.g., shoot dry weight), whereas top:G yielded the highest accuracy for traits that were already well predicted by GBLUP (e.g., growth stage). The top:G+Micro rule provided complementary information, generally resulting in the most stable and robust performance across traits. A similar trend was observed when comparing the two conditions, where GBLUP predictability was generally lower under drought, whereas MetBLUP based on top:G+Micro features achieved superior performance for major biomass traits (such as leaf dry weight), highlighting the influence of environmental factors mediated by metabolomic contributions. 3.3.2. MicroBLUP: Phenotype prediction by selected microbiome For microbiome-based models (MicroBLUP; Fig. 5 ), models using top predictable features also significantly outperformed random subsets of equal size. Under drought conditions, models using only the top 6.25% of microbiome features achieved predictive ability comparable to, or even exceeding, those using all features, and were on par with the GBLUP baseline for most biomass-related phenotypes, such as dry weight of leaves. Download figure Open in new tab Figure 5. Phenotype prediction via MicroBLUP across selection thresholds. The x-axis shows selection percentages, with the number of selected features at each threshold in parentheses. Solid lines represent models based on microbiome features selected by different methods: “top:G” (features selected using genomic information), “top:Met” (features selected using metabolomic information), and “top:G+Met” (features jointly selected using both genomic and metabolomic data). Dashed gray lines correspond to random feature selections matched by feature count. The dotted red horizontal line indicates the prediction accuracy obtained by GBLUP, used here as a genomic baseline. For each point, the marker indicates the mean prediction correlation from 5-fold CV with 5 different seeds. For the random model, this procedure was repeated 10 times; the error bars represent the standard deviation across iterations. For the other models (top:G, top:Met, top:G+Met) , the kernel is identical across iterations; accordingly, no error bars are shown. However, under control conditions, predictive ability remained generally lower than that of GBLUP. Notably, for major biomass-related traits such as shoot weight and leaf weight under control conditions, top:Met-based feature selection outperformed top:G, suggesting that microbiome features that are predictable from metabolomic data better capture relevant biological variation in these traits. In contrast, under drought conditions, the difference between the selection rules was smaller. In several cases where the microbiome was strongly filtered (6.25%), top:G achieved slightly higher accuracy than top:Met, indicating a shift in the dominant sources of microbiome predictability across environments. 4. Discussion 4.1. Predictive asymmetry between omics layers Our results reveal a clear asymmetry in predictive power between the metabolome and microbiome layers. The metabolome was highly predictable from microbiome data, while the reverse direction showed weaker and more environment-sensitive patterns ( Fig. 2 ). Moreover, the selected microbiome features were specific to each drought or control condition ( Fig. 3 ). This finding is consistent with the notion that metabolite pools primarily represent downstream outputs of plant physiology, and that these metabolic outputs shape microbial community structure. As a result, the metabolome provides relatively stable biological signals, whereas the microbiome is more sensitive to environmental fluctuations. 4.2. Genetic versus environmental contributions to phenotype prediction The predictability-based feature selection approach provides a unique perspective on how genetic and environmental factors shape phenotypic variation through different omics layers. Traits with high genomic predictability tended to show better phenotype prediction accuracy when using top:G-selected features, particularly in metabolome-based models (MetBLUP) ( Fig. 4 ). This suggests that, for highly heritable traits, the metabolome captures genotype-associated physiological variation. In contrast, traits with lower genomic predictability or stronger environmental influence benefited more from top:Micro feature sets, indicating that metabolomic data also capture environmentally responsive biochemical or physiological processes that are not directly encoded in the genome. 4.3. Environmental modulation of omics–phenotype relationships Under drought stress, genomic predictability (GBLUP accuracy) was generally reduced, whereas the phenotype prediction performance of MetBLUP and MicroBLUP models improved ( Fig. 4 , 5 ). This shift implies that environmental perturbation enhances the relevance of metabolomic and microbial responses to stress, leading to a greater contribution of environment-associated variation. Thus, while genetic control dominates under optimal conditions, environmental modulation of omic layers becomes a major driver of phenotypic variation under stress. 4.4. Hierarchical relationships among genome, metabolome, and microbiome The relationships among the three omics layers imply a hierarchical structure. The genome and metabolome exhibited strong connectivity, consistent with a causal chain from genomic variation to metabolic state to phenotypic variation in biomass-related traits ( Fig. 4 ). In contrast, the microbiome showed weaker direct associations with the genome, as evidenced by the poor performance of genome-based prediction of microbial composition, compared with metabolome-based prediction ( Fig. 3 ). We also observed the limited ability of the top:G-based MicroBLUP model to predict phenotypic traits, particularly under control conditions ( Fig. 5 ). In addition, microbial composition was explained substantially better by metabolomic data than by genomic data. Given that metabolomic features were comparably predictable from both the genome and the microbiome ( Fig. 3 ), this pattern suggests a potential hierarchical structure in which genomic variation shapes metabolomic states, which in turn influence rhizosphere microbial communities. 4.5. Biomarker selection Our reciprocal framework provides a practical approach to prioritizing omics features based on their predictability and leveraging them for phenotype prediction. In the metabolome, the filtered feature sets (25%) performed comparably to the full dataset in the MetBLUP framework, whereas in the microbiome, an even smaller subset of features (6.25%) outperformed the full dataset in the MicroBLUP framework ( Fig. 4 , 5 ). These selected subsets highlight candidate biomarkers in both metabolomic and microbiome layers, improving biological interpretability. For the metabolome in particular, reducing the number of features has practical benefits, as targeted assays can become substantially less costly when only a small set of compounds needs to be measured. 5. Conclusions Our findings highlight that integrating predictability-based feature selection enables the disentanglement of genetic and environmental contributions to complex traits. The proposed framework clarifies the distinct yet interconnected roles of the genome, metabolome, and microbiome layers. Moreover, our approach enables the selection of biologically meaningful features and enhances the performance of phenotype prediction models. Such an understanding is critical for improving multi-omics prediction models and for designing strategies to exploit both genetic and environmental variation in crop improvement. Author Contributions Conceptualization, H.Y., G.M. and H.I.; Data curation, Y.F., and Y.I.; Formal analysis, H.Y.; Supervision, H.I.; Writing – original draft, H.Y.; Writing – review & editing, G.M., H.I. Funding This study was supported by the Japan Society for the Promotion of Science (JSPS) KAKENHI (JP23KJ0506), and Bourses du Gouvernement Français (BGF; 143537P). This work was also supported by Japan Science and Technology (JST) Core Research for Evolutional Science and Technology (CREST; JPMJCR16O2), JSPS KAKENHI (JP22K21352), and the JST ALCA-Next Program 751 (JPMJAN23D1), Japan. Data Availability Statement All source codes and data are available from the repository in GitHub: https://github.com/Yoska393/ReciprocalBLUP . Conflicts of Interest The authors declare no conflicts of interest. Acknowledgments We are grateful to Yushiro Fuji, Yui Nose, and Yasunori Ichihashi at RIKEN for their support in data collection. We also thank the technical staff at the Arid Land Research Center, Tottori University, for their assistance. Appendix A Multi-omics data collection Appendix A.1 Field experiment A total of 198 diverse soybean accessions from the NARO Genebank ( https://www.gene.affrc.go.jp/ ), representing the Global Soybean Minicore Collection [ 4 , 24 ], were cultivated at the Arid Land Research Center of Tottori University in 2019. Field management followed Dang et al. [ 17 ], Toda et al. [ 25 ], Sakurai et al. [ 26 ]. Plants were grown on sandy soil under two watering regimes: a well-watered control and a drought (non-irrigated) condition. Fertilizers were applied prior to sowing, and irrigation was controlled through buried drip tubes under white mulch sheets. Biomass-related phenotypic measurements and rhizosphere root samples were collected from the same plots at the beginning of September to ensure data consistency across omics layers. Root samples were immediately frozen at − 80 ° C for metabolome and microbiome analyses. The soybean genome sequences are available from the DDBJ Sequence Read Archive (BioProject accession: PRJDB7281). Appendix A.2 Phenotypic traits Nine biomass-related traits were measured: dry and fresh weights of shoots and leaves, growth stage, plant height, stem length, number of nodes, and number of branches. Appendix A.3 Metabolome analysis Root metabolomic profiling was conducted following Sawada et al. [ 27 ] and Uchida et al. [ 28 ]. Freeze-dried and powdered root samples were extracted using 80% methanol containing internal standards, analyzed using LC–QqQ–MS (LCMS-8050, Shimadzu, Kyoto, Japan), and quantified in multiple reaction monitoring (MRM) mode. Peak areas were processed using MRMPROBS [ 29 ] and normalized against internal standards. The processed metabolomic dataset is available from the RIKEN DropMet database (ID: DM0071). Appendix A.4 Microbiome analysis Root-associated and rhizosphere bacterial communities were characterized using 16S rRNA gene sequencing. Genomic DNA was extracted following Kumaishi et al. [ 30 ] and Ichihashi et al. [ 31 ], and the V4 region of the 16S rRNA gene was amplified using a two-step PCR protocol with PNA blocking primers [ 32 , 33 ]. Amplicons were sequenced on an Illumina MiSeq platform (2 × 300 bp). Raw reads were processed using QIIME2 (ver. 2020.6.0) with Cutadapt and DADA2 [ 34 ] for quality filtering, denoising, and chimera removal. Taxonomic assignment was performed with the SILVA database (ver. 138) [ 35 , 36 ], and non-bacterial sequences were excluded. Rarefaction was performed ( n = 3000), resulting in 1,642 ASVs. Additional filtering was applied separately for each condition. The final dimensions of the data were 179 × 1193 for the control condition and 182 × 893 for the drought condition. Raw sequencing data are available at DDBJ (BioProject accession: PRJDB37881). Disclaimer/Publisher’s Note The statements, opinions and data contained in all publications are solely those of the individual author(s) and contributor(s) and not of MDPI and/or the editor(s). MDPI and/or the editor(s) disclaim responsibility for any injury to people or property resulting from any ideas, methods, instructions or products referred to in the content. Funder Information Declared Japan Society for the Promotion of Science, https://ror.org/00hhkn466 , JP23KJ0506 , JP22K21352 Japan Science and Technology Agency, https://ror.org/00097mb19 , JPMJCR16O2 , JPMJAN23D1 Footnotes https://github.com/Yoska393/ReciprocalBLUP References 1. ↵ Meuwissen , T. ; Hayes , B. ; Goddard , M. Prediction of the Total Genetic Value Using Genome-Wide Dense Marker Maps . Genetics 2001 . 2. ↵ Vahabi , N. ; Michailidis , G. Unsupervised multi-omics data integration methods: a comprehensive review . Frontiers in Genetics 2022 , 13 , 854752 . OpenUrl PubMed 3. ↵ Bi , Y. ; Yassue , R.M. ; Paul , P. ; Dhatt , B.K. ; Sandhu , J. ; Do , P.T. ; Walia , H. ; Obata , T. ; Morota , G. Evaluating metabolic and genomic data for predicting grain traits under high night temperature stress in rice . G3: Genes, Genomes, Genetics 2023 , 13 , jkad052 . OpenUrl 4. ↵ Kajiya-Kanegae , H. ; Nagasaki , H. ; Kaga , A. ; Hirano , K. ; Ogiso-Tanaka , E. ; Matsuoka , M. ; Ishimori , M. ; Ishimoto , M. ; Hashiguchi , M. ; Tanaka , H. ; et al. Whole-genome sequence diversity and association analysis of 198 soybean accessions in mini-core collections . DNA Research 2021 , 28 , dsaa032 . OpenUrl CrossRef PubMed 5. ↵ Guo , X. ; Sarup , P. ; Jahoor , A. ; Jensen , J. ; Christensen , O.F. Metabolomic-genomic prediction can improve prediction accuracy of breeding values for malting quality traits in barley . Genetics Selection Evolution 2023 , 55 , 61 . doi: 10.1186/s12711-023-00835-w . OpenUrl CrossRef PubMed 6. ↵ Marco , S. ; Loredana , M. ; Riccardo , V. ; Raffaella , B. ; Walter , C. ; Luca , N. Microbe-assisted crop improvement: a sustainable weapon to restore holobiont functionality and resilience . Horticulture research 2022 , 9 , uhac160 . OpenUrl 7. ↵ Weishaar , R. ; Wellmann , R. ; Camarinha-Silva , A. ; Rodehutscord , M. ; Bennewitz , J. Selecting the hologenome to breed for an improved feed efficiency in pigs—A novel selection index . Journal of Animal Breeding and Genetics 2020 , 137 , 14 – 22 . OpenUrl 8. ↵ Yang , Z. ; Zhao , T. ; Cheng , H. ; Yang , J. Microbiome-enabled genomic selection improves prediction accuracy for nitrogen-related traits in maize . G3 Genes|Genomes|Genetics 2023 , 14 , jkad286 . doi: 10.1093/g3journal/jkad286 . OpenUrl CrossRef 9. ↵ Kristin , A. ; Miranda , H. The root microbiota—a fingerprint in the soil? Plant and Soil 2013 , 370 , 671 – 686 . OpenUrl CrossRef 10. Xiong , W. ; Song , Y. ; Yang , K. ; Gu , Y. ; Wei , Z. ; Kowalchuk , G.A. ; Xu , Y. ; Jousset , A. ; Shen , Q. ; Geisen , S. Rhizosphere protists are key determinants of plant health . Microbiome 2020 , 8 , 27 . doi: 10.1186/s40168-020-00799-9 . OpenUrl CrossRef PubMed 11. ↵ Alemu , A. ; Åstrand , J. ; Montesinos-Lopez , O.A. ; y Sanchez , J.I. ; Fernandez-Gonzalez , J. ; Tadesse , W. ; Vetukuri , R.R. ; Carlsson , A.S. ; Ceplitis , A. ; Crossa , J. ; et al. Genomic selection in plant breeding: Key factors shaping two decades of progress . Molecular Plant 2024 , 17 , 552 – 578 . OpenUrl CrossRef PubMed 12. ↵ Ulrich , D.E.M. ; Sevanto , S. ; Ryan , M. ; Albright , M.B.N. ; Johansen , R.B. ; Dunbar , J.M. Plant-microbe interactions before drought influence plant physiological responses to subsequent severe drought . Scientific Reports 2019 , 9 , 249 . doi: 10.1038/s41598-018-36971-3 . OpenUrl CrossRef PubMed 13. ↵ Wen , W. ; Li , D. ; Li , X. ; Gao , Y. ; Li , W. ; Li , H. ; Liu , J. ; Liu , H. ; Chen , W. ; Luo , J. ; et al. Metabolome-based genome-wide association study of maize kernel leads to novel biochemical insights . Nature Communications 2014 , 5 , 3438 . doi: 10.1038/ncomms4438 . OpenUrl CrossRef PubMed 14. ↵ Montesinos-López , O.A. ; Montesinos-López , A. ; Brandon , A. ; Delgado-Enciso , I. ; Chavira-Flores , M. ; Crossa , J. ; Breseghelo , F. ; Dreisigacker , S. ; Sun , J. ; Gerard , G.S. ; et al. Genomic prediction powered by multi-omics data . Frontiers in Genetics 2025 , 16 , 1636438 . OpenUrl PubMed 15. ↵ Christensen , O. ; Börner , V. ; Varona , L. ; Legarra , A. Genetic evaluation including intermediate omics features . Genetics 2021 , 219 . 16. ↵ Zhao , T. ; Zeng , J. ; Cheng , H. Extend mixed models to multilayer neural networks for genomic prediction including intermediate omics data . Genetics 2022 , 221 , iyac034 . OpenUrl PubMed 17. ↵ Dang , T. ; Fuji , Y. ; Kumaishi , K. ; Usui , E. ; Kobori , S. ; Sato , T. ; Narukawa , M. ; Toda , Y. ; Sakurai , K. ; Yamasaki , Y. ; et al. I-SVVS: integrative stochastic variational variable selection to explore joint patterns of multi-omics microbiome data . Briefings in Bioinformatics 2025 , 26 . 18. ↵ R Core Team . R: Language and Environment for Statistical Computing , 2023 . 19. ↵ Hamazaki , K. ; Iwata , H. RAINBOW: Haplotype-based genome-wide association study using a novel SNP-set method . PLoS Computational Biology 2020 , 16 , 1 – 17 . OpenUrl 20. ↵ Wickham , H. ggplot2: Elegant Graphics for Data Analysis ; Springer-Verlag New York , 2016 . 21. ↵ Okutani , F. ; Hamamoto , S. ; Aoki , Y. ; Nakayasu , M. ; Nihei , N. ; Nishimura , T. ; Yazaki , K. ; Sugiyama , A. Rhizosphere modelling reveals spatiotemporal distribution of daidzein shaping soybean rhizosphere bacterial community . Plant, Cell & Environment 2020 , 43 , 1036 – 1046 . OpenUrl 22. Kim , D.R. ; Jeon , C.W. ; Cho , G. ; Thomashow , L.S. ; Weller , D.M. ; Paik , M.J. ; Lee , Y.B. ; Kwak , Y.S. Glutamic acid reshapes the plant microbiota to protect plants against pathogens . Microbiome 2021 , 9 , 244 . OpenUrl CrossRef PubMed 23. ↵ Yang , Q. ; Wang , G. Isoflavonoid metabolism in leguminous plants: an update and perspectives . Frontiers in Plant Science 2024 , 15 , 1368870 . OpenUrl PubMed 24. ↵ Kaga , A. ; Shimizu , T. ; Watanabe , S. ; Tsubokura , Y. ; Katayose , Y. ; Harada , K. ; Vaughan , D.A. ; Tomooka , N. Evaluation of soybean germplasm conserved in NIAS genebank and development of mini core collections . Breeding science 2012 , 61 , 566 – 592 . OpenUrl CrossRef PubMed 25. ↵ Toda , Y. ; Sasaki , G. ; Ohmori , Y. ; Yamasaki , Y. ; Takahashi , H. ; Takanashi , H. ; Tsuda , M. ; Kajiya-Kanegae , H. ; Lopez-Lozano , R. ; Tsujimoto , H. ; et al. Genomic prediction of green fraction dynamics in soybean using unmanned aerial vehicles observations . Frontiers in Plant Science 2022 , 13 , 828864 . OpenUrl PubMed 26. ↵ Sakurai , K. ; Toda , Y. ; Kajiya-Kanegae , H. ; Ohmori , Y. ; Yamasaki , Y. ; Takahashi , H. ; Takanashi , H. ; Tsuda , M. ; Tsujimoto , H. ; Kaga , A. ; et al. Time-series multispectral imaging in soybean for improving biomass and genomic prediction accuracy . The Plant Genome 2022 , 15 , e20244 . OpenUrl CrossRef PubMed 27. ↵ Sawada , Y. ; Akiyama , K. ; Sakata , A. ; Kuwahara , A. ; Otsuki , H. ; Sakurai , T. ; Saito , K. ; Hirai , M.Y. Widely targeted metabolomics based on large-scale MS/MS data for elucidating metabolite accumulation patterns in plants . Plant and Cell Physiology 2009 , 50 , 37 – 47 . OpenUrl CrossRef PubMed Web of Science 28. ↵ Uchida , K. ; Sawada , Y. ; Ochiai , K. ; Sato , M. ; Inaba , J. ; Hirai , M.Y. Identification of a unique type of isoflavone O-methyltransferase, GmIOMT1, based on multi-omics analysis of soybean under biotic stress . Plant and Cell Physiology 2020 , 61 , 1974 – 1985 . OpenUrl CrossRef PubMed 29. ↵ Tsugawa , H. ; Kanazawa , M. ; Ogiwara , A. ; Arita , M. MRMPROBS suite for metabolomics using large-scale MRM assays . Bioinformatics 2014 , 30 , 2379 – 2380 . OpenUrl CrossRef PubMed 30. ↵ Kumaishi , K. ; Usui , E. ; Suzuki , K. ; Kobori , S. ; Sato , T. ; Toda , Y. ; Takanashi , H. ; Shinozaki , S. ; Noda , M. ; Takakura , A. ; et al. High throughput method of 16S rRNA gene sequencing library preparation for plant root microbial community profiling . Scientific Reports 2022 , 12 , 19289 . OpenUrl PubMed 31. ↵ Ichihashi , Y. ; Fukushima , A. ; Shibata , A. ; Shirasu , K. High impact gene discovery: simple strand-specific mRNA library construction and differential regulatory analysis based on gene co-expression network . In Plant Transcription Factors: Methods and Protocols ; Springer , 2018 ; pp. 163 – 189 . 32. ↵ Caporaso , J. ; Lauber , C. ; Walters , W. ; Berg-Lyons , D. ; Lozupone , C. ; Turnbaugh , P. ; Fierer , N. ; Knight , R. Global patterns of 16S rRNA diversity at a depth of millions of sequences per sample . Proceedings of the National Academy of Sciences 2011 , 108 , 4516 – 4522 . OpenUrl Abstract / FREE Full Text 33. ↵ Toju , H. ; Baba , Y. DNA metabarcoding of spiders, insects, and springtails for exploring potential linkage between above- and below-ground food webs . Zoological Letters 2018 , 4 , 1 – 12 . OpenUrl CrossRef PubMed 34. ↵ Callahan , B. ; McMurdie , P. ; Rosen , M. ; Han , A. ; Johnson , A. ; Holmes , S. DADA2: High-resolution sample inference from Illumina amplicon data . Nature Methods 2016 , 13 , 581 – 583 . OpenUrl PubMed 35. ↵ Quast , C. ; Pruesse , E. ; Yilmaz , P. ; Gerken , J. ; Schweer , T. ; Yarza , P. ; Peplies , J. ; Glöckner , F.O. The SILVA ribosomal RNA gene database project: improved data processing and web-based tools . Nucleic acids research 2012 , 41 , D590 – D596 . OpenUrl CrossRef PubMed Web of Science 36. ↵ Yilmaz , P. ; Parfrey , L.W. ; Yarza , P. ; Gerken , J. ; Pruesse , E. ; Quast , C. ; Schweer , T. ; Peplies , J. ; Ludwig , W. ; Glöckner , F.O. The SILVA and “all-species living tree project (LTP)” taxonomic frameworks . Nucleic acids research 2014 , 42 , D643 – D648 . OpenUrl CrossRef PubMed Web of Science View the discussion thread. Back to top Previous Next Posted December 01, 2025. Download PDF Data/Code Email Thank you for your interest in spreading the word about bioRxiv. NOTE: Your email address is requested solely to identify you as the sender of this article. Your Email * Your Name * Send To * Enter multiple addresses on separate lines or separate them with commas. You are going to email the following Reciprocal BLUP: A Predictability-Guided Multi-Omics Framework for Plant Phenotype Prediction Message Subject (Your Name) has forwarded a page to you from bioRxiv Message Body (Your Name) thought you would like to see this page from the bioRxiv website. Your Personal Message CAPTCHA This question is for testing whether or not you are a human visitor and to prevent automated spam submissions. Share Reciprocal BLUP: A Predictability-Guided Multi-Omics Framework for Plant Phenotype Prediction Hayato Yoshioka , Gota Morota , Hiroyoshi Iwata bioRxiv 2025.11.27.690927; doi: https://doi.org/10.1101/2025.11.27.690927 Share This Article: Copy Citation Tools Reciprocal BLUP: A Predictability-Guided Multi-Omics Framework for Plant Phenotype Prediction Hayato Yoshioka , Gota Morota , Hiroyoshi Iwata bioRxiv 2025.11.27.690927; doi: https://doi.org/10.1101/2025.11.27.690927 Citation Manager Formats BibTeX Bookends EasyBib EndNote (tagged) EndNote 8 (xml) Medlars Mendeley Papers RefWorks Tagged Ref Manager RIS Zotero Tweet Widget Facebook Like Google Plus One Subject Area Bioinformatics Subject Areas All Articles Animal Behavior and Cognition (7635) Biochemistry (17697) Bioengineering (13895) Bioinformatics (41951) Biophysics (21456) Cancer Biology (18594) Cell Biology (25520) Clinical Trials (138) Developmental Biology (13381) Ecology (19903) Epidemiology (2067) Evolutionary Biology (24323) Genetics (15612) Genomics (22510) Immunology (17737) Microbiology (40401) Molecular Biology (17183) Neuroscience (88622) Paleontology (667) Pathology (2833) Pharmacology and Toxicology (4825) Physiology (7644) Plant Biology (15158) Scientific Communication and Education (2046) Synthetic Biology (4296) Systems Biology (9825) Zoology (2271)