BiMSGP: A Bidirectional Mamba-Based Model for Scalable and Accurate Genomic Prediction in Plants

1. Introduction

Accurately predicting complex phenotypic traits from genomic data is a central challenge in plant breeding and genetics, with significant implications for accelerating genetic gain and improving crop productivity [1,2,3,4,5,6]. Genomic prediction (GP) uses genome-wide markers to estimate breeding values or phenotypic outcomes before all candidates are phenotyped. Its practical application is complicated by the high-dimensional, low-sample-size structure of many breeding datasets, in which the number of markers greatly exceeds the number of individuals [7,8,9,10].

Traditional statistical approaches, such as genomic best linear unbiased prediction (GBLUP) and ridge-regression BLUP, are computationally efficient and remain strong baselines [11,12,13]. Their predominantly additive formulations, however, can have limited capacity to represent non-linear effects and interactions that contribute to complex traits [14,15]. Comparative studies have consequently explored machine- and deep-learning alternatives, while also showing that performance depends strongly on trait architecture, sample size, and model tuning [9,16,17]. Convolutional models such as DeepGS and DLGWAS learn local marker patterns, and newer systems including DNNGP, SoyDNGP, and DeepCCR extend deep genomic prediction to multi-omics, soybean, and large rice populations [6,18,19,20,21]. Fixed local receptive fields and the cost of modeling long marker sequences nevertheless remain important constraints.

Long-range sequence models have demonstrated the value of integrating distant genomic context, from hybrid convolutional–recurrent networks and long-context regulatory models to genomic language models [22,23,24]. State space models provide another route to long-context modeling: HiPPO and S4 established efficient structured-memory formulations, and Mamba introduced input-dependent selective state spaces with linear sequence complexity [25,26,27,28]. Building on this foundation, we propose BiMSGP (Bidirectional Mamba-based model for Scalable and Accurate Genomic Prediction), a bidirectional Mamba architecture specifically designed for plant genomic prediction.

By processing genomic marker sequences in both forward and backward directions, BiMSGP enables each marker to incorporate contextual information from both sides of the sequence. This bidirectional scan can help the model use long-range marker dependencies; it is not direct evidence of biological epistasis between distant loci. The model retains the linear-complexity scan of Mamba.

The main contributions of this work are as follows:

1.

We propose BiMSGP, a bidirectional Mamba-based architecture tailored for genomic prediction. Ablation studies on Wheat599 show that shallow bidirectional stacks outperform deeper and unidirectional variants. The main comparison tables report the residual-on $L=1$ configuration; a residual-off ablation reached a higher mean PCC on that panel and was not substituted into the main tables.

2.

We evaluate BiMSGP on three datasets of increasing scale using within-dataset 10-fold cross-validation: Wheat599 (599 accessions and 1,279 markers), Wheat2000 (2,000 accessions and 33,709 markers), and SoyBase (trait-specific n up to 15,878 accessions and 42,195 SNPs in the processed matrix). We report numerical ranks together with paired fold-wise tests and do not interpret small PCC differences as significant advantages.

3.

We analyze the accuracy–efficiency trade-off of BiMSGP. The model has far more parameters than lightweight baselines (42M versus 81K for DeepGS on Wheat599) and a higher training cost on Wheat2000 and SoyBase; millisecond-scale inference is the operational advantage after a model has been trained.

The remainder of this paper is organized as follows. Section 2 details the BiMSGP architecture, evaluation metrics, and experimental setup. Section 3 presents the main experimental results, ablation studies, and computational efficiency analysis. Section 4 discusses the key findings, limitations, and implications. Finally, Section 5 concludes the paper and outlines future research directions.

2. Materials and Methods

This section details the datasets, data preprocessing procedures, model architecture, implementation details, and evaluation framework used in the development and assessment of the proposed BiMSGP model.

2.1. Datasets and Preprocessing

BiMSGP was evaluated on three publicly available genomic datasets (Table 1):

• Wheat599: 599 wheat lines genotyped with 1,279 Diversity Arrays Technology (DArT) markers after quality control and grain yield (GY) measured in four environments [6,29].
• Wheat2000: 2,000 Iranian bread wheat lines with 33,709 DArT markers and eight quantitative traits (Grain Length (GL), Grain Width (GW), Grain Hardness (GH), Thousand Kernel Weight (TKW), Test Weight (TW), Sodium Dodecyl Sulphate Sedimentation (SDS), Grain Protein (GP), and Plant Height (PHT)) [18,30].
• SoyBase: a soybean panel derived from the public SoyDNGP/SoyBase resources [20,31]. Gao et al. reported 13,784 accessions and 32,032 SNPs after their complete-case and marker-intersection filters. The processed genotype table used in this study contained 15,878 genotyped accessions and 42,195 SNP columns; we did not restrict the main tables to that 13,784 × 32,032 subset. After excluding accessions with missing values of the target phenotype, sample sizes ranged from 14,485 (Yield) to 15,878 (oil and protein). Seven quantitative traits (Plant Height (Hgt), Oil Content (Oil), Protein Content (Protein), Hundred-seed Weight (SdWgt), Flowering Date (R1), Maturity Date (R8), Yield) and four categorical traits (Flower Color (FC), Pubescence Density (P_DENS), Pod Color (POD), Stem Termination (ST)) were analyzed.

