An integrated sequential deep learning framework with experience replay for large-scale species distribution modeling under data scarcity

Large-scale species distribution models (SDMs) across many taxa are constrained by data scarcity, sampling bias, and environmental heterogeneity. We present an integrated framework combining environmentally informed sampling with sequential deep learning for large-scale multispecies SDMs, demonstrated using 1, 359 food-relevant plant species in South America. A deep feedforward neural network is trained sequentially across species using experience replay to stabilize learning under species imbalance. Sampling integrates historical climate, soil, and land-cover predictors through geographic thinning, ecoregion-constrained augmentation, and background stratification based on multivariate environmental similarity. Performance is evaluated against established SDM algorithms using independent continental and national test datasets, while a future climate projection is used to benchmark model robustness. Results show stable performance across species niche breadths, improved sensitivity for taxa with restricted environmental tolerances, and consistent rankings across scales and climate conditions. However, this stability yields smoother suitability surfaces that may obscure fine-scale niche contrasts. Furthermore, the analyses reveal that algorithmic behavior under environmental novelty is driven primarily by data structure and sampling design rather than the specific type of algorithmic architecture. The framework is suited for large-scale multispecies SDMs under heterogeneous data when predictive stability and cross-scale consistency are prioritized over detailed single-species inference.

Species Distribution Models (SDMs) are essential tools in ecology, biogeography, and conservation (Vasconcelos et al., 2024). SDMs estimate species-environment relationships by contrasting presence records with background points, using environmental predictors to generate spatial predictions of habitat suitability, which may correlate with abundance or indicate persistence (Lee-Yaw et al., 2022). Specifically, predicting the spatial and temporal dynamics of food plants under historical and future climate scenarios is critical for assessing vulnerabilities in global and regional food security (Wheeler and Von Braun, 2013). Furthermore, large-scale modeling of hundreds of rare and native crops is essential to identify resilient candidate taxa for detailed studies on food security and community adaptation to adverse environmental conditions (Guzzon et al., 2021). In this context, South America provides an ideal but highly demanding scenario for evaluating scalable modeling frameworks. The subcontinent encompasses major biogeographic domains, including the Amazon rainforest, the Andes, and the Atlantic Forest, which generate steep climatic, edaphic, and elevational gradients over relatively short geographic distances (Oliveira-Filho et al., 2021). These conditions drive remarkable levels of plant richness and endemism (Stein et al., 2014) but simultaneously exacerbate geographic and environmental sampling biases due to uneven accessibility and research effort, demanding highly robust modeling approaches (Tydecks et al., 2018).

While traditional SDM applications predominantly focus on single species or localized taxa to address specific ecological hypotheses (Tümer et al., 2026), scaling up to multispecies frameworks significantly enhances predictive power by capturing complex macroecological patterns, leveraging shared environmental responses across taxa, and facilitating robust community-level inference (Altwegg and Nichols, 2019). To achieve this, large-scale studies have traditionally relied on global occurrence repositories, such as the Global Biodiversity Information Facility (GBIF), and broad environmental datasets to apply machine learning algorithms like MaxEnt or random forests. However, these conventional methods typically require independent model calibration for each taxon, limiting their ability to transfer ecological information among species and increasing computational overhead when scaled (Vega et al., 2017). Furthermore, their reliability can vary systematically among species with different degrees of environmental specialization, and these models must overcome distinct challenges associated with data availability (Qazi et al., 2022), expert supervision (Boyd et al., 2023), and computational demands (Doser et al., 2023). Recent advances in deep neural networks (DNNs), particularly convolutional neural networks (CNNs), have expanded this capacity by enabling automated feature extraction from high-dimensional remote-sensing data while offering architectural flexibility and transfer learning capabilities for multispecies modeling scenarios (Deneu et al., 2021).

However, scaling DNNs to hundreds of species often necessitates a continual learning approach, where the network incrementally acquires information from new taxa without retraining from scratch on all previously processed data. While this method enables scalability, it introduces additional challenges, such as catastrophic forgetting (Xu and Zhu, 2018), as well as computational overhead and complex spatial dependencies. Although techniques like shared weights reduce computational demands (Alzubaidi et al., 2021) and specific network modules improve spatial invariance (Jaderberg et al., 2016), scaling these models to large-scale multispecies SDMs still requires explicit strategies to isolate the learning process of each species and prioritize overall model stability. These challenges demand tractable approaches that isolate the effects of continual learning and prioritize training stability over pure architectural complexity.

Beyond these architectural and regional considerations, modeling hundreds of species increases the likelihood of sampling and data-driven challenges (Soley-Guardia et al., 2024). Among them, geographic and taxonomic sampling bias is critical, as global biodiversity datasets often contain uneven spatial sampling effort, positional uncertainty, and context-dependent recording practices in occurrence records (Anderson et al., 2020; Rocchini et al., 2011; Schmidt-Lebuhn et al., 2013). A second issue concerns within-model class imbalance, which arises when presences are contrasted against large background or pseudo-absence sets (Benkendorf et al., 2023). A third issue is cross-species imbalance, where common species contribute many more records than rare species (Kellenberger et al., 2024). These imbalances disproportionately affect narrow-niche species, whose restricted environmental tolerances increase sensitivity to spatial bias and extrapolation. Although conventional corrections, including geographic thinning, environmental filtering, and background weighting have been widely applied (Fourcade et al., 2014), large multispecies models are still affected by default parametrization (Valavi et al., 2022) and overly broad calibration regions (Cooper and Soberón, 2018). DNNs can partially address some of these limitations by learning non-linear representations from complex predictors, reducing reliance on manual feature engineering, and enabling scalable training across many species when sufficient data coverage is available (Hu et al., 2025).

However, effective large-scale SDM applications increasingly depend on integrating sampling design, calibration area definition, training strategy, and evaluation diagnostics within a unified modeling workflow, rather than on model architecture alone. While weighting schemes and adaptive loss functions in DNNs may improve optimization for rare species (Zbinden et al., 2024), they cannot compensate for sparse, geographically biased, or environmentally unrepresentative occurrence data, making explicit bias mitigation essential in large-scale SDM applications.

To mitigate these inherent data limitations, DNN frameworks facilitate the transfer and retention of shared representations across species, as mentioned before, through continual learning. As successfully demonstrated across various network architectures including CNNs and RNNs (Yang et al., 2017; Hussain et al., 2018), transfer and continual learning mechanisms are highly effective for scalable multispecies modeling. In this context, pretrained networks can reuse representations learned from densely sampled species to enhance predictions for data-scarce ones (Hirn et al., 2024). Continual learning, implemented through sequential training, provides a principled framework for incrementally updating a single model while retaining previously learned representations (Parisi et al., 2019). In the context of SDMs, this concept translates into progressively incorporating new species data without eroding ecological representations learned from earlier ones, thereby enabling scalable modeling across hundreds of species. Although sequential machine learning approaches have been explored in ecological data analysis (Civantos-Gómez et al., 2021), their application to DNNs for species distribution modeling remains limited, to the best of our knowledge.

As previously noted, sequential training suffers from catastrophic forgetting, which occurs when new data partially overwrite previously learned relationships (Wickramasinghe et al., 2024). To mitigate this effect, different strategies have been proposed, including regularization-based constraints, optimization and representation learning approaches, and experience replay mechanisms (Wang et al., 2024). Among these strategies, experience replay offers a simple yet effective solution inspired by biological memory consolidation (Hayes et al., 2021) that preserves previously learned representations without requiring architectural modifications. For SDMs, this approach involves selectively reusing information from previously modeled species to limit forgetting and maintain stable species–environment relationships during large-scale sequential training. Furthermore, current SDM approaches also face difficulties in controlling extrapolation (Mesgaran et al., 2014) and integrating environmentally informed sampling across the environmental space of occurrences (Varela et al., 2014). These challenges are particularly pronounced under spatial and temporal transfer, where algorithmic behavior and species niche breadth shape model reliability. While environmentally informed sampling operates during data construction, diagnostic tools such as the Multivariate Environmental Similarity Surface (MESS; Elith et al., 2010) are commonly used post hoc to evaluate extrapolation under novel environmental conditions. In large multispecies settings, however, extrapolation diagnostics can be leveraged to iteratively refine sampling strategies, reducing uneven environmental coverage and mitigating the limitations of random data splits and data scarcity, which are known to compromise model reliability and ecological realism (Guillera-Arroita et al., 2015). Without methodological advances that better integrate environmentally informed sampling with diagnostic evaluation and strengthen continual learning, large-scale deep networks remain prone to biased, unstable, and computationally demanding outcomes (Estopinan et al., 2022).