SNP markers (SoyBase) were encoded as {0, 1, 2}, representing the number of minor-allele copies. DArT markers in Wheat599 and Wheat2000 were encoded as binary values (0 = absence, 1 = presence). No additional marker selection was applied beyond the source quality-control procedures described below. Within each training fold, marker values were standardized with scikit-learn StandardScaler fitted on the training partition only and applied to the held-out partition (bimsgp/train.py). Phenotypes for regression were likewise standardized to zero mean and unit variance within each training fold; metrics were reported after inverse transformation to the original scale.

The Wheat599 matrix used for both the main tables and the ablation studies contains 1,279 markers. The source panel comprised 1,447 DArT markers; markers with minor allele frequency below 0.05 were removed and missing genotypes were imputed by sampling from the marginal Bernoulli distribution, following Crossa et al. [29] and Wang et al. [6]. No second, different imputation (for example, population-mean imputation) was applied at training time. Wheat2000 retained 33,709 DArT markers after the quality control of Ma et al. [18]. For SoyBase, genotypes were read from the processed marker table (geno-big.parquet; 15,878 rows and 42,195 SNP columns) without applying the 32,032-SNP intersection or 13,784-accession complete-case filter of Gao et al. Samples with missing values of the target phenotype were then excluded, giving the trait-specific n in Table 1. No additional genotype imputation was applied in the training scripts. The public code expects these already-processed matrices (Geno.csv / wheat_GS_v2.RData / geno-big.parquet); scripts that rebuild the 1,279-marker Wheat599 table or the 42,195-column SoyBase table from raw public downloads are not included in the released tree.

Markers were arranged as a one-dimensional sequence in the column order of the source genotype matrices. Wheat599 columns are DArT identifiers (for example, wPt and c names) and do not include chromosome or physical-position annotation, so the sequence is not guaranteed to follow genomic coordinates. Wheat2000 followed the column order of the genotype matrix in the public wheat_GS_v2.RData object. SoyBase SNPs followed the column order of the processed genotype table. Concatenating markers across chromosomes, or using an unsorted marker list, places the last marker of one group next to the first marker of the next; such junctions are artifacts of serialization rather than biological adjacency. Figure 1 is an abstract schematic of this marker sequence and does not depict chromosome boundaries.

Class counts and proportions for the four SoyBase categorical traits are given in Table 2. Stem termination (ST) is the most imbalanced (class S: 8.1%), which is why balanced accuracy is reported alongside overall accuracy.

All datasets are publicly accessible from their original publications or repositories [6,18,20].

2.2. BiMSGP Model Architecture

The proposed BiMSGP model adopts a bidirectional selective state space modeling approach to effectively capture long-range dependencies among genomic markers while maintaining computational efficiency.

2.2.1. Background on Mamba Architecture

The Mamba model [27] is a recent state space model (SSM) designed as an efficient alternative to Transformers for modeling long sequences. Unlike Transformers, which suffer from quadratic computational complexity with respect to sequence length, Mamba achieves linear time and memory complexity through a selective state space mechanism.

Specifically, Mamba is based on a structured state space model (S4) [26] with the following key innovations:

• Selective State Space: The model parameters (A, B, C matrices) are made input-dependent, allowing the network to dynamically focus on relevant parts of the sequence while filtering out noise. This is particularly beneficial for high-dimensional genomic marker data where only a subset of markers may be strongly associated with the target trait.
• Hardware-aware Algorithm: Mamba employs a parallel scan algorithm that enables efficient GPU-parallel computation, significantly improving training and inference speed compared to standard SSMs.
• Linear Complexity: For a sequence of M markers, Mamba has $\mathcal{O}\left(M\right)$ time and memory complexity, making it particularly suitable for genomic data with tens of thousands of markers ($M=1,279$–$42,195$ in our experiments).

In standard Mamba, information flows only in one direction (forward). To better capture long-range bidirectional dependencies among SNP markers, we propose BiMSGP, which extends Mamba into a bidirectional architecture.