Therefore, this work aims to advance beyond standard DNNs by integrating experience replay for continual learning and MESS-derived metrics to actively guide background sampling, thereby mitigating the high computational overhead associated with sequential CNN architectures (Delange et al., 2021). The proposed framework further seeks to address geographic and environmental biases, providing a tractable and stable alternative for large-scale species distribution modeling. To achieve this, the architecture operationalizes sequential training within a deep feedforward neural network, combining it with species-wise environmental sampling to scale across hundreds of species. The framework incorporates experience replay as a sequential training strategy, enabling the retention of learned distributions while generalizing to new taxa. Additionally, MESS is used to quantify environmental extrapolation and to stratify background sampling across levels of environmental novelty, supporting environmentally informed sampling in large multispecies models. It is applied to a continental dataset from GBIF comprising native and non-native food plants in South America, with species records originating from Ecuador and predictors derived from climate, soil, and land cover variables. These species represent ecologically relevant components of food systems, whose contrasting environmental responses and management systems introduce a non-trivial challenge for large-scale SDMs. The framework is evaluated for its capacity to improve predictive performance, spatial generalization, and ecological plausibility under both historical and future climate projections, relative to conventional SDM algorithms. Independent testing relies on continental-scale out-of-sample GBIF occurrences and an independent national-scale dataset derived from agronomic census and ethnobotanical records. Model evaluation emphasizes discrimination performance and spatial reliability under extrapolative conditions across gradients of species niche breadth and climate projections. This enables trade-offs between predictive accuracy and robustness to be assessed across algorithms. Crucially, demonstrating how environmentally informed sampling and sequential deep learning with experience replay can be combined provides a highly transferable framework. This integrated approach offers practical tools for ecologists across broader fields to overcome sampling bias, catastrophic forgetting, and extrapolation in diverse geographical contexts.

To guide this study, the following research questions were formulated:

How does the proposed sequential DNN framework perform compared with established SDM algorithms in predictive performance and robustness across species?

How do experience replay and architectural design affect training stability and knowledge retention during sequential multispecies learning?

How do SDM algorithms differ in extrapolation behavior and spatial response across species niche breadth gradients under historical and future climate conditions?

Species occurrence and environmental datasets

The processed GBIF dataset (GBIF.org, 2024) comprised 0.97 million occurrence records for 1, 359 selected food-relevant species, spanning 847 datasets across South America (Figure 1A). These records were retained after standardizing scientific names against the GBIF Backbone Taxonomy and filtering for georeferenced preserved specimens and human observations. These species were compiled from agronomic census data and ethnobotanical registers, most of which had representation at parish-level units in Ecuador and across its four continental bioregions (Figure 1B, de la Torre et al., 2008; INEC, 2020; Santos and Acosta, 2023). These sources correspond to the most recent nationally available compilations at the time of data assembly. Dataset-level metadata were compiled to derive citation counts as a proxy for data usage and curation, and datasets with zero citations were excluded prior to modeling (Khan et al., 2021). This curated dataset included 1, 094 native, 189 non-native, 73 endemic, and 3 species of unknown origin, with more than 50 species belonging to Fabaceae, Arecaceae, Melastomataceae, Solanaceae, and Malvaceae families, most of them observed in the Andean and Amazon bioregions (Figure 1C). Records with observation dates prior to 1995, before widespread GPS adoption, were excluded. Records flagged for spatial inconsistencies were also removed, including those with invalid or duplicated coordinates and those identified by the CoordinateCleaner package (Zizka et al., 2019) as country centroids, urban footprints, or permanent surface water bodies. The species exceeding 19, 000 records included: bay cedar (Guazuma ulmifolia), ball moss (Tillandsia recurvata), macaw leaf (Geonoma macrostachys), açai palm (Euterpe precatoria), stilt palm (Iriartea deltoidea) and walking palm (Socratea exorrhiza) (Figure 1D). Occurrence counts per species varied widely, from 1 to 97, 350, with a median of 123. Data contributions came from 298 organizations, with the largest shares from Aarhus University, Medellín Public Utilities, and the Missouri Botanical Garden. The highest record numbers were from Colombia, Peru, Brazil, Ecuador, and Bolivia, each exceeding 10, 000 records and including areas with more than 50 distinct species.

The environmental predictors used for SDMs were downscaled or aggregated to the target pixel resolution (~4.6 km) and included 19 bioclimatic variables (Fick and Hijmans, 2017) and 9 land-use and land-cover variables (Hurtt et al., 2020) for historical (1970–2000) and future (2060–2080) periods. Future bioclimatic variables were derived from the BCC-CSM2-MR model (Wu et al., 2019), under the Shared Socioeconomic Pathway SSP5-8.5, with the same scenario also applied to land-cover variables. This pathway is commonly used in sensitivity analyses because it represents a high-emission scenario associated with strong projected changes in temperature and precipitation, particularly in tropical regions (Eyring et al., 2016). Additional static predictors comprised 9 soil parameters (Hengl et al., 2017), 6 topographic derivatives from a digital elevation model (CGIAR - CSI, 2008; Theobald et al., 2015), and 3 global surface water layers (Pekel et al., 2016). Although non-climatic variables are difficult to project and thus held constant, including them was considered essential to capture baseline physical constraints, despite imposing a structural limitation under future scenarios (Stanton et al., 2012). Specifically, these environmental predictors were chosen to represent abiotic variation across complementary ecological scales. Bioclimatic variables characterize temperature and precipitation regimes that define species climatic niches, while land-use and land-cover variables capture habitat availability and anthropogenic landscape modification (Thuiller, 2024). Soil properties describe edaphic conditions influencing nutrient availability and water retention (Xiong et al., 2017), whereas topographic and surface water variables account for terrain-driven microclimatic variation and hydrological conditions affecting plant establishment and persistence (Malone et al., 2025). Together, these combined predictor groups capture the major environmental gradients shaping the potential distributions of food-relevant plant species across the continent. For both historical and future periods, the 46 environmental predictors were reduced to 10 components using principal component analysis (PCA) with standarization. This procedure aimed to balance model complexity and environmental data structure (Zhang et al., 2023), while facilitating large-scale multispecies SDM modelling. By reducing dimensionality and multicollinearity, it addressed computational constraints without compromising data integrity, retaining 86% of the original variance. PCA was trained on the historical dataset and the resulting transformation was applied to both periods to ensure comparability across scenarios.

Full details of the environmental input dataset, including the specific soil, topographic, land-cover, and water-related predictors, as well as their downscaling procedures and PCA interpretation, are provided in the Supplementary Information (Section S1). Details on computing and software resources, as well as data and programming workflows availability are described in Section S2.

Environmentally informed sampling

This section introduces the environmentally informed sampling framework adopted in this study, outlining its main components and overall rationale. The schematic overview in Figure 2 provides a conceptual description of the workflow.

Occurrence sampling and augmentation

To obtain spatially independent and environmentally representative samples from species occurrences, a two-step sampling framework was implemented, combining geographic thinning and environmental proximity strategies. Occurrence points were intersected with reduced environmental predictors represented by ten PCA dimensions derived from the full environmental dataset, and sampling was conducted separately for each species.

In the first step, a greedy spatial thinning algorithm was applied to reduce point density while maintaining a minimum separation distance between selected records (Kiedrzyński et al., 2017). For each species, the algorithm starts from a randomly selected occurrence and iteratively adds the point that maximizes its minimum geographic distance to all previously selected records, while discarding candidate points located within the specified distance threshold, promoting spatial balance and increased environmental variability (Aiello-Lammens et al., 2015). A minimum separation distance of 5 km was enforced, exceeding the target pixel resolution, and the sample size was limited to 50 records per species, based on sensitivity analyses (van Proosdij et al., 2016) and computational baseline to ensure equal batch representation across all taxa. Although thinning improved spatial independence, it constrained sample availability, particularly for rare species.

The second step addressed species with insufficient records by applying an environmentally constrained feature-space augmentation using the RESOLVE Ecoregions dataset (Dinerstein et al., 2017). In this context, augmentation functions as an ecological data enrichment process that synthesizes novel representative points within analogous environmental spaces, a strategy recommendation in SDMs to resolve severe data scarcity (Robinson et al., 2018). Rather than assigning false taxonomic labels, we selected neighbor-species, defined as environmentally analogous taxa sharing similar niches, or in their absence, background cells within ecologically analogous ecoregions. These were identified via a k-nearest neighbor search in the PCA-transformed environmental space of the target species (Arya et al., 1998). This strategy acts as a structural prior to mitigate sampling bias rather than biological augmentation, ensuring that species–environment response surfaces remain consistent and representative of the full breadth of the species’ realized niche, even when the density of observed records is low (Broussin et al., 2024).

This two-step framework yields a final set of normalized occurrences labeled as original, neighbor-species, or cell-based samples, offering an ecologically consistent extension of target species data for subsequent modeling steps and evaluation.

Background sampling strategy

To generate background samples across different levels of environmental similarity, we adopted the logic of model-agnostic extrapolation diagnostics (Velazco et al., 2024), which prioritize the detection of non-analog environmental conditions. Consequently, this procedure was implemented to first derive MESS for each target species, where each raster cell stores the minimum similarity score across all environmental predictors. We retained only negative MESS values, corresponding to extrapolative environmental conditions, while positive values indicating interpolative space were masked out. The remaining cells were then classified into four ordinal bins based on empirical quantiles of the negative MESS distribution. Lower quantiles represented stronger extrapolation, whereas higher quantiles corresponded to progressively less dissimilar environmental conditions, with the lowest quantile bin capturing the most extrapolative environments.

Based on the MESS stratification, candidate points from other species were intersected with the classified MESS raster to assign them to extrapolation bins. To prioritize sampling in the most extrapolative regions, a total of 5, 000 background samples were selected to provide sufficient spatial coverage without overwhelming the optimization process, following recommendations for large-scale, presence-only modeling (Phillips et al., 2009). These samples were allocated proportionally across the four MESS bins, with decreasing shares from the most to the least extrapolative bin, assigning 40%, 30%, 20%, and 10% of samples to bins one through four, respectively. When available candidates within a bin were insufficient, remaining samples were filled by randomly selecting raster cells from the same bin. This sampling approach allowed models to discriminate between the core interpolative niche and increasingly novel environments, improving model stability across non-analog climate scenarios and preventing unconstrained extrapolation (Brodie et al., 2022).

Together, geographic thinning, ecoregion-constrained augmentation, and MESS-based background sampling define a species-specific accessible calibration area. This reflects the set of locations and environments that were plausibly accessible to the species within its biogeographic context, as described by Barve et al. (2011).

Sequential deep neural network framework

The implementation of the sequential deep neural network (sDNN) framework is introduced here, encompassing model architecture, training strategy, transfer learning, and spatial inference. Figure 3 provides an overview of the sDNN components, which are detailed in the following subsections.

To configure the sDNN and mitigate catastrophic forgetting across species, we conducted preliminary training experiments. The purpose of these tests was to determine the optimal hyperparameter balance for transfer learning stability. The evaluation pipeline involved systematically varying architectural capacity, experience replay buffer sizes, learning rates, and regularization parameters, as detailed in the following subsections. By iteratively training across randomized species sequences, we established the specific constraints required to maintain generalization without degrading previously acquired representations.

sDNN architecture

The sDNN was implemented with a fully connected architecture using the Keras API (Chollet, 2015). The model received PCA-transformed environmental predictors as input (excluding geographic coordinates) and consisted of five fully connected hidden layers with decreasing widths: 320, 160, 80, 40, and 20 units. These architectural values, totaling 71, 841 parameters, along with other configurations described in this section, were empirically selected based on preliminary experiments in which each parameter was evaluated independently during training while keeping the others fixed (see Section 3.3.). This design reflects common practice in SDMs using tabular environmental predictors, where species–environment relationships are hierarchically structured from broad climatic gradients to finer-scale factors (Albert et al., 2010), and implemented here through progressively narrower network layers while limiting over-parameterization under moderate sample sizes. In addition, this configuration balances model complexity and performance, while keeping the sDNN comparable to the benchmark algorithms (see Section 2.4.1) and other configurations under a common training protocol.

Across this architecture, a consistent set of activation and regularization strategies was applied throughout. All layers applied the rectified linear unit (ReLU) activation function, given its computational efficiency and effectiveness in avoiding vanishing gradients (Ramachandran et al., 2017) but also for its suitability for approximating non-linear and threshold-like ecological responses commonly observed in species–environment relationships. To reduce overfitting, the first two layers used L2 regularization (λ = 0.005 and 0.001), and all layers were followed by dropout layers with progressively decreasing rates (0.3, 0.2, 0.1, 0.05), forming a regularization funnel. This combination of regularization methods is widely used to control overfitting and improve generalization in deep neural networks (Moradi et al., 2020). Complementing this regularization, the output layer consisted of a single sigmoid-activated unit, consistent with a binary presence–background classification setting in SDMs. This configuration acts as a structural bottleneck, ensuring a constant parameter footprint and forcing the network to synthesize generalized environmental features across all taxa (Ven and Tolias, 2019). These shared representations serve as a foundational feature extractor for the subsequent independent fine-tuning phase (Section 2.3.3), where the output is specifically calibrated for species-specific predictions. Binary cross-entropy was used as the corresponding loss function, which is the standard choice for classification tasks (Hui and Belkin, 2021) and is widely adopted in presence–background SDM formulations. It is defined as:

where, y i ∈ { 0, 1 } is the true label and y i ^ ∈ (0, 1) is the predicted probability for the positive class (i.e. predicted distribution).

Model parameters were then optimized under this architectural and regularization configuration, using the Adam optimizer for its proven effectiveness across a wide range of neural network architectures (Li et al., 2022). Learning rates of 1.0 × 10– 5 and 1.0 × 10– 4 were selected from a range of tested values for initial training and fine-tuning, respectively. These lower learning rates were preferred to ensure stable convergence during sequential training across species and to avoid abrupt shifts in learned species–environment relationships. Together, these settings enabled stable convergence during replay-based sequential learning while allowing faster adaptation during final-layer adjustments.

Sequential model training and experience replay strategy

To prepare the training and validation datasets for each species, occurrence records and background samples were merged and labeled as binary classes (1 = presence, 0 = background). Because the sDNN was trained as a supervised classifier, strong class imbalance between presence and background samples can lead to unstable gradients and biased optimization. Therefore, occurrence records were randomly resampled strictly for computational class balancing by applying Gaussian jitter with a mean of 0 ± 0.1 (hereafter, the ± symbol denotes standard deviation) exclusively in the PCA-transformed environmental space. This replication was continued until the number of presence samples matched the background sample size (i.e. 5000 samples), with the sole purpose of stabilizing network training (Buda et al., 2018) rather than performing a second ecological augmentation. By confining this process strictly to the environmental feature space, we prevent the alteration of spatial structure, the introduction of incompatible ecological conditions, and the inflation of autocorrelation. Data were then split into training (80%) and validation (20%) sets using random sampling stratified by class (presence and background), with spatial structure handled during data preparation (see Section 2.2). The validation subset was used solely as a diagnostic set to monitor the area under the receiver operating characteristic (ROC) curve (AUC) and to inform decisions during the construction of the sDNN, including the detection of catastrophic forgetting.