2.2.2. Overall Framework

As illustrated in Figure 1(A), BiMSGP consists of an embedding layer, bidirectional Mamba stacks, and a decoder head. The marker sequence is embedded, processed in the forward and reverse directions, concatenated, and decoded to a trait prediction. Figure 1(B) shows the same pipeline in more detail: a linear embedding maps each marker to $d_{\mathrm{model}}$, residual Mamba blocks (when used) run on the original and flipped sequences, the two representations are concatenated, and a lightweight MLP decoder produces the output. Chromosome boundaries are not marked; the figure is an abstract marker sequence in source-file column order.

2.2.3. Input and Embedding Layer

The input to BiMSGP is a matrix of genomic markers with shape $[B,M,1]$, where B is the batch size and M is the number of markers (sequence length). Each genomic marker is initially treated as a univariate feature.

A linear embedding layer projects the input from 1 dimension to the model dimension $d_{\mathrm{model}}$:

$${\mathbf{X}}_{\mathrm{emb}}=\mathrm{Linear}(1\to d_{\mathrm{model}})·{\mathbf{X}}_{\mathrm{input}}$$

(1)

This produces an embedded representation of shape $[B,M,d_{\mathrm{model}}]$ that serves as input to the subsequent bidirectional Mamba stacks. The embedding layer allows the model to learn a richer representation of the genomic markers before processing them through the state space layers.

2.2.4. Bidirectional Mamba Stacks

The core of BiMSGP is the bidirectional Mamba architecture, which processes the marker sequence in both forward and backward directions to capture long-range dependencies from both ends.

For the forward path, the sequence passes through a stack of Mamba layers with residual connections and layer normalization:

$${\mathbf{X}}_{\mathrm{fwd}}^{\left(l\right)}=\mathrm{LayerNorm}\left({\mathbf{X}}_{\mathrm{fwd}}^{(l-1)}+\mathrm{Mamba}\left({\mathbf{X}}_{\mathrm{fwd}}^{(l-1)}\right)\right),l=1,\dots ,L$$

(2)

where L is the number of layers in the stack. In our experiments, we typically set L to 1 or 2, as deeper stacks did not yield further improvements according to the ablation studies.

For the backward path, the sequence is first flipped, processed by an identical Mamba stack, and then flipped back:

$${\mathbf{X}}_{\mathrm{bwd}}=\mathrm{Flip}\left(\mathrm{MambaStack}\left(\mathrm{Flip}\left({\mathbf{X}}_{\mathrm{emb}}\right)\right)\right)$$

(3)

In configurations with residual connections, a residual connection and layer normalization are applied after each Mamba layer to stabilize training; the ablation study also evaluates configurations without residual connections.

2.2.5. Feature Fusion and Decoder Head

The forward and backward features are concatenated along the feature dimension:

$${\mathbf{X}}_{\mathrm{fused}}=\mathrm{Concat}\left([{\mathbf{X}}_{\mathrm{fwd}},{\mathbf{X}}_{\mathrm{bwd}}]\right)\in {\mathbb{R}}^{B\times M\times 2d_{\mathrm{model}}}$$

(4)

After a final layer normalization, the fused representation is flattened and passed through a lightweight decoder head:

$$\widehat{\mathbf{y}}=\mathrm{Decoder}\left(\mathrm{Flatten}\left({\mathbf{X}}_{\mathrm{fused}}\right)\right)$$

(5)

Two decoder heads are used (bimsgp/model.py). Wheat599 and Wheat2000 use the wide head (three linear layers; first dropout at least 0.15). SoyBase uses the compact head (two linear layers; dropout $0.1$). Let $d_{\mathrm{in}}=2d_{\mathrm{model}}M$ for bidirectional models.

The Wheat (wide) decoder is

$$\begin{array}{cc}\hfill & \mathrm{Linear}(d_{\mathrm{in}}\to 4d_{\mathrm{model}})\to \mathrm{GELU}\to \mathrm{LayerNorm}\to \mathrm{Dropout}\left(0.15\right)\hfill \\ \hfill & \to \mathrm{Linear}(4d_{\mathrm{model}}\to 2d_{\mathrm{model}})\to \mathrm{GELU}\to \mathrm{Dropout}\left(0.1\right)\hfill \\ \hfill & \to \mathrm{Linear}(2d_{\mathrm{model}}\to d_{\mathrm{output}}).\hfill \end{array}$$

(6)

The SoyBase (compact) decoder is