Model training used the historical PCA-transformed predictors and was conducted in sequential batches using experience replay, a technique that reuses previously stored samples to decorrelate training data and mitigate catastrophic forgetting (Moradi et al., 2020). In an ecological context, experience replay acts as a regularization mechanism rather than a biological assumption of niche similarity, allowing the model to retain stable species–environment response functions at the shared environmental feature level across sequential species. This process was implemented by iterating sequentially over each species in a randomized order, which was performed once and kept constant across all experimental runs to ensure comparability. While presentation order can influence knowledge retention and the severity of catastrophic forgetting (Bell and Lawrence, 2022), curriculum learning orderings were excluded to maintain the scope of this research. For each iteration i ∈ { 1, …, T }, where T is the total number of species, the species-specific dataset was denoted as D i. To retain information from previously trained species, a replay buffer R was maintained. At each iteration, the effective training dataset was constructed by combining the current species data with the replay buffer, except for the first iteration:

After training on D i t r a i n, a random subset of 500 samples from D i was added to R. If the buffer exceeded a size of 4000 samples, it was randomly trimmed to maintain this limit. This buffer size (i.e., 40% of the training set) was selected empirically after testing different sizes and also guided by the findings of Fedus et al. (2020), who showed that larger replay capacities can improve performance by moderating fluctuations in the data distribution and enhancing training stability.

Because the incorporation of replay data can reintroduce class imbalance across iterations, class weights were computed dynamically at each training step as the inverse of class frequencies:

where, w c ​ denotes the weight for class c ∈ { 0, 1 } and n c is the number of samples belonging to class in D i t r a i n. Under this weighting scheme, each model was trained for up to 50 epochs with a batch size of 256, a relatively small-batch regime given the tabular structure of the input data, and recommended to improve generalization (Keskar et al., 2017). Early stopping was applied based on validation accuracy to prevent overfitting and to stabilize training during sequential updates. Training was halted when validation accuracy failed to improve by at least 0.001 over five consecutive epochs or when it reached 0.95. Accuracy was computed by binarizing predicted probabilities at a fixed threshold of 0.5, assigning class 1 to probabilities ≥ 0.5.

Model fine-tuning and spatial prediction

Following the initial phase, a second model run was performed across all species using the previously trained model. In this step, the model was fine-tuned independently for each species without sequential learning. This explicit separation aligns with decoupled representation learning frameworks (Kang et al., 2020), mitigating the bias toward abundant taxa inherent in joint optimization by isolating the extraction of shared features from species-specific calibration. While keeping the same batch size and number of epochs, specific adjustments were introduced to maximize accuracy, including a stricter early stopping criterion to halt training when validation accuracy reached 0.99 and extending the patience to ten consecutive epochs. After training again with the D i t r a i n set, the fine-tuned model was used to generate spatial predictions for each species using the PCA-transformed historical and future environmental raster stacks. Following inference, the predicted suitability values were rescaled to the [0, 1] range, as the output from the sigmoid activation function did not consistently occupy the full scale. This behavior is mainly related to experience replay, which favors conservative cross-species representations and prevents extreme probability values. Therefore, a per species min–max normalization was applied to align the effective output distribution with the standardized probabilistic outputs generated by the benchmark classifiers, ensure consistency across models, and facilitate subsequent validation across benchmark algorithms and species groups.

Comparative validation across algorithms and datasets

Benchmark algorithms

In addition to the sDNN, a suite of benchmark SDM algorithms was evaluated for comparison. These included Maximum Entropy (MaxEnt; Phillips et al., 2006) for its robustness under data scarcity, and Generalized Linear Models (GLM) as an interpretable parametric baseline. We also included high-performance tree-based models, namely Random Forest (RF; Breiman, 2001), Gradient Boosting (GBM; Friedman, 2001), and Extreme Gradient Boosting (XGBoost; Chen and Guestrin, 2016), given their capacity to handle complex, non-linear ecological data (Elith et al., 2008). This selection ensures a robust evaluation of the sDNN against both established ecological benchmarks and state-of-the-art machine learning methods.

Unlike the single sequential architecture of the sDNN, these benchmark algorithms were trained independently to construct a distinct classifier for each species. To ensure comparability, all models were trained using the same 80/20 data split and the standardized, PCA-transformed environmental predictors detailed in Section 2.1. Utilizing orthogonal principal components intrinsically addressed multicollinearity assumptions, particularly for parametric approaches like GLM. Applying this identical feature space to the remaining non-parametric algorithms (MaxEnt and tree-based models) did not violate their statistical premises, as they lack formal requirements for data normality or linearity. The AUC on the internal validation subset was used exclusively for hyperparameter configuration. Subsequently, the optimized models generated spatial predictions, which were evaluated using sensitivity and complementary metrics (see Section 2.4.3).

For the supervised classification models used in this comparison (GLM, RF, GBM, and XGBoost), the exact same computational class-balancing resampling applied to the sDNN was implemented. Gaussian jitter was utilized for these baseline classifiers strictly to maintain a standardized training protocol across all algorithms, ensuring that performance comparisons reflect algorithmic differences rather than discrepancies in data delivery. This standardized resampling also mitigates algorithmic bias toward the majority background class during model fitting (Sidumo et al., 2022; Zhang et al., 2022). Model hyperparameters were optimized using five-fold cross-validation, with one, four, and seven tunable parameters for RF, GBM, and XGBoost, respectively, while GLMs did not involve hyperparameter tuning beyond standard model fitting. Optimization targeted parameters controlling model complexity, learning dynamics, and subsampling or regularization, depending on the algorithm.

MaxEnt was treated separately, as it is a presence–background method that does not rely on class-balanced optimization. Therefore, it was trained using the original occurrence data following standard sampling practices (Ahmadi et al., 2023), reflecting the classical SDM performance. We adopted this approach because MaxEnt’s excessive computational time with large background datasets, particularly across hundreds of species, necessitates a distinct implementation from our scalable deep learning and tree-based pipelines. Models were calibrated using the ENMeval R package (Kass et al., 2021) through the implementation with spatial block partitions, optimizing two feature class settings (L, LQ) and regularization multipliers ranging from 1 to 5 through their combined exploration.

Species differentiation by niche breadth

To strengthen validation and test the robustness of SDM algorithms, species were differentiated by their environmental tolerances, allowing performance to be evaluated across groups from ecological broad to narrow niche breadths. The latter was estimated by assessing environmental coverage with PCA-transformed predictors, following Zhigila et al. (2025), using training occurrence records after normalizing the sample size as described in Section 2.2.1 to quantify environmental variability consistently across species. Variability was calculated as the sum of variances across all principal components extracted at occurrence locations:

where, X j and p represent the values of the j t h principal component for a given species and the number of retained components (here p = 10), respectively. Based on this continuous metric, species were grouped into four niche breadth quantiles (Q1–Q4). Higher values reflected broader environmental variability and a greater number of represented ecoregions (Q3–Q4), whereas lower values indicated narrower tolerances and more restricted environmental conditions (Q1–Q2).

Test datasets and evaluation metrics