$$\begin{array}{cc}\hfill & \mathrm{Linear}(d_{\mathrm{in}}\to 2d_{\mathrm{model}})\to \mathrm{GELU}\to \mathrm{LayerNorm}\to \mathrm{Dropout}\left(0.1\right)\hfill \\ \hfill & \to \mathrm{Linear}(2d_{\mathrm{model}}\to d_{\mathrm{output}}).\hfill \end{array}$$

(7)

Here, $d_{\mathrm{output}}$ denotes the task-specific output dimension: one output is used for regression, whereas classification uses the configured number of output logits. This design enables BiMSGP to effectively model complex non-linear marker–trait associations from both directions while keeping the model computationally tractable.

2.2.6. Implementation Details

The model was implemented in PyTorch using the official mamba-ssm library. All experiments were conducted on Ubuntu Linux systems with NVIDIA GPUs, including A100 GPUs with 40 GB or 80 GB of memory, an RTX 3090, and an RTX 4090.

Unless otherwise stated, the main-table configuration (Table 3, Table 4, Table 5 and Table 6) used BiMSGP with $d_{\mathrm{model}}=64$, $d_{\mathrm{state}}=64$, $L=1$ bidirectional Mamba layer per direction, residual connections and layer normalization after each Mamba block, AdamW (learning rate $5\times {10}^{-4}$, weight decay $0.01$), ReduceLROnPlateau (patience 10, factor 0.5), and a maximum of 100 epochs. Wheat599 and Wheat2000 used the wide decoder and batch size 8. SoyBase quantitative traits used the compact decoder and batch size 32, except oil and R8, which used $L=2$ and 20 epochs. SoyBase classification used the compact decoder and batch size 1. Decoder dropout is 0.15 then 0.1 in the wide head and 0.1 in the compact head. The released trainer (bimsgp/train.py) runs in FP32; mixed precision (FP16) was used in the original internal runs but is not implemented in the public script (no autocast/GradScaler). This residual-on $L=1$ Wheat599 setting is the configuration reported in Table 3 (42.0M parameters; Environment 1 mean PCC $0.6198$). Ablation experiments on Wheat599 Environment 1 additionally evaluated residual-off and depth/dimension variants; removing residual connections at $L=1$ increased mean PCC to $0.6387$ (Table 8), but that variant was not substituted into the main comparison tables. For regression, continuous phenotypes were standardized to zero mean and unit variance within each training fold.

2.3. Competing Methods

To evaluate BiMSGP, we compared it against four widely used genomic prediction methods: GBLUP, DeepGS, DNNGP, and SoyDNGP. Brief descriptions of each method are provided below.

2.3.1. GBLUP

Genomic best linear unbiased prediction (GBLUP) is one of the most widely adopted linear mixed model approaches for genomic selection [11,12]. It replaces the traditional pedigree-based relationship matrix with a genomic relationship matrix (GRM) constructed from genome-wide markers. The model can be expressed as:

$$\mathbf{y}=\mathbf{X}\beta +\mathbf{Z}\mathbf{g}+ϵ$$

(8)

where $\mathbf{y}$ is the vector of phenotypes, $\beta$ represents fixed effects, $\mathbf{g}$ is the vector of genomic breeding values with $\mathbf{g}\sim N(\mathbf{0},\mathbf{G}{\sigma }_g^2)$, $\mathbf{G}$ is the GRM, and $ϵ$ is the residual error. GBLUP assumes additive genetic effects and has low computational cost, making it a strong baseline for genomic prediction across various species.

2.3.2. DeepGS

DeepGS is a deep learning-based method for genomic prediction that utilizes a deep convolutional neural network (CNN) to extract hidden features directly from raw genotypic data [18]. The model treats SNP markers as a one-dimensional sequence and applies convolutional layers, pooling, and fully connected layers to learn complex non-linear patterns. DeepGS was one of the early attempts to apply deep learning to genomic selection and has demonstrated improved performance over traditional linear models on several crop datasets by capturing non-linear and non-additive marker effects.

2.3.3. DNNGP

Deep Neural Network Genomic Prediction (DNNGP) is a deep learning framework specifically designed for integrating multi-omics data in plant genomic prediction [6]. It employs a hierarchical multilayer neural network with batch normalization, ReLU activation, and early stopping to mitigate overfitting. DNNGP can handle various types of omics data and dynamically learns feature representations from raw inputs. It has shown competitive or superior performance compared to GBLUP and other deep learning methods, particularly on large-scale breeding datasets, while maintaining reasonable computational efficiency.

2.3.4. SoyDNGP

SoyDNGP is a soybean-specific deep learning framework for genomic prediction [20]. It converts standard VCF files into 3D tensor inputs and employs a convolutional neural network (CNN) architecture optimized for soybean breeding traits. The model is designed to be web-accessible and user-friendly, focusing on capturing complex patterns in high-dimensional SNP data for traits such as yield, oil content, and protein content. SoyDNGP has demonstrated strong performance on large soybean panels and serves as a specialized competitor for legume genomic prediction tasks.