For final testing purposes, two independent datasets were prepared. The first consisted of a continental-scale GBIF out-of-sample dataset, including 833 species with a mean of 293 ± 301 original occurrence records per species. These data were strictly independent from model calibration, as only non-augmented records not used during training or validation diagnostics were retained, with overlapping geographic coordinates explicitly excluded. The second dataset consisted of a national-scale dataset from Ecuador, constructed at the parish level by integrating agricultural census information and census-derived ethnobotanical records, covering 847 species (see Figure 1B). Parish-level presences were defined using independent evidence from these sources. Model predictions were evaluated by randomly sampling up to ten raster cells within each presence-designated parish to extract suitability values, thereby standardizing contributions across spatial units and complementing point-based occurrence testing. Quantitative test metrics were calculated only for species represented in these external datasets, while species lacking independent test records were excluded from metric-based evaluation but retained for spatial prediction analyses. Historical scenario predictions were extracted at corresponding locations from both datasets to compute probabilistic and threshold-based evaluation metrics:

Sensitivity, evaluated at fixed probability thresholds τ p r o b = { 0.6, 0.7, 0.8, 0.9 }, was used to quantify the proportion of observed presences correctly identified as suitable across a range of binarization decisions, from permissive (0.6) to more restrictive (0.9) distribution limits. Evaluating across this range ensures standardized algorithmic comparability and emphasizes ecologically critical omission errors.

Boyce Index (Boyce et al., 2002), was used to evaluate the consistency between predicted suitability values and observed presences across environmental gradients. This metric is particularly suitable for presence-only and presence–background SDMs, as it assesses whether areas predicted as more suitable are disproportionately associated with observed occurrences, without relying on an explicit absence definition.

Together, these metrics capture both continuous and binary aspects of model performance. Sensitivity and the Boyce Index were prioritized to avoid mischaracterizing highly suitable, unsampled habitats as false positives, which is a known limitation of absence-reliant metrics in presence-only models (Leroy et al., 2018). τ p r o b The Boyce Index ranges from −1 (counter-predictions) to 1 (perfect consistency), with 0 indicating random expectation. For some species, the Boyce Index was undefined due to low contrast or near-uniform suitability predictions. These cases were excluded from Boyce-based analyses and reported separately as a diagnostic indicator of non-informative prediction surfaces rather than model failure. Sensitivity ranges from 0 to 1, indicating the proportion of observed occurrences correctly predicted as suitable.

Additional custom metrics based on the MESS were applied to both historical and future periods using their corresponding spatial predictions:

Environmental Reliability Index (ERI) was defined as the proportion of predicted presence pixels occurring under environmental conditions similar to those observed during model training. This metric was designed to quantify the extent to which predicted suitable areas fall within environmentally known conditions, thereby providing an explicit measure of extrapolation risk in spatial projections. Two MESS-based extrapolative thresholds were considered: τ M E S S > 0, restricting predictions to conditions fully represented in the training data, and τ M E S S > − 10, an empirical buffer optimized to the minimum sufficient threshold, ensuring all evaluated algorithms can express representative distribution surfaces even when the strict zero threshold yields severely truncated areas. This allows moderate extrapolation into proximate novel conditions (Frans et al., 2022). Spatial predictions were binarized using probability thresholds τ p r o b, and ERI was computed as:

where, N MESS > τ M E S S is the number of predicted presence pixels in environmental reliable areas, and N τ p r o b is the total number of predicted presence pixels after probability thresholding.

Extrapolative presence frequency (EPF), was defined as a spatial index quantifying how often a pixel was classified as extrapolative predicted presence across probability thresholds τ p r o b and climatic periods. Extrapolative predicted presence was defined as thresholded predictions occurring under τ M E S S < 0, enabling the identification of predictions outside the training domain. The index captures the spatial persistence of extrapolative predictions and was computed as:

where F h i s (x; τ p r o b) and F f u t (x; τ p r o b) represent the proportion of probability thresholds for which predicted presence at pixel x occurred under extrapolative environmental conditions, as identified by the MESS criterion, in the historical and future periods, respectively.

ERI and EPF were calculated using the corresponding prediction layers and a common MESS surface derived from the training data. For ERI, values close to 1 indicate predictions made under environmentally reliable conditions, reflecting higher confidence and lower extrapolation, whereas values close to 0 indicate the opposite. Similarly, EPF ranges from 0 to 1, where values near 0 indicate pixels that are rarely extrapolative, and values near 1 indicate pixels that are consistently extrapolative across probability thresholds and climatic periods.

Comparisons across algorithms were performed as a function of probability thresholds or niche breadth quantiles using median sensitivity values. Statistical differences were evaluated through iterative pairwise comparisons, where all possible two-sample combinations of algorithms were independently tested using the Wilcoxon signed-rank test for paired samples (Demšar, 2006). To rigorously control for Type I error rates across these multiple comparisons, all resulting p-values were adjusted using the Holm method. Non-significant results (adjusted p > 0.05) were explicitly noted to indicate cases where no statistical differences were detected between the sDNN and other algorithms.

Visualization and future scenarios change

To visualize spatial predictions of suitability change and extrapolative areas, ensembles were derived separately by averaging algorithm suitability probability layers and EPF maps across niche breadth quantiles. In this context, mean aggregation was used to enhance comparability across algorithms, as it retains information from models with more extreme predictions (Hao et al., 2019). For suitability probability ensembles, layers were binarized using an adjusted threshold of τ p r o b = 0.3. This adjustment corrects the statistical smoothing inherent in mean aggregation, which reduces prediction variance and would otherwise inflate omission errors if stricter validation thresholds were applied (Thuiller et al., 2009). The resulting binary historical and future layers were then combined to generate spatial maps of predicted distribution change, distinguishing areas of gain, loss, and persistence. These ensemble predictions were further compared against the mean ensemble using Spearman correlation on continuous EPF values and Cohen’s kappa on categorical distribution change classes to assess deviations from the consensus pattern.

To highlight the magnitude and direction of changes in future predictions relative to historical conditions, a relative change index was calculated as:

where P f u t u r e and P h i s t o r i c a l represent the future and historical suitability-probabilities averages derived for each niche breadth quantiles, and ϵ is a constant with a value of 0.01 used to avoid division by zero. The result of this operation is a dimensionless ratio, where values between −0.1 and 0.1 indicate little or no change between future and historical predictions. Negative ratios between −0.5 and −0.1 correspond to decreases of up to 50% relative to the historical value, whereas values below −0.5 indicate decreases greater than 50%. Positive ratios between 0.1 and 0.5 represent increases of up to 50%, and values above 0.5 indicate increases exceeding 50%. In the results section, these ratios are reported as aggregated positive and negative change rates, which were mapped using a divergent color palette and compared with the mean ensemble representative through Spearman correlation.

Niche breadth stratification and sampling composition

Across the niche breadth gradient defined by Q1–Q4, species ranged from those with environmentally narrow tolerances to those with broadly distributed environmental profiles. Within this gradient, native species dominated all quantile classes (> 72%), with endemic species primarily concentrated in Q1 (12.3%), whereas non-native species were more prevalent in Q4 (24.9%). The proportion of augmented samples declined progressively from narrow to broad niche breadth quantiles, reflecting higher data scarcity and spatial clustering characteristic of species with restricted environmental tolerances. As a result, predictions for species in Q1 were based on a larger contribution of ecoregion-constrained augmented samples relative to original occurrences. Table 1 summarizes the distribution of species across niche breadth classes, including species counts, native species prevalence, environmental variability, ecoregion coverage, and the proportion of augmented occurrence samples.

Predictive performance across benchmark algorithms and validation scales

Evaluation using the continental-scale test dataset revealed clear differences among algorithms. sDNN achieved the highest median sensitivity (0.86). GBM, GLM, and XGBoost showed comparable intermediate performance (0.80 to 0.85), whereas RF and MaxEnt consistently ranked lowest, both below 0.70 (Figure 4A).

At the national scale, sensitivity declined, and rankings shifted moderately. GLM achieved the highest sensitivity (0.75), followed by an intermediate group including sDNN and GBM (0.61 to 0.70), while XGBoost, RF, and MaxEnt again occupied the lowest positions (Figure 4B). Despite these changes, sDNN remained among the top-performing algorithms across both spatial scales.

Across individual probability thresholds, algorithm rankings were largely consistent with aggregated results. At the continental scale, sDNN and GBM showed high and stable sensitivity (>0.79) with no statistically significant differences (adjusted p > 0.05) at the 0.8 and 0.9 thresholds, whereas GLM and XGBoost displayed greater variability and RF and MaxEnt declined sharply at restrictive thresholds (Figure 5A).

At the national scale, performance ranges overlapped among sDNN, GBM, and GLM (0.44 to 0.90), but differences among thresholds were more variable. No statistically significant differences (adjusted p > 0.05) were observed between GLM and sDNN at the 0.8 threshold. In both datasets, sensitivity increased at lower thresholds and declined as thresholds became more conservative (Figure 5B).

When performance was stratified by niche breadth categories, algorithm rankings remained consistent across narrow to broad quantiles. At the continental scale, sDNN and GBM achieved the highest sensitivities for narrow quantiles (>0.85), followed by GLM and XGBoost, while RF and MaxEnt consistently ranked lowest. Sensitivity declined moderately for broad quantiles across all algorithms (Figure 6A).

At the national scale, GLM and sDNN performed best for narrow (>0.70), and broad (>0.65) quantiles, maintaining a decreasing sensitivity pattern across these categories and algorithms, although this trend was less pronounced. sDNN showed no statistically significant differences (adjusted p > 0.05) with GLM in Q3–Q4 and GBM in Q4 (Figure 6B).

In summary, sensitivity rankings remained stable across thresholds and niche breadth categories, while testing at the continental scale produced higher and less variable sensitivities than national-scale evaluation. Despite reduced absolute performance at the finer scale, the comparative hierarchy among algorithms was maintained.

Evaluation using the Boyce index highlighted a contrasting pattern. At the continental scale, GLM and MaxEnt achieved the highest median Boyce values (both ≈ 0.80), whereas sDNN and RF showed intermediate values (0.64, 0.56) and GBM and XGBoost exhibited lower ecological consistency (Figure 7A). At the national scale, GLM and MaxEnt ranked highest (>0.6), with all other algorithms showing median Boyce values below 0.5.

Regarding undefined Boyce values, GLM, sDNN, and GBM were particularly prone to this behavior for a substantial number of species (>50), reflecting flat or low-variance suitability surfaces that limit rank-based evaluation (Figure 7B). These undefined cases were more frequent at the continental scale, particularly for narrow quantile (Q1: 73 species) and decreasing toward broad quantile (Q4: 35 species), compared to the national scale (Q1: 55 species; Q4: 39 species). This indicates that Boyce-based diagnostics are primarily driven by low contrast predictions and limited discriminatory power, rather than by spatial aggregation effects.

To complete this analysis, Table 2 summarizes the mean sensitivity and Boyce Index values across algorithms and spatial scales. While the comparative hierarchy of algorithms remained consistent with the median-based evaluations shown previously, mean values were systematically lower than their corresponding medians across all algorithms and metrics. This divergence indicates a left-skewed distribution in model performance, reflecting a subset of species with markedly low predictive accuracy or ecological consistency that disproportionately reduced the overall averages. Nevertheless, algorithmic rankings remained stable, confirming that observed performance trade-offs are robust to the central tendency metric. Furthermore, sample composition for species with poor discrimination (continental sensitivity < 0.5) predominantly retained original occurrences, with cell-based samples virtually absent. This confirms that low accuracy is not an artifact of extreme environmental augmentation. However, XGBoost, sDNN and GBM included slightly higher neighbor-species counts than MaxEnt, GLM and RF for these specific taxa, suggesting a greater sensitivity to the structural variance introduced by neighbor-based sampling.

Effects of replay, architecture, and regularization in sequential training

Sequential training experiments revealed that architectural capacity and experience replay jointly controlled training stability and knowledge retention across species during large-scale sequential learning. Deeper network configurations consistently achieved stable performance, reaching median AUC values of 0.94 to 0.95 throughout sequential learning, whereas simplified architectures such as the 64–32–16–8–4 configuration exhibited rapid performance degradation indicative of catastrophic forgetting (Figure 8A). In these reduced architectures, experience replay partially mitigated forgetting but did not fully restore stable performance across species.

Learning rate played a critical role in maintaining stable sequential updates. Rates equal to or higher than 1.0 × 10 -5 preserved performance across the first 500 species, whereas lower values such as 1.0 × 10 -6 induced gradual forgetting, even under replay (Figure 8B).

Regularization strategies further modulated sequential training dynamics. L2 regularization stabilized early learning and reduced performance fluctuations, whereas dropout showed a largely invariant effect, comparable to the no-regularization control during sequential learning. Their combination improved early convergence and maintained stable performance across species, indicating complementary control over overfitting and gradient drift (Figure 8C).

Replay buffer size strongly affected training stability, with stable performance at small to intermediate capacities. Performance stabilized for buffer sizes ranging from 100 to 500 samples, whereas very large buffers exceeding 5000 samples (matching training sample size) reintroduced instability. This behavior is likely driven by the over-representation of earlier species in the replay memory and reduced adaptability to newly introduced species data. Notably, the unbuffered configuration remained stable and comparable to intermediate buffer sizes (Figure 8D).

Across training configurations evaluated on the test datasets using the best-performing settings (see Section 2.3.1), sequential configurations showed more stable sensitivity at both continental and national scales. Among these, the buffered setting achieved the highest median sensitivity (0.86 at the continental scale) and the lowest dispersion, while non-sequential training exhibited lower sensitivity and higher variability (Figure 9A). In contrast, Boyce Index values indicated a trade-off across training configurations, with non-sequential training achieving a higher median but exhibiting frequent negative outliers, whereas buffered and unbuffered sequential settings showed lower median values (0.4-0.5) and broader dispersion across species (Figure 9B). The undefined cases exceeded 100 species across all settings, with the non-sequential setting peaking at 185 species. The pattern wherein species within the narrow quantile exhibited a higher frequency of undefined cases (Q1: 77) than the broad quantile (Q4: 37) was again observed at the continental scale, while at the national scale, the trend persisted but was more pronounced in Q4 (Q1: 71; Q4: 51).

Extrapolation behavior and spatial response across niche breadth classes under historical and future climate conditions

Under the extrapolative threshold τ M E S S > − 10, ERI values showed limited differentiation among algorithms under both historical and future climate conditions. RF, GBM, XGBoost, and MaxEnt consistently achieved high ERI values above 0.95 in both periods, indicating low levels of extrapolation. sDNN showed slightly lower ERI values (0.91–0.94), while GLM clearly ranked lowest (0.71–0.76), separating from the remaining algorithms.

When applying the stricter extrapolation threshold τ M E S S > 0, differences among algorithms became more pronounced across both climate periods. RF achieved the highest ERI values (>0.89). XGBoost, MaxEnt, GBM, and sDNN formed an intermediate group with ERI values roughly between 0.71 and 0.87, whereas GLM clearly separated with low ERI values below 0.49. For most algorithms, ERI values were slightly lower under future climate projections, with a mean decrease of approximately 0.05 relative to historical conditions. Figure 10 presents ERI values across algorithms under historical and future climate conditions.

Across niche breadth quantiles, the spatial extent of EPF differed markedly among algorithms while largely following the ranking observed for ERI (Figure 11). EPF increased from narrow to broad quantiles, indicating modulation by niche breadth. GLM consistently exhibited the largest extrapolative areas across all EPF classes and quantiles, whereas RF showed minimal EPF. XGBoost, MaxEnt, and GBM displayed intermediate EPF patterns characterized by continuous and coherent regions, a pattern also observed for sDNN but with lower extrapolative extent than GLM. High EPF values were predominantly concentrated in tropical humid regions and climatic transition zones.