2.4. Evaluation Metrics

To rigorously assess the predictive performance of BiMSGP and competing methods, we employed multiple complementary evaluation metrics suitable for both regression and classification tasks in genomic prediction.

For all quantitative traits (regression tasks), the primary metric was the Pearson correlation coefficient (PCC) between predicted and observed phenotypic values, which measures the strength of the linear relationship and is widely used in genomic selection studies because it is invariant to scale. We also reported the root mean squared error (RMSE) and mean absolute error (MAE) to evaluate the accuracy of absolute predictions:

$$\mathrm{PCC}={\displaystyle \frac{{\sum }_{i=1}^n(y_i-\overline{y})({\widehat{y}}_i-\overline{\widehat{y}})}{\sqrt{{\sum }_{i=1}^n{(y_i-\overline{y})}^2{\sum }_{i=1}^n{({\widehat{y}}_i-\overline{\widehat{y}})}^2}}}$$

(9)

$$\mathrm{RMSE}=\sqrt{{\displaystyle \frac{1}{n}}\sum _{i=1}^n{(y_i-{\widehat{y}}_i)}^2},\mathrm{MAE}={\displaystyle \frac{1}{n}}{\displaystyle \sum _{i=1}^n}|y_i-{\widehat{y}}_i|$$

(10)

where n is the number of evaluated samples, $y_i$ and ${\widehat{y}}_i$ are the observed and predicted values, respectively, and $\overline{y}$, $\overline{\widehat{y}}$ denote their means.

All regression metrics were computed on the original (unstandardized) phenotypic scale after inverse transformation.

For the four categorical traits in the SoyBase dataset (flower color (FC), pubescence density (P_DENS), pod color (POD), and stem termination (ST)), we used classification accuracy and balanced accuracy as the main metrics. Accuracy is the proportion of correctly classified samples:

$$\mathrm{Accuracy}={\displaystyle \frac{1}{n}}{\displaystyle \sum _{i=1}^n}\mathbf{1}({\widehat{y}}_i=y_i),$$

(11)

where $\mathbf{1}(·)$ is the indicator function. Balanced accuracy is particularly suitable for potentially imbalanced classes and is calculated as the average of recall across all classes:

$$\mathrm{Balanced}\mathrm{Accuracy}={\displaystyle \frac{1}{C}}{\displaystyle \sum _{c=1}^C}{\displaystyle \frac{{\mathrm{TP}}_c}{{\mathrm{TP}}_c+{\mathrm{FN}}_c}}$$

(12)

where C is the number of classes, ${\mathrm{TP}}_c$ and ${\mathrm{FN}}_c$ are the true positives and false negatives for class c, respectively.

All performance metrics were evaluated using 10-fold cross-validation. Regression used scikit-learn KFold (10 splits, shuffle, random_state=42). SoyBase classification used StratifiedKFold with the same split count and seed. Within each task, BiMSGP and the baseline methods shared those fold assignments so that methods could be compared on identical held-out partitions. For neural-network models, the held-out test partition of each fold was scored at every epoch, and the checkpoint with the highest test PCC (regression) or balanced accuracy (classification) was retained. No inner validation split was created inside the training partition; hyperparameters listed above were fixed rather than re-tuned on the test fold, but epoch selection did use the test partition. This historical protocol can produce an optimistic bias relative to nested cross-validation, in which early stopping would use only data from the training partition. Reported neural-network accuracies should therefore be interpreted as an upper bound under this protocol. GBLUP has no epoch-wise checkpointing and is not subject to the same selection step.

Pairwise differences in PCC were tested on the ten shared folds with a two-sided paired t-test and a two-sided Wilcoxon signed-rank test. We report the mean difference, a 95% confidence interval, Cohen’s paired effect size $d_z$, and both p-values (Appendix Table A1–Table A3). Boldface in the main tables marks the numerically highest mean only and does not imply statistical significance. P-values are unadjusted for the number of trait-by-method contrasts. Fold-level scores were not archived for Wheat599 SoyDNGP or for SoyBase SoyDNGP, so those comparisons remain descriptive (mean ± SD only).

3. Results

In this section, we evaluate the proposed BiMSGP model against GBLUP [12], DeepGS [18], DNNGP [6], and SoyDNGP [20] on Wheat599 [6], Wheat2000 [18,30], and SoyBase [20,31]. All three datasets were used for regression, whereas classification was evaluated only on SoyBase; all evaluations used 10-fold cross-validation with shared fold assignments. We distinguish numerically highest mean PCC from statistically significant paired differences (Appendix Table A1–Table A3). We then report ablation studies and computational cost.