This ranking was consistent with correlations against the mean ensemble, which were highest for GBM and GLM (> 0.9) and decreased across MaxEnt, RF, sDNN, and XGBoost (0.65 – 0.83).

With respect to distribution change, niche breadth quantiles showed marked jumps in area, particularly between Q1 and Q4, where loss areas depleted persistence, compared to Q2 and Q3, where losses were substantial but less extensive and gains more pronounced (Figure 12). This pattern coincides with the marginal placement of Q1 and Q4 along the environmental gradient. Among algorithms, tree-based and boosting approaches (GBM, XGBoost, and MaxEnt), together with sDNN, exhibited loss-dominated patterns with localized gains, whereas RF showed a similar behavior but with losses consistently dominating across all quantiles. In contrast, GLM showed smoother transitions, characterized by extensive persistence and relatively larger gain areas, particularly across Q1–Q3.

Across quantiles, loss areas were predominantly concentrated in the Amazon Basin bioregion and its transition zones. The mean ensemble exhibited a similar non-linear behavior across niche breadth classes. Cohen’s kappa agreement for distribution change classes was highest for GBM and sDNN (>0.7), intermediate for MaxEnt and XGBoost (0.55–0.66), and lowest for RF and GLM (<0.34).

Following the distribution change analysis, the relative change index showed a pattern broadly comparable across niche breadth quantiles; however, responses from Q1 to Q3 remained largely similar, with a clear departure emerging in Q4 (Figure 13). In Q1–Q3, losses were concentrated across the Amazon Basin, with additional hotspots in the Chocó and Andean foothill regions, while gains were mainly observed in the Atlantic Forest.

In contrast, Q4 exhibited more spatially homogeneous and less intense losses extending across the Amazon Basin, including the Cerrado and Chaco regions. Algorithm-specific patterns showed that, except in Q4, RF, XGBoost, GBM, and sDNN exhibited coherent relative change signals dominated by losses, whereas MaxEnt displayed relatively more gain than loss areas, except in Q4. GLM was characterized by smoother spatial transitions, with gain areas prevailing.

The mean ensemble summarized these responses, exhibiting highest correlation for GBM and GLM (>0.86), similarly high values for MaxEnt and RF (~0.83–0.84), decreased for sDNN (0.76), and a sharp drop for XGBoost (0.19).

Bias and imbalance in large-scale SDMs and their implications for sampling and evaluation

Curation effects varied systematically along the niche breadth quantiles. Specifically, species in broad quantiles retained most of their original records, whereas those in narrow quantiles required additional sampling via an environmentally constrained feature-space augmentation approach due to structural data scarcity. Consequently, narrow quantiles (Q1 and Q2) relied more strongly on neighbor-based samples than broad quantiles (Q3 and Q4), creating a structural dependency on synthetic information that complicates direct algorithm comparison unless results are stratified or aggregated by niche breadth. Furthermore, these steps introduced a deliberate spatial alignment across both occurrence and background sampling, ensuring that background points shared the same spatial sampling structure as the occurrence data, thus mitigating sampling effort bias (Phillips et al., 2009). This strategy impacted SDM algorithms differently, exerting a stronger structural constraint on tree-based classifiers, which typically exhibit limited extrapolation capabilities (Heikkinen et al., 2012). Nevertheless, a persistent spatial density bias toward Ecuador, Colombia, Peru, and Brazil remained, reflecting unequal sampling effort across the region (Pinkert et al., 2022). Since species exhibiting limited continental-scale sensitivity predominantly relied on original occurrences (Table 2), this spatial artifact is likely a primary driver of lower predictive accuracy by restricting niche representation. While our approach provides a stable baseline for complex biodiversity datasets, alternative sampling methods could further refine predictions in unevenly surveyed regions.

This structural imbalance between narrow and broad quantiles, also conditioned performance evaluation. While sensitivity captures discrimination ability, it is largely insensitive to suitability contrast and spatial structure. Accordingly, low or undefined Boyce values reflect flat or weakly contrasted predictions that are common for data-scarce or environmentally restricted species, rather than model failure. For instance, shared representations in sDNNs stabilize discrimination but artificially smooth probability surfaces for narrow-range species, reducing the gradient contrast required for rank-based metrics. This behavior occurs because the Boyce Index primarily responds to suitability ranking and contrast rather than classification accuracy (Hirzel et al., 2006). Without explicit niche-breadth stratification and complementary metrics, ecological realism cannot be reliably assessed under niche-dependent sampling bias in large-scale SDMs, a limitation that was revealed by our results.

The use of independent continental- and national-scale test datasets confirmed that evaluation outcomes are conditioned as much by sampling and evaluation design as by algorithm choice. This reinforces that SDM performance at large scales is fundamentally shaped by data structure rather than algorithmic sophistication alone, a conclusion well-established in the SDM literature regarding the critical impact of sampling design on model evaluation (Ploton et al., 2020).

Performance trade-offs and structural implications of transfer learning for large-scale SDMs

In this study, transfer learning represented a major advantage of the proposed sDNN framework for multispecies SDMs. Shared representations allowed information to accumulate across taxa and improved prediction stability for species with restricted environmental tolerances. Conventional algorithms are trained independently for each species and show limited capacity to generalize across large species collections (Norberg et al., 2019). In contrast, the advantage of the proposed framework arises from training a single model sequentially across species. In our framework, effective transfer depended on sufficient architectural capacity, moderately constrained experience replay, and appropriately tuned learning rate and regularization parameters, as demonstrated by the sequential training experiments. Reduced architectures, excessively large replay buffers, or overly high learning rates consistently impaired adaptability. In some cases, these factors led to catastrophic forgetting, a behavior also reported in other domains of continual learning (Kirkpatrick et al., 2017). Together, these design choices ensured training stability and supported robust generalization across evaluation thresholds, niche breadth quantiles, and both continental and national scales.

Despite these advantages, sDNN performance revealed trade-offs inherent to shared representations. While discrimination performance remained consistently high, predictions exhibited smoother suitability surfaces, particularly within narrow quantiles, alongside an increase in extrapolation toward broader quantiles. This behavior reflects information sharing across species that stabilizes predictions but yields smoother suitability responses. Consequently, this mechanism can reduce niche contrast, especially for species with restricted environmental tolerances, rather than indicating model overfitting (Merow et al., 2014). This interpretation is supported by stable performance under independent continental- and national-scale evaluation, strong correlation with the mean ensemble, and the absence of fragmented or extreme suitability patterns typically associated with overfitted SDMs (Araújo and New, 2007; Radosavljevic and Anderson, 2014). Future work should explore bioregion-level training strategies to strengthen species-specific responses for rare taxa and better leverage shared representations.

From an applied perspective, the proposed sDNN framework is not intended for detailed single-species inference. Instead, it is better suited for large-scale multispecies SDMs characterized by heterogeneous and imbalanced data. Although evaluated here using food-relevant plant species, the framework provides a methodological basis that could be extended to other taxonomic groups under similar data constraints. Such extensions would require explicit consideration of taxon-specific processes, including mobility, dispersal, and life-history traits (Franklin, 2013), as these are expected to influence sampling bias, calibration areas, and extrapolation behavior. Under these conditions, when predictive stability and cross-scale consistency are prioritized, this approach provides clear advantages. However, it cannot replace classical SDM alternatives, which remain better suited for detailed model interpretation and capturing fine-scale niche contrasts.