3.1. Performance on Wheat599 Dataset

The Wheat599 dataset consists of 599 historical wheat lines genotyped with 1,279 DArT genomic markers (after quality control and imputation) across four agroclimatic mega-environments [6]. Phenotypic values in each environment were standardized to a zero mean and unit variance. Predictive performance was assessed using 10-fold cross-validation, with the Pearson correlation coefficient (PCC) as the primary metric and root mean squared error (RMSE) as a secondary metric.

BiMSGP had the numerically highest mean PCC in all four environments (Table 3 and Figure 2). In Environment 1, the mean PCC was 0.6198 ± 0.0708 (RMSE 0.8038 ± 0.0591), compared with 0.5720 ± 0.0763 for DeepGS, 0.5543 ± 0.1015 for DNNGP, 0.5234 ± 0.0798 for GBLUP, and 0.5225 ± 0.1209 for SoyDNGP. Paired tests on the shared folds (Table ) supported a BiMSGP advantage over GBLUP in all four environments ($p_t\le 0.006$). Versus DeepGS, the Environment 1 and Environment 3 differences were significant ($p_t=0.039$ and $0.027$), whereas Environments 2 and 4 were not ($p_t=0.083$ for both). Versus DNNGP, Environments 1, 2, and 4 were significant, but Environment 3 was not ($p_t=0.141$). SoyDNGP fold-level scores were not archived for Wheat599, so that comparison remains descriptive.

3.2. Performance on Wheat2000 Dataset

To assess scalability with a larger sample size and more markers, BiMSGP was evaluated on the Wheat2000 dataset, which consists of 2,000 Iranian bread wheat landraces genotyped with 33,709 DArT markers retained after quality control. Eight agronomic traits were analyzed: grain hardness (GH), grain length (GL), grain protein (GP), grain width (GW), plant height (PHT), sodium dodecyl sulphate sedimentation (SDS), thousand kernel weight (TKW), and test weight (TW).

As shown in Table 4 and Figure 3, BiMSGP had the numerically highest mean PCC for six of eight traits (GH, GL, GP, SDS, TKW, and TW). DeepGS had a slightly higher mean PCC for grain width (0.7500 ± 0.0331 versus 0.7464 ± 0.0384) and plant height (0.3279 ± 0.0617 versus 0.3212 ± 0.0674); neither difference was significant (GW $p_t=0.325$; PHT $p_t=0.473$; Table ). None of the six numerical BiMSGP ranks versus DeepGS was statistically significant (all $p_t\ge 0.093$), including close contrasts such as TKW ($\Delta$PCC $=+0.0008$, $p_t=0.837$). Versus GBLUP, BiMSGP was significantly higher for GH, GL, GW, TKW, and TW, but not for GP, SDS, or PHT. Versus DNNGP and SoyDNGP, BiMSGP was significantly higher for all eight traits.

3.3. Performance on Soybean Dataset (SoyBase)

We further evaluated BiMSGP on the larger SoyBase panel using the same within-dataset 10-fold protocol (trait-specific n of 14,485–15,878 and 42,195 SNPs in the processed matrix). Seven quantitative traits (plant height (Hgt), oil content (Oil), protein content (Protein), hundred-seed weight (SdWgt), flowering date (R1), maturity date (R8), and yield (Yield)) and four categorical traits (flower color (FC), pubescence density (P_DENS), pod color (POD), and stem termination (ST)) were analyzed. For quantitative traits, performance was measured by PCC and RMSE; for classification, accuracy and balanced accuracy were reported.

As shown in Table 5 and Figure 4, BiMSGP had the numerically highest mean PCC on five of seven quantitative traits: Hgt (0.8141 ± 0.0135), SdWgt (0.9201 ± 0.0073), R1 (0.6821 ± 0.0114), R8 (0.7846 ± 0.0176), and yield (0.7900 ± 0.0093). GBLUP had a higher mean PCC for oil (0.8825 versus 0.8816) and protein (0.7182 versus 0.7171). Paired tests (Table ) showed that GBLUP was significantly better for oil ($\Delta$PCC $=-0.0010$, $p_t=0.025$), whereas the protein difference was not significant ($p_t=0.706$). BiMSGP was significantly higher than GBLUP for Hgt, SdWgt, R1, and R8; the yield difference versus GBLUP and DeepGS was not significant ($p_t=0.074$ and $0.193$). Versus DeepGS, only Hgt, R1, and R8 were significant; SdWgt, oil, protein, and yield were not.

For the classification tasks (Table 6 and Figure 5), BiMSGP was competitive but was not the top method on every metric. SoyDNGP attained the highest accuracy and balanced accuracy on the binary traits FC and P_DENS, while DeepGS performed best on the three-class trait ST. BiMSGP had the highest balanced accuracy on POD (0.8497 ± 0.0123). Class proportions (Table 2) show substantial imbalance, especially for ST (class S: 8.1%), which is consistent with the gap between overall accuracy (0.8179) and balanced accuracy (0.6695) for that trait. Fold-level classification scores were not archived for all methods, so these ranks are descriptive.

Overall, BiMSGP was competitive for both regression and classification under within-dataset 10-fold cross-validation. Numerical rank over GBLUP or DeepGS should not be read as a uniform significant advantage, particularly for SoyBase oil and protein and for most Wheat2000 traits versus DeepGS.

3.4. Ablation Studies on Model Architecture

To systematically evaluate the contribution of key components in the proposed Mamba-based genomic prediction framework, we performed a comprehensive set of ablation experiments on the Wheat599 dataset. This dataset comprises 599 historical wheat lines from the CIMMYT Global Wheat Program [29], genotyped with 1,447 DArT markers. After removing markers with minor allele frequency below 0.05 and imputing missing values with samples from the marginal Bernoulli distribution (as described in [29]), 1,279 markers remained. This is the same processed matrix used in Table 3; no alternative population-mean imputation was applied at training time. Phenotypic values for Environment 1 were standardized to zero mean and unit variance (mean ± standard deviation: $0.0000$±$0.9992$). All ablation configurations were trained under identical optimization settings: learning rate $0.0005$, weight decay $0.01$, 100 epochs, and batch size 16. Performance was assessed via 10-fold cross-validation using Pearson correlation coefficient (PCC, higher is better), root mean squared error (RMSE), and mean absolute error (MAE).

The ablation study examined five main factors: (i) presence/absence of residual connections, (ii) number of layers (L), (iii) model dimension ($d_{\mathrm{model}}$), (iv) state dimension ($d_{\mathrm{state}}$), and (v) architectural variants (standard bidirectional Mamba vs. forward-only, backward-only, and sum-pooling variants).

Table 7 summarizes the computational cost of different configurations. Trainable parameter counts ranged from approximately 21M (smallest model) to 673M (largest), while average training time for 100 epochs varied from 0.96 minutes (lightweight forward Mamba) to 8.28 minutes (large $d_{\mathrm{model}}=256$ model). Inference time per sample remained efficient (1.64–1.98 ms) across all configurations.

Table 8 reports results for the effect of residual connections and layer depth (fixed $d_{\mathrm{model}}=64$, $d_{\mathrm{state}}=64$). Removing residual connections consistently yielded a modest but stable improvement in PCC across all layer counts, with no degradation in RMSE or MAE. Performance was optimal with shallow architectures ($L=1$–4), while increasing depth to $L=10$ produced a slight decline, suggesting that deeper Mamba stacks may introduce unnecessary complexity for this genomic prediction task.

Table 9 isolates the impact of model dimension, state dimension, and directional variants (all with $L=2$). Increasing $d_{\mathrm{model}}$ from 32 to 128 improved PCC, with the peak at $d_{\mathrm{model}}=128$ (PCC $0.6387$±$0.0835$). Further scaling to $d_{\mathrm{model}}=256$ yielded comparable performance at substantially higher computational cost (parameters increased ∼4×, from 168M to 673M). The state dimension $d_{\mathrm{state}}=32$ performed marginally better than 16 or 64. Among directional variants, the standard bidirectional Mamba had a higher mean PCC than forward-only, backward-only, and sum-pooling configurations, indicating that bidirectional scanning is useful for using marker context from both directions.

Overall, the Wheat599 Environment 1 ablations support shallow bidirectional stacks over deeper or unidirectional variants. Residual-off $L=1$ reached mean PCC $0.6387$ versus $0.6191$ with residuals (Table 8); a similar mean PCC of $0.6387$ also appeared for residual-on $d_{\mathrm{model}}=128$ at $L=2$ (Table 9). Neither peak was used as the main-table configuration. The published comparisons use residual-on $L=1$ and $d_{\mathrm{model}}=64$ (Environment 1 mean PCC $0.6198$).

3.5. Computational Efficiency

We evaluated the computational efficiency of BiMSGP on the Wheat599 dataset ($n=599$, $m=1,279$), which represents a typical small-to-medium-scale genomic prediction scenario. Training time is reported for 100 epochs for the neural-network models, whereas GBLUP time is reported per cross-validation fold.