Furthermore, unlike complex DNN architectures that demand high-performance computing clusters, the proposed sDNN is highly computationally efficient and can be executed on standard workstations. A comparative assessment of computational resource consumption demonstrated that the sDNN maintains stable execution times during the fine-tuning and inference phases, remaining highly competitive with benchmark algorithms without requiring specialized infrastructure (see Supplementary Information, Section 2). However, from a computational perspective, sequential training is time-consuming and restricts parallel execution, resulting in lower GPU utilization compared with fully independent model training (Xu et al., 2025). Future work should address these constraints by exploring layer and sequence parallelization strategies (Lim et al., 2024) to improve scalability on modern GPU architectures, while also systematically evaluating the impact of species presentation order, as the single randomized sequence used here highlights order dependency as an unexplored limitation.

Extrapolation and robustness of multispecies SDMs across niche breadth gradients under climate change

Our results indicate that extrapolation is a central constraint in large-scale multispecies SDMs, with consistent differences among algorithms and modulation by niche breadth. Algorithms achieving higher sensitivity, such as GBM and sDNN, exhibited moderate extrapolative areas. Conversely, more conservative approaches, particularly RF, restricted predictions to environmentally reliable space, contrasting with the extensive extrapolation observed in GLMs. This pattern highlights systematic differences in extrapolative behavior across algorithms, which are further structured by species niche breadth. Species with broad environmental tolerances occupy larger extrapolative areas, while those with restricted tolerances, despite more limited extrapolation, produce less reliable projections under environmental novelty. This behavior is consistent with the effects of incomplete niche representation on response curve behavior as described by Thuiller et al. (2004). It is further reinforced here by the already established greater reliance of taxa with restricted tolerances on augmented samples within this framework.

Under future climate conditions, increased environmental novelty did not induce broad instability across algorithms, as most methods preserved comparable spatial change patterns between historical and future periods. Nevertheless, departures from this consistency were mainly associated with species with restricted environmental tolerances. In this context, the sDNN exhibited higher spatial coherence, reflected by the highest Cohen’s kappa agreement (>0.7) with the ensemble mean. Conversely, the GLM exhibited the most anomalous behavior across niche breadth classes, characterized by the largest extrapolative areas (highest EPF) and the lowest consensus with the ensemble mean (Cohen’s kappa <0.34), resulting in gain-dominated patterns that diverged from the loss-dominated signals characteristic of the remaining algorithms. This contrast reflects not only the advantages of sequential training discussed in Section 4.2, but also the greater flexibility of sDNN and related algorithms in modeling non-linear species-environment relationships compared with parametric approaches such as GLMs (Elith et al., 2006). These limitations become particularly relevant when using PCA-transformed predictors, where linear responses can lead to diffuse or unstable projections under environmental novelty (Dormann et al., 2013). While these projections are inevitably conditioned by the underlying data augmentation required for data-scarce taxa, the overall consistency across algorithms remains robust. Consequently, future work in this context should primarily focus on assessing how alternative dimensionality reduction strategies affect extrapolation and robustness in large-scale multispecies SDMs under climate change.

Despite these methodological constraints, a consistent non-linear pattern of distributional change emerged across the niche breadth gradient, largely independent of algorithmic choice. This pattern suggests that climate-driven redistribution does not scale monotonically with niche breadth, but instead concentrates at the extremes of the gradient, corresponding to taxa with highly restricted and broad environmental tolerances. The vulnerability of species with restricted tolerances to climate change has been widely reported (Slatyer et al., 2013; Thuiller et al., 2005). In contrast, evidence for climate-driven suitability contraction in broad-tolerance taxa is more contradictory and often depends on how niche breadth is defined, the taxonomic scope, and the spatial context (Sweeney and Jarzyna, 2022). Within this context, niche breadth in this study was derived from PCA-transformed predictors integrating climatic and land-cover gradients. This operationalization frames niche breadth as an emergent property of the predictor structure and its combination, emphasizing context dependence and constraints under environmental novelty. Consequently, it avoids assuming any inherent advantage for species with broad tolerances. This implies that the way niche breadth is defined and operationalized can materially influence ecological inference in large-scale SDMs, underscoring the need for careful definition and standardization.

Furthermore, it is important to acknowledge that spatial shifts in food plants represent changes in baseline resource availability rather than definitive biological outcomes for dependent species. Species distribution models fundamentally project abiotic habitat suitability and do not incorporate complex biotic interactions, demographic dynamics, or physiological adaptability (Araújo et al., 2019). Consequently, inferring direct biological contexts solely from these spatial differences has inherent limitations. Nevertheless, when evaluating hundreds of species lacking detailed mechanistic data, quantifying shifts in potential food plant distribution provides a crucial first-order approximation of landscape vulnerability and resource bottlenecks. Although the predictor set used in this study captures major climatic, topographic, hydrological, and land-cover gradients, other ecological processes may also influence the comprehended distributions of food-relevant plant species and contribute to prediction uncertainty. Local biotic interactions, including competition, mutualisms, herbivory, and pathogen pressure, together with dispersal limitations, can modify species occurrences independently of environmental suitability (Yates et al., 2018). Likewise, ecological responses may depend on temporally lagged environmental conditions rather than concurrent climatic variables, particularly for processes such as recruitment, productivity, and phenology that integrate environmental conditions over preceding weeks or months (Chen and Peñuelas, 2025). In addition, differences in occurrence data completeness, taxonomic accuracy, and spatial sampling intensity remain important sources of uncertainty in continental biodiversity datasets. Beyond spatial coverage, sampling bias and predictive performance are strongly influenced by data quantity, such as sample size, and data quality, including species detectability. These constraints particularly affect rare or range-restricted species, where complex interactions between local abundance and occurrence probability depend heavily on sampling effort and detection accuracy (Bouchet et al., 2019; Wang and Jackson, 2023). Future multispecies SDMs could therefore benefit from incorporating temporally explicit predictors, biotic information, detection-aware frameworks, and process-informed modeling approaches to better represent ecological mechanisms underlying species distributions (Urban et al., 2016). Similarly, projections in this study were based on a single global circulation model and a high-emission scenario. These findings illustrate scenario-specific algorithmic behavior under extreme environmental shifts rather than providing absolute forecasts of species distributions. Future research should incorporate multi-model ensembles and broader emission pathways to systematically quantify climatic uncertainty. This enhanced analysis will improve predictive accuracy and deliver a clearer understanding of climate change implications, transitioning from algorithmic benchmarking to actionable ecological forecasting.

This study addressed the challenge of modeling food-relevant plant species distributions at large spatial scales under conditions of data scarcity, sampling bias, and environmental heterogeneity. By integrating environmentally informed sampling, extrapolation diagnostics, and sDNN training with experience replay, the proposed framework advances large-scale SDM workflows beyond architecture-centric solutions and reliance on massive computing resources. Across independent continental- and national-scale test datasets, the framework showed stable predictive performance across species with contrasting niche breadths. The results demonstrated improved sensitivity for taxa with restricted environmental tolerances and robustness across binarization thresholds and spatial scales. While this performance gain reflects the high plasticity of sDNNs, it underscores that scalable multispecies SDMs require a balanced integration of predictive accuracy, ecological realism, and computational feasibility.

Specifically, our evaluation revealed key structural trade-offs among the methodologies. Higher discrimination was associated with smoother suitability surfaces and increased extrapolation for species with broader tolerances, whereas more conservative algorithms restricted predictions to environmentally reliable space. Extrapolation analyses indicated that algorithmic behavior under environmental novelty was driven primarily by niche breadth and sampling design rather than by algorithm choice alone. Despite moderate extrapolation compared to parametric methods, the sDNN retained coherent spatial responses and stable rankings across historical and future climate conditions, consistent with controlled generalization. Ultimately, these findings reinforce that robust large-scale modeling depends fundamentally on managing data structure and training design rather than relying solely on algorithmic sophistication.

The proposed framework establishes a methodological baseline for sDNN-based architectures and the development of large-scale SDMs under climate change, with demonstrated applicability to food-relevant plant species. Its design is not taxon-specific and can potentially be extended across taxonomic groups, contributing to spatially explicit biodiversity assessments under heterogeneous data availability and environmental constraints.