Dashes for GBLUP indicate that it is not a neural network: there is no comparable trainable-parameter count, training is reported per cross-validation fold rather than per 100 epochs, and a separate per-sample neural inference time was not recorded. As shown in Table 10, BiMSGP used approximately 42 million parameters and 1.52 minutes for 100 epochs on Wheat599, compared with 81K parameters and 0.40 minutes for DeepGS and 10.4M parameters and 0.48 minutes for DNNGP. The Environment 1 mean PCC was numerically higher (0.6198 versus 0.5720 for DeepGS). Inference remained fast at 1.7 ms per sample.

Training times on Wheat2000 and SoyBase were much longer (Appendix B, Table A1): about 174 minutes per 100 epochs on Wheat2000 (1.105 billion parameters) and 462 minutes on SoyBase (691 million parameters), versus minutes for DeepGS. This training cost is a real practical constraint. Millisecond inference is the advantage after a model has been trained, not a claim that BiMSGP is cheaper to fit.

4. Discussion

The results indicate that BiMSGP is a competitive bidirectional Mamba-based model for genomic prediction under within-dataset 10-fold cross-validation. Model ranking depended on the crop, trait, and baseline, which is consistent with broader benchmarking evidence that no single algorithm is universally best [8,9,16,17]. Bidirectional scanning lets each marker use context from both directions and may help the model use long-range marker dependencies; the experiments do not demonstrate biological epistasis between distant loci.

On Wheat599, BiMSGP had the numerically highest mean PCC in all four environments, with the largest contrast in Environment 1 (0.6198 versus 0.5720 for DeepGS). Paired tests supported an advantage over GBLUP in every environment, but not every contrast with DeepGS or DNNGP was significant. With $n=599$, fold-to-fold variation is large, so overlapping intervals should be interpreted cautiously.

On Wheat2000 ($n=2,000$, $m=33,709$), BiMSGP had the numerically highest mean PCC on six of eight traits. Versus DeepGS, none of those six differences was significant, and DeepGS was numerically higher for GW and PHT. Versus GBLUP, several grain-related traits were significant, whereas GP, SDS, and PHT were not. These results support scalability of the architecture to tens of thousands of markers, but they do not support a claim of uniform superiority over strong convolutional or additive baselines.

On SoyBase, GBLUP remained a strong additive baseline. It was significantly better than BiMSGP for oil ($p_t=0.025$) and numerically better for protein ($p_t=0.706$). BiMSGP was significantly higher than GBLUP for Hgt, SdWgt, R1, and R8, while yield was not significant versus GBLUP or DeepGS. Classification performance was competitive rather than uniformly best, and ST in particular is imbalanced. These evaluations are within-dataset cross-validation, not prediction into an independent population, year, or environment.

Ablation studies on Wheat599 Environment 1 showed that bidirectional processing had a higher mean PCC than unidirectional and sum-pooling variants, and that shallow stacks ($L=1$ or $L=2$) with moderate $d_{\mathrm{model}}$ were sufficient. Removing residual connections gave a modest further gain (mean PCC 0.6387 versus 0.6191 at $L=1$); the main tables used the residual-on $L=1$ configuration.

Computationally, BiMSGP is expensive to train relative to DeepGS (42M versus 81K parameters on Wheat599; 174 and 462 minutes per 100 epochs on Wheat2000 and SoyBase). That training cost is a genuine limitation for routine refitting. Millisecond inference is the benefit after a model has been trained.

Several additional limitations should be noted. First, epoch selection used the test partition of each fold, which can inflate neural-network accuracies relative to nested cross-validation. Second, marker order followed source-file column order; chromosome junctions in the serialized sequence are not biological adjacencies. Third, the study does not include multi-omics inputs or independent-population validation. Future work should include nested evaluation, independent breeding populations or environments, model compression, and multi-trait or genotype-by-environment extensions [14,32].

5. Conclusions

BiMSGP is a bidirectional Mamba-based model for genomic prediction of complex plant traits. Bidirectional scanning lets each marker use sequence context from both directions while retaining linear complexity; this is a modeling device for long-range marker dependencies, not a demonstration of biological epistasis. Under within-dataset 10-fold cross-validation, BiMSGP was numerically strongest for many quantitative traits and competitive for SoyBase categorical traits, but several close comparisons with DeepGS or GBLUP were not statistically significant, and GBLUP was significantly better for SoyBase oil. Shallow bidirectional configurations with moderate model dimensions had higher mean PCC than deeper or unidirectional alternatives.

Training BiMSGP is substantially more expensive than lightweight baselines; millisecond inference is the operational advantage after training. Independent-population or multi-environment evaluation, nested cross-validation, and model compression are important next steps.