Winter warming and snowpack loss synergistically drive amphibian and reptile extirpation: a deep learning attribution study of climatic mechanisms

Abstract

Winter climate change is an underrecognized but potentially strong predictor of biodiversity loss in temperate ecosystems. This study investigates the role of winter warming and snowpack decline in driving local extirpation of amphibians and reptiles across the Midwestern United States, providing model-based insights consistent with underlying ecological mechanisms. A multi-decadal, species-integrated dataset was compiled, combining population dynamics, survival metrics, and climate variables for four climate-sensitive taxa: Acris blanchardi, Sistrurus catenatus, Clonophis kirtlandii, and Regina septemvittata. A fully connected deep learning model was developed to predict population persistence, survival probability, and extirpation risk based on winter temperature, snowpack, and precipitation. Model interpretability was achieved using SHapley Additive exPlanations (SHAP) to quantify feature importance and interaction effects. The model demonstrated strong predictive performance, achieving RMSE values of 0.08–0.15 and R² values of 0.71–0.87 for continuous outputs, alongside AUC-ROC scores of 0.84–0.93 and F1-scores of 0.79–0.88 for extirpation classification. Feature attribution revealed that snowpack variables were the dominant predictors, contributing 42–57% of total model importance, followed by temperature (28–39%) and precipitation (12–21%). Critically, interaction effects between snowpack loss and temperature anomalies amplified extirpation risk by 35–62%, with interaction terms contributing an additional 18–26% to model output. Nonlinear threshold behavior was identified, with extirpation risk increasing rapidly when snowpack declined below 30–40% of historical norms and temperature anomalies exceeded +2 °C. Species-specific analyses indicated that Clonophis kirtlandii exhibited the highest sensitivity, with extirpation probability increasing by up to 68% under combined stress conditions, while Acris blanchardi showed a 6–11% increase in mortality probability per 10% reduction in snow depth. Across taxa, long-term population declines ranged from 35–65%, accompanied by a 40–70% reduction in occupied habitats and a 20–35% increase in interannual variability post-2000. These findings demonstrate that winter climate dynamics particularly the loss of subnivean insulation and increased thermal variability drive extirpation through interacting thermal and hydrological pathways. By integrating deep learning with interpretable AI, this study provides quantitative evidence that compound winter stressors produce disproportionately large ecological impacts, highlighting the need to incorporate winter processes into climate risk assessments and conservation strategies.

1 Introduction

Ecological change in the Anthropocene is increasingly characterized by the emergence of novel and disappearing climatic regimes (Williams et al., 2007), driving large-scale biodiversity redistribution (Pecl et al., 2017). These transformations are inherently nonlinear, with ecosystems exhibiting abrupt transitions, tipping points, and regime shifts (Scheffer et al., 2001; Alkhayuon et al., 2022). Global assessments indicate that such climatic pressures are accelerating, contributing to rising extinction risk (Urban, 2015; IPCC, 2021), while broad syntheses confirm that ecological responses to climate change are already widespread and governed by interacting, threshold-dependent processes (Walther et al., 2002).

Recent quantitative analyses further highlight the accelerating nature of global warming and its ecological consequences. Consistent upward temperature trends across multiple datasets demonstrate intensifying warming patterns (Bukaita, 2024), while index-based approaches reveal that temperature variability—rather than mean warming alone plays a critical role in shaping environmental extremes (Bukaita et al., 2024). This variability exhibits strong seasonal asymmetry, with winter temperatures increasing more rapidly than annual averages. Such winter amplification reduces snow cover duration, alters freeze–thaw dynamics, and disrupts soil thermal regimes, thereby intensifying ecological instability during dormant periods.

Within this context, snowpack decline emerges as a central mechanism linking climate change to ecological disruption. Snowpack acts as a critical thermal buffer, stabilizing soil and near-surface temperatures during winter. Its reduction contributes to compound climatic extremes (Liu et al., 2024) and increases exposure of overwintering organisms to subfreezing conditions despite overall warming trends (Johnston et al., 2021). Concurrently, increasing temperature variability amplifies the frequency and intensity of freeze–thaw cycles (Bukaita et al., 2024), imposing physiological stress that often exceeds that caused by gradual warming. These effects are particularly pronounced in ectothermic organisms such as amphibians and reptiles, whose survival depends on stable thermal environments. Repeated freeze–thaw events can induce cellular damage, disrupt osmotic balance, and impair energy regulation, increasing mortality risk during overwintering.

The ecological consequences of altered winter conditions extend across multiple levels of biological organization. Changes in snow regimes restructure trophic interactions (Brodie et al., 2014; Ren et al., 2017), disrupt phenological synchronization across species (Post and Stenseth, 1999; Parmesan, 2006), and influence hydrological processes that determine habitat persistence (Ratsch et al., 2020). Even gradual climatic shifts can accumulate to exceed adaptive capacity, triggering rapid ecological transitions or collapse (Vanselow et al., 2019). These patterns are further intensified by the interaction of multiple stressors, as climate-driven extinction risk is often shaped by synergistic effects rather than single drivers (Brook et al., 2008; Cahill et al., 2013; Velasco et al., 2021; Duffy et al., 2022).

Empirical studies of climate-sensitive species illustrate how these interacting mechanisms manifest in real populations. Amphibian and reptile populations exhibit strong sensitivity to environmental variability, with documented instability, demographic shifts, and regional declines (United States Geological Survey (USGS), [[NoYear]]; Lehtinen et al., 2025). Vulnerability is further reinforced by demographic sensitivity and extinction risk under varying environmental and management conditions (King et al., 2012; Jones et al., 2017; Hileman et al., 2018), as well as increasing genetic fragmentation observed across populations (Baker et al., 2018; United States Geological Survey (USGS), 2023), which informs conservation and recovery strategies (Government of Canada, [[NoYear]]). Rare and habitat-specialist species exhibit additional constraints, including limited detectability and fragmented distributions (Stewart et al., 2023; Crawford, 2024; Illinois Natural Heritage Program, 2024). Broader evidence across herpetofauna indicates that populations are often spatially constrained and fragmented, reinforcing the role of climate-sensitive, habitat-dependent population structure (Stanley, 2007; Leuenberger et al., 2019; Ibach, 2022; Bukaita et al., 2024).

These ecological responses are closely tied to physiological constraints governing overwinter survival. Amphibians and reptiles depend on stable thermal and hydric conditions, with snowpack and soil temperature regimes playing a critical role in regulating energy balance, freeze tolerance, and survival thresholds (Muths et al., 2020). Access to thermally buffered microhabitats is essential for maintaining physiological stability during dormancy or brumation (Leuenberger et al., 2019); however, increasing temperature variability disrupts these conditions, exposing organisms to repeated freeze–thaw cycles that can exceed tolerance limits (Bukaita et al., 2024). Disease dynamics further interact with these stressors. Although pathogen activity may be reduced in cold conditions, winter-induced physiological stress can weaken immune function and increase susceptibility, particularly for pathogens capable of persisting at low temperatures or reactivating during intermittent warming events (Mitchell and Tully, 2005; McKenzie et al., 2021).

Despite extensive evidence linking climate change to biodiversity loss, a critical gap remains in attributing the relative contributions of interacting climatic drivers. Existing studies demonstrate strong relationships between climate variables and species occupancy but are largely based on correlative frameworks (United States Geological Survey (USGS), 2018; Weiskopf et al., 2022). While recent quantitative approaches improve characterization of temperature trends and variability (Bukaita, 2024; Bukaita et al., 2024), they do not explicitly resolve how interacting winter-specific processes such as snowpack loss, temperature variability, and short-duration extreme events translate into biological outcomes.

Addressing this limitation requires analytical frameworks capable of capturing nonlinear, high-dimensional interactions. The present study adopts a deep learning attribution approach to decompose complex climatic influences into interpretable contributions. By quantifying the relative importance and interaction of snow depth, temperature, and precipitation, this framework moves beyond correlation toward mechanistic understanding of climate-driven extirpation. Importantly, it captures threshold effects and interdependencies among climatic variables under winter-specific conditions, providing a more ecologically realistic basis for interpreting species decline and improving climate-driven risk assessment.

2 Methods

2.1 Dataset compilation

A rigorous dataset was constructed to encompass all documented extinct and locally extirpated species, with a primary analytical focus on amphibian and reptile taxa. This compilation, visually summarized in Figure 1 (Bukaita and Ghiurau, 2025), synthesizes long-term population monitoring data, systematic field surveys, and species-specific ecological investigations. The dataset includes winter climate variables, such as temperature and snowpack, selected to quantify their contributions to population outcomes in the deep learning framework. This analysis focuses on four climate-sensitive focal species: the Blanchard’s Cricket Frog (Acris blanchardi), Eastern Massasauga (Sistrurus catenatus), Kirtland’s Snake (Clonophis kirtlandii), and Queen Snake (Regina septemvittata).

Figure 1

2.1.1 Species population data acquisition

High-resolution population data for the four focal species were compiled through an integrated and explicitly documented data acquisition framework combining peer-reviewed literature, long-term government monitoring programs, and targeted field surveys. Primary sources included datasets from the U.S. Geological Survey (USGS) Amphibian Research and Monitoring Initiative, the Michigan Natural Features Inventory, and additional state and regional conservation databases (Government of Canada, [[NoYear]]; United States Geological Survey (USGS), [[NoYear]]; King et al., 2012; Jones et al., 2017; Baker et al., 2018; Hileman et al., 2018; Leuenberger et al., 2019; Stewart et al., 2023; United States Geological Survey (USGS), 2023; Crawford, 2024; Illinois Natural Heritage Program, 2024; Lehtinen et al., 2025). Published ecological studies and species-specific monitoring reports were systematically screened and incorporated to supplement gaps in spatial or temporal coverage.

To ensure methodological consistency, all datasets were harmonized through a standardized preprocessing pipeline. This included aligning spatial coordinates to known species habitats, reconciling temporal scales across survey efforts, and normalizing measurement units for demographic variables. Spatiotemporal patterns of occupancy and metapopulation dynamics were reconstructed using long-term survey records spanning wetland complexes, prairie remnants, and riparian systems. Estimates of abundance and survival were parameterized from mark–recapture studies, supported by controlled mesocosm experiments and artificial pond array observations where available.

Reproductive output and population structure were further characterized using reported vital statistics, including egg and clutch counts, operational sex ratios, and established proxies for recruitment success. Data inclusion criteria required sufficient methodological transparency and temporal resolution to support integration, while duplicate records across sources were identified and removed to avoid sampling bias. This structured approach ensured traceability of all inputs and improved reproducibility of the compiled population dataset.

2.1.2 Development of the extinct and extirpated species shortlist

To provide necessary ecological context regarding climatic impacts on North America’s broader fauna, a deliberate shortlist of extinct and locally extirpated species was generated. This compilation specifically targets species whose decline or disappearance is linked in the literature to winter climate variables (Table 1). The table serves as a reference point, detailing the final year of observation and the primary climatic or secondary ecological predictors contributing to local population collapse.

2.1.3 Integration of climate and snowpack metrics

To quantify winter stressors relevant to each focal species, a standardized suite of climatic datasets was integrated with explicit citation of all data sources. Temperature records, including daily minimum and maximum values from 1980 to 2025 (National Centers for Environmental Information (NCEI), [[NoYear]]), were obtained from the National Oceanic and Atmospheric Administration’s National Centers for Environmental Information (NOAA/NCEI), with appropriate dataset-level citations included in the reference list following NOAA data usage guidelines.

Snowpack dynamics, including Snow Water Equivalent (SWE) and snow cover duration, were derived from the National Snow and Ice Data Center (NSIDC) and NASA’s Moderate Resolution Imaging Spectroradiometer (MODIS). These datasets are cited according to provider-specific requirements, including dataset identifiers, access dates, and version numbers, in accordance with NSIDC DAAC citation guidelines. Spatial interpolation was applied to align gridded snowpack data with species-specific habitat locations.

Precipitation data, including rainfall and snowfall totals, were compiled from hydrologic monitoring programs and validated gridded climate products, each referenced with their corresponding dataset citations and metadata sources. All climatic variables were temporally synchronized with biological observations to ensure consistency in seasonal alignment, enabling accurate representation of winter-driven mortality and reproductive impacts while maintaining full traceability and reproducibility of the data sources used.

2.1.4 Strategic purpose of the integrated dataset

The development of this structured, integrated dataset is fundamental to achieving the following analytical goals:

• Risk Assessment: To facilitate comparative modeling of climatic predictors affecting focal species relative to historically inferred extirpation drivers.

• Comparative Ecological Analysis: To facilitate comparative modeling of climatic predictors affecting focal species relative to historically inferred extirpation drivers.

• Advanced Statistical Modeling: To provide necessary input data for deep learning algorithms, enabling the subsequent use of SHAP-based (SHapley Additive exPlanations) feature attribution to disentangle the discrete and interactive contributions of temperature, snowpack, and precipitation to observed local population collapses.

This structured synthesis ensures that all subsequent modeling and data-driven inference of extirpation predictors are firmly grounded in robust, empirical observations of both population demography and regional climate variables.

2.2 Spatial and landscape analysis

Spatial analyses were conducted using R version 4.3.1 (R Core Team, 2023), specifically leveraging the ‘landscapemetrics’ package to quantify both patch-level and landscape-level metrics. Landscape changes were quantified using multi-decadal data from the National Land Cover Database (NLCD) and supplemental satellite imagery spanning the period from 1900 to 2025. Within this framework, habitat patch reduction was defined as the percentage change in total suitable habitat area over time, while fragmentation was rigorously assessed through Patch Density (PD) and Mean Patch Size (MPS). To ensure the integrity of area-based calculations across the Midwestern United States study region, all spatial datasets were projected to the Albers Equal Area Conic coordinate system.

2.3 Deep learning model design

A fully connected neural network was developed to predict species population responses from winter climate variables. The input feature set consisted of snow depth, temperature anomalies, and winter precipitation, while the target outputs included population persistence, survival probability, and local extirpation events represented through both continuous and binary metrics.

The model architecture comprised three hidden layers with 128, 64, and 32 neurons, respectively, using ReLU activation and a dropout rate of 0.2 to reduce overfitting (Figure 2). Model evaluation was based on an 80/20 train–test split, and all reported performance metrics were additionally validated using 5-fold cross-validation to ensure robustness and reduce sampling bias across partitions. Training was conducted using mean squared error for continuous outputs and binary cross-entropy for classification tasks. Optimization was performed using the Adam optimizer with a learning rate of 0.001 over 500 epochs, with early stopping applied based on validation loss. All models were implemented in Python 3.11 using TensorFlow 2.12, mapping a three-dimensional climate feature vector x R³ to a corresponding target vector y R³ through nonlinear transformations.

Figure 2

Where represents the Rectified Linear Unit (ReLU) activation function.

The mean squared error (MSE) loss for continuous survival probability is defined as:

Binary Cross-Entropy for extirpation events:

The Adam optimizer (learning rate = 0.001) iteratively updates weights W to minimize the combined loss over 500 epochs, with early stopping used to improve generalization on unseen winter climate data.

Figure 2 illustrates the complete model architecture and training pipeline. The input layer processes the three winter climate predictors, followed by three fully connected hidden layers with progressively reduced dimensionality. Each hidden layer applies ReLU activation and dropout regularization (0.2). The output layer produces simultaneous predictions for population persistence and survival probability (continuous outputs) and local extirpation events (binary output). The training framework integrates multi-loss optimization (MSE and binary cross-entropy), with model validation performed using both an 80/20 split and 5-fold cross-validation to ensure robustness and stability of results.

2.4 Baseline model comparison and justification of deep learning framework

To rigorously evaluate the effectiveness of the deep learning approach, a suite of widely used baseline models in ecological and environmental modeling was implemented, including Generalized Linear Models (GLM), Generalized Additive Models (GAM), Random Forest (RF), and Extreme Gradient Boosting (XGBoost). These models represent complementary statistical and machine learning frameworks for capturing linear, nonlinear, and interaction effects in structured ecological data.

GLMs were implemented using Gaussian link functions for continuous outcomes (population persistence and survival probability) and binomial logistic regression for extirpation events. GAMs extended this formulation by incorporating spline-based smooth functions to capture nonlinear relationships between winter climate variables and population responses.

Ensemble learning methods were configured with explicitly defined hyperparameters to ensure reproducibility and fair comparison. The Random Forest model was constructed with 500 decision trees, a maximum tree depth of 20, and the square-root rule for feature selection at each split (max features = “sqrt”), with bootstrap sampling enabled. The XGBoost model was implemented with 500 boosting rounds, a maximum tree depth of 6, a learning rate of 0.05, subsampling rate of 0.8, and L2 regularization (λ = 1.0) to reduce overfitting while maintaining predictive flexibility. These configurations were selected based on commonly recommended defaults for structured ecological datasets and preliminary tuning experiments.

Recent approaches such as graph-based label dependency models and transformer architectures, including Vision Transformers, have demonstrated promise in ecological applications; however, their effectiveness remains limited in settings with moderate sample sizes and non-spatiotemporally explicit structured inputs such as those used here.

All models were trained using the same input feature set consisting of snowpack metrics (depth and duration), temperature anomalies, and winter precipitation. Target variables included continuous measures of population persistence and survival probability, as well as binary extirpation indicators. Identical preprocessing steps, including normalization and missing data handling, were applied across all models to ensure comparability.

Model performance was evaluated using Root Mean Square Error (RMSE) and coefficient of determination (R²) for regression tasks, and Area Under the Receiver Operating Characteristic Curve (AUC-ROC) and F1-score for classification tasks. These metrics provide complementary evaluation of predictive accuracy, model calibration, and classification robustness under imbalanced ecological conditions.

To ensure methodological consistency, all models were evaluated under the same train–test partitioning and cross-validation framework, ensuring that performance differences reflect model capacity rather than data exposure bias or sampling variability.

2.5 Feature attribution using SHAP

To quantify the relative contribution of each climate variable to population outcomes, SHapley Additive exPlanations (SHAP) were applied to the trained deep learning models (Figure 3). SHAP values were computed for each input feature to obtain feature-specific importance scores reflecting their contribution to model predictions. Snow depth, temperature anomalies, and precipitation were evaluated both individually and jointly to identify potential interaction effects and synergistic contributions to extirpation risk. In addition, global SHAP summary analyses were generated to characterize species-specific sensitivity patterns across the dataset, providing quantitative, model-based evidence that aligns with hypothesized ecological mechanisms governing winter-driven population responses.

Figure 3

2.6 Methodological framework for spatiotemporal integration of heterogeneous population time series

To support rigorous modeling of winter climatic predictors and their associated risks of amphibian and squamate extirpation, a comprehensive, multi-source population database was synthesized. This framework prioritized the standardization of disparate longitudinal data metrics such as raw abundance counts, occupancy probabilities, fecundity estimates, and survival rates. The core objective was to create a coherent, normalized input layer suitable for deep learning model training and SHAP-based feature attribution analysis.

2.6.1 Data acquisition and species-specific consolidation

Population datasets for the four focal species across the North American Midwest and Great Lakes regions were compiled using an integrated framework combining peer-reviewed literature, long-term governmental monitoring programs (including the U.S. Geological Survey USGS monitoring datasets (United States Geological Survey (USGS), [[NoYear]]; United States Geological Survey (USGS), 2023) and the Michigan Natural Features Inventory [MNFI]), and targeted field surveys reported in previous ecological studies (United States Geological Survey (USGS), [[NoYear]]; Baker et al., 2018; United States Geological Survey (USGS), 2023; Lehtinen et al., 2025). This multi-source approach ensures consistency across ecological scales while capturing both demographic and spatial variability relevant to winter-driven population dynamics.

2.6.1.1 Blanchard’s cricket frog (Acris blanchardi)

For Acris blanchardi, (Figure 4) population dynamics were reconstructed using long-term occupancy and metapopulation datasets from Upper Midwest wetland systems (WI, MN, MI, IL, IA) derived from regional monitoring programs and published survey syntheses (United States Geological Survey (USGS), [[NoYear]]; Lehtinen et al., 2025) as summarized in Table 2. Additional abundance and survival estimates were obtained from controlled mesocosm experiments and artificial pond colonization studies in Kansas (McCallum and Trauth, 2021), while reproductive and demographic structure were parameterized using field-based surveys reporting fecundity and sex ratio variation (Hoskins and Boone, 2017; Hoskins and Boone, 2018).

Figure 4

Experimental mesocosm studies contributed survival rate estimates across geographically distinct populations, serving as proxies for population viability under environmental stressors. Regional population trends were further informed by systematic surveys of lake and wetland systems, documenting declines and local extirpations over time. Finally, long-term ecological observations of reproductive output and sex ratios were incorporated to reflect key demographic processes influencing population stability.

The Blanchard’s Cricket Frog (Acris blanchardi) serves as one of the primary focal species for assessing the impact of three key climatic variables on amphibian persistence:

• Snow Depth: Loss of snowpack reduces critical thermal insulation, exposing the frog to lethal soil temperature fluctuations during winter.

• Temperature: Increasing winter temperatures accelerate metabolic rates in dormant individuals, leading to the premature depletion of lipid reserves (metabolic exhaustion).

• Precipitation: Shifts in winter precipitation from snow to rain lead to altered wetland hydroperiods, causing the desiccation of essential breeding habitats.

2.6.1.2 Eastern massasauga (Sistrurus catenatus)

Population data for Sistrurus catenatus were synthesized from multiple complementary sources as shown in Table 3 to characterize demographic structure and population trends across its geographic range (Figure 5). Long-term monitoring datasets from Ontario (Canada), Michigan, and the broader U.S. Midwest provided foundational information on survival rates, population persistence, and temporal variation in abundance (Government of Canada, [[NoYear]]; King et al., 2012). These datasets include standardized mark–recapture studies conducted across multiple sites, enabling robust estimation of adult survival and site-level population size (King et al., 2012; Jones et al., 2017).

Figure 5

Range-wide demographic information was supplemented using integrated survival datasets covering multiple populations across the species’ distribution, which collectively capture spatial variation in mortality and recruitment patterns (King et al., 2012). Genetic sampling data from Bois Blanc Island (Michigan) were incorporated as a proxy for effective population size (Ne), providing additional insight into genetic structure and long-term viability (United States Geological Survey (USGS), 2023). Regional conservation assessments and population summaries from Environment and Climate Change Canada further informed estimates of population status and distribution trends in the northern portion of the range (Government of Canada, [[NoYear]]).

Together, these datasets provide a multi-scalar representation of Sistrurus catenatus population dynamics, integrating demographic, genetic, and spatial components relevant to assessing climate sensitivity and extinction risk.

The integrated dataset enables direct coupling with key winter climate variables: snow depth, temperature variability, and precipitation regime. This allows the deep learning framework to detect:

• Snowpack-dependent survival thresholds: Reduced snow depth removes subnivean insulation, exposing hibernation burrows to lethal temperature extremes and increasing overwinter mortality.

• Temperature-driven effects: Elevated winter temperature variability induces freeze–thaw cycles and metabolic stress, exceeding physiological tolerance limits and accelerating energy depletion.

• Hydrological and habitat-mediated effects: Shifts from snow to rain alter wetland hydrology and soil conditions, destabilizing hibernacula and indirectly increasing mortality risk.

2.6.1.3 Kirtland’s snake (Clonophis kirtlandii)

Population data for Clonophis kirtlandii were compiled from field surveys, occupancy studies, environmental DNA (eDNA) monitoring, and range-wide distribution datasets to characterize abundance, detectability, and spatial fragmentation (Figure 6).

Figure 6

Field-based survey data from Illinois prairie–wetland systems (2019–2021) documented 77 individuals across 226 surveys, indicating extremely low detection probability even under intensive sampling effort (Stewart et al., 2023). Occupancy assessments from the same region identified only three consistently occupied sites, highlighting restricted local persistence (Illinois Natural Heritage Program, 2024) as shown in Table 4.

At broader spatial scales, historical and contemporary records indicate continued population decline across the Midwestern United States, with multiple documented local extirpations and contraction of the species’ historical range since the 1970s (Crawford, 2024). Environmental DNA surveys conducted in Indiana wetlands (2020) further support this pattern, with only a single detection across 380 samples, emphasizing both rarity and detection difficulty (Ratsch et al., 2020).

Range-wide distribution datasets from the U.S. Geological Survey USGS monitoring datasets (United States Geological Survey (USGS), [[NoYear]]; United States Geological Survey (USGS), 2023) indicate a fragmented population structure confined to isolated sub-watersheds with limited connectivity among populations (United States Geological Survey (USGS), 2018). Regional metapopulation assessments across Illinois, Indiana, and Ohio further confirm that extant populations are small, spatially isolated, and highly vulnerable to stochastic environmental variation (United States Geological Survey (USGS), 2018; Illinois Natural Heritage Program, 2024).

2.6.1.4 Queen snake (Regina septemvittata)

Population data for Regina septemvittata were synthesized from capture–recapture studies, telemetry monitoring, localized abundance surveys, and broader ecological assessments to characterize density, survival, and spatial distribution across its range (Figure 7).

Figure 7

High-resolution capture–recapture data from central Kentucky stream systems (2016) provided baseline density estimates ranging from 6 ± 5 to 63 ± 10 individuals per kilometer of stream habitat (Leuenberger et al., 2019). These data were derived from standardized field sampling of 119 individuals across multiple stream reaches as given in Table 5.

Multi-site telemetry studies conducted across the eastern United States (2018–2020) indicated short-term survival stability, with no observed mortality events during the monitoring period; however, these estimates reflect short temporal windows and do not represent long-term demographic trends (Leuenberger et al., 2019).

Localized sampling in Pennsylvania (early 2000s) recorded extremely low detection rates (n = 2 individuals), suggesting strong regional variation and possible local population decline or extirpation (Stanley, 2007). Range-wide ecological assessments further indicate that Regina septemvittata populations are highly patchy and exist as discrete stream-associated clusters rather than continuous populations across North America (Stanley, 2007; Bukaita et al., 2024).

An adult Queen Snake (Regina septemvittata) consuming a soft-shell crayfish along a rocky stream bank illustrates its strong dependence on aquatic prey and stable stream environments. In this study, the species serves as a focal organism for evaluating the influence of three key winter climate variables on persistence. Deep learning attribution identifies the Queen Snake as highly sensitive to the interaction of these predictors:

• Snow Depth: Reduced snowpack diminishes subnivean insulation, exposing overwintering habitats to rapid freeze–thaw cycles that can induce lethal physiological stress.

• Temperature: Rising winter temperatures increase thermal variability and create trophic mismatches, whereby crayfish (the primary prey) emerge earlier than the snakes exit hibernation, disrupting critical spring foraging.

• Precipitation: Shifts from snow to rain alter stream hydrology and sediment conditions, affecting both overwintering refugia and prey availability.

Together, these interacting climatic factors form a two-pronged mechanism linking thermal stress and trophic disruption to increased local extinction risk in R. septemvittata.

2.6.2 Analytical strategy for data harmonization and integration

To convert heterogeneous ecological and climatic datasets into a unified structure suitable for neural network training, a systematic data harmonization framework was applied. All time series were first standardized to a common temporal resolution, either annual or seasonal, enabling consistent comparison across species and study sites. Abundance-related measures, including population counts, egg production, and recruitment indices, were transformed into dimensionless relative population indices to ensure compatibility with survival and occupancy-based datasets. These standardized biological metrics were then temporally and spatially aligned with corresponding winter climate variables, including minimum surface temperature (Tmin), snowpack depth and duration, and winter precipitation anomalies.

To account for variability in data quality and sampling design, a composite indexing approach was used in which spatially distributed population records were aggregated into representative range-wide time series where appropriate. Each dataset was weighted according to survey effort, spatial coverage, and methodological consistency, with greater influence assigned to long-term and systematically collected monitoring records. For cryptic species such as Clonophis kirtlandii, observational uncertainty was explicitly incorporated by representing detection variability as statistical variance within the input matrix. In these cases, presence–absence records and detection-based observations were integrated using probabilistic formulations to estimate latent demographic variables, including population density and survival probability.

Temporal continuity was preserved through conservative interpolation of limited data gaps using linear or spline-based methods. This approach maintained continuity in time series structure while minimizing the risk of introducing artificial trends or bias into the reconstructed ecological signals.

2.6.3 Structuring for deep learning model development

The integrated spatiotemporal matrix was structured to enable effective deep learning model training and validation. The chronological organization of the datasets allowed for temporal partitioning, where earlier decades were used for model training to capture historical climate–population relationships, while more recent observations were reserved for independent validation to evaluate predictive performance and generalization capability.

To enhance robustness in model interpretation and attribution, the dataset integrates multiple studies, survey designs, and experimental sources, thereby capturing intraspecific variability, spatial heterogeneity, and measurement uncertainty. This comprehensive structure provides a stable foundation for applying deep learning methods in combination with SHAP-based attribution, enabling reliable quantification of the contribution of climatic predictors to local extirpation risk.

2.6.4 Composite population index construction and validation

To address heterogeneity across population metrics (occupancy, abundance, survival, detection probability, and presence/absence), all variables were transformed into a standardized, dimensionless Population State Index (PSI) bounded within the interval (0,1). This transformation ensures comparability across heterogeneous data types while preserving relative ecological structure.

For each metric x_i(t), normalization was performed as:

where represents the normalized value at time t.

Binary variables (presence/absence and eDNA detection) were directly encoded as:

The composite Population State Index (PSI) was then computed as a weighted aggregation:

Where weights w_i​ were assigned based on two explicit criteria: (i) empirical reliability of the data source and (ii) ecological informativeness of the metric. Specifically, higher weights were assigned to long-term, systematically monitored variables such as occupancy and survival because these metrics are derived from repeated standardized sampling and directly reflect demographic persistence. Intermediate weights were assigned to abundance-related proxies due to moderate sensitivity to sampling design, while lower weights were assigned to detection-based measures such as eDNA because of their higher false-negative rates, stochastic detection probability, and dependence on environmental conditions.

This weighting strategy is therefore not arbitrary but grounded in measurement reliability and ecological interpretability, ensuring that more stable and biologically informative indicators contribute more strongly to the latent population state representation.

Ecologically, this transformation is justified because all included variables represent observable manifestations of an underlying latent population condition. Occupancy reflects spatial persistence, abundance captures demographic magnitude, survival reflects stability through time, and detection-based measures provide probabilistic evidence of presence. Mapping these heterogeneous observations onto a unified scale allows the PSI to approximate a continuous latent representation of population viability, enabling integrated analysis across disparate datasets.

To assess the sensitivity of results to the weighting scheme, a robustness analysis was conducted using alternative formulations, including equal weighting, reliability-weighted (baseline) configuration, and metric-exclusion scenarios (e.g., removing eDNA or occupancy variables). Model performance was evaluated using RMSE, AUC-ROC, and SHAP-based feature ranking stability. Across all tested configurations, variations in predictive performance remained within 5–8%, and the relative importance of key climatic drivers (snowpack, temperature, and precipitation) remained unchanged.

These results confirm that the weighting strategy is both ecologically justified and methodologically robust, and that model outcomes are not driven by arbitrary or unstable weighting assumptions.

2.7 Data preprocessing and model validation (revised and aligned)

To ensure methodological consistency and robust comparison across modeling approaches, a standardized preprocessing and validation pipeline was implemented for all models, including the deep learning framework and baseline models (GLM, GAM, Random Forest, and XGBoost).

2.7.1 Data normalization and transformation

All continuous input features, including snow depth, temperature anomalies, and winter precipitation, were standardized using z-score normalization to ensure comparability across variables and to eliminate scale-dependent bias during model training.

To address heterogeneity among population metrics, all biological response variables were transformed into a dimensionless population persistence index (PPI) bounded between 0 and 1. Occupancy and presence–absence data were represented as binary values, survival rates were retained as proportional measures, abundance data were normalized using species-specific min–max scaling, and detection-based proxies such as eDNA and low-detection counts were converted into probabilistic estimates based on detection likelihood.

This unified transformation framework enables consistent model learning across diverse ecological indicators while preserving their relative biological significance.

2.7.2 Handling missing data

This conservative strategy minimizes bias while maintaining temporal continuity. Missing observations were handled using a conservative, tiered strategy. Small gaps representing less than 10% of a time series were imputed using seasonal smoothing methods, including moving averages or spline interpolation. Larger gaps were excluded to prevent the introduction of artificial temporal patterns. This approach minimizes bias while preserving the integrity and continuity of the data.

2.7.3 Train–test splitting and cross-validation

To ensure fair and reproducible comparison across all modeling approaches:

• A consistent 80/20 train–test split was applied across all models

• A 5-fold cross-validation strategy was implemented using identical folds for each model

• Temporal structure was preserved where applicable, ensuring that training data preceded validation data to avoid information leakage

This unified framework guarantees that differences in performance reflect model capability rather than data partitioning artifacts.

2.7.4 Evaluation metrics

Model performance was evaluated using identical metrics across all approaches:

• Continuous outputs (population persistence, survival probability):

  • Root Mean Square Error (RMSE)

  • Coefficient of Determination (R²)

• Binary classification (extirpation events):

  • Area Under the Receiver Operating Characteristic Curve (AUC-ROC)

  • F1-score

These metrics enable direct and objective comparison between statistical, ensemble, and deep learning models.

2.7.5 Sensitivity analysis and robustness checks

To assess the stability of model outputs under alternative preprocessing assumptions, sensitivity analyses were conducted:

• Alternative normalization strategies (min–max scaling vs. z-score)

• Variation in composite index construction (species-specific vs. global normalization)

• Inclusion/exclusion of detection-based proxies

Results indicated that model performance varied by less than ±5% across preprocessing configurations, and the relative importance of snowpack and temperature variables remained consistent.

This demonstrates that the primary conclusions are robust to reasonable variations in data transformation and integration.

3 Results

3.1 Model performance and predictive accuracy

The fully connected neural network (FCNN) demonstrated strong predictive performance across both regression and classification tasks, indicating its effectiveness in modeling complex climate–ecological relationships. For continuous population predictions, the model achieved root mean square error (RMSE) values ranging from 0.08 to 0.15 and coefficients of determination (R²) between 0.71 and 0.87, reflecting high accuracy in capturing population dynamics across species and regions.

For extirpation classification, the FCNN exhibited robust discriminatory ability, with area under the receiver operating characteristic curve (AUC-ROC) values ranging from 0.84 to 0.93 and F1-scores between 0.79 and 0.88. These results indicate reliable identification of high-risk conditions and strong overall classification performance.

Comparative model evaluation shown in Figure 8 further highlights the advantages of the FCNN over traditional statistical approaches. As shown in Table 6, the FCNN consistently outperformed generalized linear models (GLM) and generalized additive models (GAM) across all evaluation metrics. RMSE values were substantially lower for the FCNN (0.08–0.15) compared to 0.18–0.28 for GLM and GAM models, while R² values improved from 0.61–0.73 in baseline models to 0.71–0.87 in the FCNN. These improvements demonstrate that deep learning more effectively captures nonlinear responses and interaction effects among winter climate variables, which are critical for accurately modeling extirpation risk.

Figure 8

To assess model robustness, bootstrap resampling (n = 1,000 iterations) was conducted to estimate uncertainty in predictive performance. The resulting confidence intervals were narrow (± 0.01–0.03 RMSE), indicating stable and consistent model behavior across resampled datasets. Overall, these results confirm that the FCNN provides both high predictive accuracy and strong generalization capability for modeling climate-driven population dynamics.

3.2 Quantitative feature attribution and relative importance

SHAP-based feature attribution illustrated in Figure 9 revealed clear differences in the contribution of winter climate variables:

Figure 9

• Snowpack (depth and duration): 42–57% mean absolute SHAP contribution

• Temperature anomalies: 28–39% contribution

• Winter precipitation: 12–21% contribution

Snowpack consistently emerged as the dominant predictor, contributing ~1.3× more than temperature and 1.4–2.1× more than precipitation alone. Global SHAP distributions indicate that snow depth below ~30% of the historical mean corresponds to sharply increased extirpation probability, illustrating strong nonlinear thresholds.

3.3 Interaction effects and synergistic amplification

Interaction analysis presented in Figure 10 showed that combined climatic stressors amplified extirpation risk beyond additive expectations. The snowpack vs temperature interaction increased SHAP magnitude by 18–26% relative to independent effects. Low snowpack combined with high temperature anomalies elevated predicted extirpation probability by 35–62%, demonstrating threshold and nonlinear escalation in risk.

Figure 10

3.4 Species-specific quantitative sensitivity

3.4.1 Acris blanchardi

Snowpack contributed 48–55% of total SHAP importance. A 10% reduction in snow depth increased predicted mortality by 6–11%. Temperature anomalies contributed ~30%, with higher winter temperatures increasing variability in survival predictions, while precipitation contributed ~15–20%, mainly affecting hydroperiod-related variability as illustrated in Figure 11.

Figure 11

3.4.2 Sistrurus catenatus

Snowpack contributed 44–52%, temperature 31–37%. High winter temperature variability (>3 °C SD) increased predicted mortality risk by 22–38%, and combined stressors increased extirpation probability by up to 55% as shown in Figure 12.

Figure 12

3.4.3 Clonophis kirtlandii

Snowpack contributed 50–57%, with strong interaction effects. Combined stress scenarios increased extirpation probability by up to 68%, amplified by low baseline density and fragmentation (prediction uncertainty +25%) shown in Figure 13.

Figure 13

3.4.4 Regina septemvittata

Snowpack contributed 42–48%; snowpack and temperature interactions increased SHAP values by ~24%, while precipitation effects contributed 18–22%. Populations were most sensitive under combined reduced snowpack and increased winter rainfall, with persistence decreasing by 30–47% as shown in Figure 14.

Figure 14

3.5 Spatiotemporal population trends

Standardized indices revealed multi-decadal declines:

• Amphibians: 35–65%

• Reptiles: consistently low densities (<150 individuals/site)

• Habitat patch reduction: 40–70%

Post-2000 interannual variance increased by 20–35%, consistent with rising climate variability (Figure 15).

Figure 15

3.6 Threshold behavior and nonlinear climate response

The fully connected neural network (FCNN) identified clear ecological thresholds associated with rapid increases in extirpation risk. Snowpack reductions below approximately 30–40% of the historical mean were associated with a sharp rise in extirpation probability, indicating a critical loss of subnivean thermal buffering. Similarly, temperature anomalies exceeding +2 °C were linked to amplified metabolic stress during overwintering periods.

When these thresholds occurred simultaneously, their effects were not additive but strongly nonlinear, resulting in a two- to threefold acceleration in extirpation risk compared to linear model expectations. These findings highlight the presence of interacting climatic tipping points, where relatively small additional changes in winter conditions can trigger disproportionately large ecological responses. Importantly, the magnitude and sensitivity of these thresholds varied across species and were influenced by ecological traits such as habitat specialization and population structure (Figure 16).

Figure 16

3.7 Model-based attribution of extirpation risk

High-risk predictions (>0.7 probability) were typically associated with:

• Snowpack reduction >50% of SHAP contribution

• Temperature variability 25–35%

• Precipitation 10–20%, primarily indirect

Combined interactions contributed up to 25% of total model output, highlighting multi-variable coupling (Figure 17).

Figure 17

3.8 Integrated climate–population response dynamics

FCNN outputs indicate that winter climate predictors influence extirpation risk through multiple interacting pathways. Thermal stress amplification contributes an estimated effect size increase of approximately +20–40%, while hydrological alterations account for an additional +15–30% contribution. In addition, compound interactions among climatic variables produce a further amplification of +30–60%, highlighting the nonlinear nature of these effects. Collectively, these results demonstrate that combined winter climate changes exert a disproportionately strong impact on population persistence, as illustrated in Figure 18.

Figure 18

4 Discussion

4.1 Winter climate as a primary driver of extirpation risk

This study provides strong quantitative evidence that winter climate dynamics, particularly snowpack loss and increasing temperature variability, are dominant predictors of amphibian and reptile extirpation risk across the Midwestern United States. Snowpack variables contributed the largest proportion of model importance (42–57%), exceeding both temperature and precipitation effects. These findings are consistent with prior work demonstrating the ecological importance of snowpack as a thermal buffer in overwintering systems (Muths et al., 2020; Johnston et al., 2021). For example, Johnston (Johnston et al., 2021) documented rapid population declines in a subnivean hibernator following snow drought conditions, highlighting the critical role of snow cover in maintaining stable winter microclimates.

More broadly, these results extend existing climate–biodiversity frameworks, which have traditionally emphasized mean temperature increases, by demonstrating that winter instability is a key driver of population decline. This aligns with recent studies showing that increased temperature variability and climate “reddening” can elevate extinction risk beyond the effects of gradual warming alone (Alkhayuon et al., 2022; Duffy et al., 2022). Together, these results position winter climate instability as a central, yet underrecognized, driver of biodiversity loss in temperate ecosystems.

4.2 Mechanistic pathways linking climate to population decline

The integration of deep learning with SHAP-based attribution provides quantitative support for key ecological mechanisms linking winter climate variability to population decline. These pathways are consistent with established ecological theory (Cahill et al., 2013; Urban, 2015) and empirical observations across amphibian and reptile systems.

Thermal stress and metabolic disruption arise from increased winter temperatures and variability, which elevate metabolic demands during dormancy and accelerate the depletion of energy reserves. This mechanism has been widely documented in amphibians, where altered thermal regimes directly affect survival rates and physiological condition (Vanselow et al., 2019). The strong sensitivity observed in Acris blanchardi aligns with prior studies showing that environmental stressors significantly influence both survival and developmental stability in this species (Mitchell and Tully, 2005; Muths et al., 2020).

A second pathway involves the loss of subnivean insulation due to reduced snowpack, which eliminates the thermal buffering capacity of overwintering habitats and exposes organisms to repeated freeze–thaw cycles and extreme temperature fluctuations. Similar processes have been observed in subnivean and soil-dwelling species, where reduced snow cover increases mortality risk through heightened thermal instability (Ren et al., 2017; Lehtinen et al., 2025). This mechanism is particularly relevant for species such as Sistrurus catenatus and Clonophis kirtlandii, which depend on stable overwintering microhabitats for survival.

Hydrological and habitat alteration represents a third major pathway, driven by shifts from snow-dominated to rain-dominated winter precipitation regimes. These changes modify hydrological processes, including wetland hydroperiods and stream flow dynamics, which in turn influence breeding habitat availability, overwintering refugia, and prey dynamics. Semi-aquatic species such as Regina septemvittata are especially affected by these alterations, consistent with broader ecological evidence showing that climate-driven hydrological shifts can restructure aquatic and semi-aquatic ecosystems (Illinois Natural Heritage Program, 2024; Liu et al., 2024).

Importantly, these mechanisms do not operate independently but instead interact to produce compounded ecological stress, consistent with multi-driver extinction theory (Brook et al., 2008). This interaction structure is supported by SHAP-derived effect sizes, where thermal variables account for 28–39% of model-attributed mortality risk, while interaction effects further amplify total impacts by 18–26%.

Top of Form.

Bottom of Form.

4.3 Synergistic effects and nonlinear climate responses

A key contribution of this study is the identification of nonlinear interactions between climate variables, particularly between snowpack and temperature. The results show that combined stressors increase extirpation risk by up to 62%, far exceeding additive expectations. This supports prior work demonstrating that synergistic interactions among environmental stressors can accelerate biodiversity loss (Brook et al., 2008; Velasco et al., 2021).

The observed threshold behavior, where extirpation risk increases rapidly beyond specific climatic limits, is consistent with the concept of ecological tipping points (Scheffer et al., 2001; Vanselow et al., 2019). Similar nonlinear dynamics have been reported in climate-driven extinction models, where gradual environmental change leads to abrupt population collapse once critical thresholds are exceeded (Alkhayuon et al., 2022). These findings reinforce the idea that species responses to climate change are often discontinuous rather than gradual.

Beyond interactions among climatic variables, the results have important implications for broader multi-driver extinction dynamics. Winter climate stress does not operate in isolation but can amplify the effects of non-climatic pressures such as habitat fragmentation, contamination, and disease. For example, fragmented habitats reduce access to thermally stable overwintering refugia, increasing exposure to freeze–thaw stress. Similarly, pesticide exposure prior to winter can reduce physiological resilience, lowering the capacity of amphibians to tolerate metabolic and thermal stress during dormancy. Disease dynamics may also interact with winter conditions, as physiologically stressed individuals emerging from energetically costly winters may exhibit reduced immune function, increasing susceptibility to infection. These combined effects suggest that winter climate variability acts as a stress multiplier, intensifying the impact of co-occurring environmental pressures and accelerating extirpation risk beyond climate-only predictions.

4.4 Species-specific vulnerability and ecological traits

Although all focal species exhibited sensitivity to winter climate dynamics, the magnitude and nature of responses varied according to ecological traits. Species with low population density and high fragmentation, such as Clonophis kirtlandii, exhibited the highest sensitivity and variability in predicted outcomes. This is consistent with empirical studies showing that fragmented populations are more vulnerable to environmental stochasticity and local extinction (Ibach, 2022; Stewart et al., 2023).

Habitat specialists, including wetland-dependent and stream-associated species, showed stronger responses to hydrological changes, aligning with research demonstrating that ecological specialization increases sensitivity to environmental change (Baker et al., 2018). In contrast, species with broader ecological tolerance, such as Regina septemvittata, exhibited lower sensitivity to individual variables but stronger responses to combined stressors, supporting findings that generalist species may be more resilient to single stressors but still vulnerable to compound effects (Cahill et al., 2013).

4.5 Implications for amphibian and reptile decline

This study contributes to a growing body of evidence linking climate change to amphibian and reptile declines (Walther et al., 2002; Parmesan, 2006; Urban, 2015). However, unlike many previous studies that emphasize summer temperature increases or precipitation changes, this work highlights the critical role of winter processes. Recent studies have begun to recognize the importance of winter severity and snowpack in shaping amphibian population dynamics (Muths et al., 2020; Weiskopf et al., 2022), but these factors remain underrepresented in broader climate–biodiversity research.

Importantly, the winter-driven mechanisms identified in this study are likely to interact with additional ecological stressors that were not explicitly modeled but are well documented in amphibian and reptile decline. Habitat fragmentation can limit dispersal and restrict access to suitable overwintering microhabitats, thereby amplifying the effects of snowpack loss and thermal instability. Likewise, sublethal exposure to contaminants such as pesticides may impair physiological condition prior to winter, reducing survival under energetically demanding overwinter conditions. Disease processes, while often temperature-dependent, may exert indirect effects through winter-mediated immune suppression or post-winter susceptibility, particularly following years with high climatic variability. These interactions highlight that the observed climate effects likely represent a conservative estimate of real-world extinction risk, where multiple stressors act simultaneously.

The mechanisms identified here particularly snowpack loss and freeze–thaw instability provide a plausible explanation for observed patterns of local extirpation in temperate regions. These findings are consistent with long-term observational studies reporting dynamic population declines and turnover in species such as Acris blanchardi (Lehtinen et al., 2025). Collectively, this study suggests that winter climate change is an underrecognized but critical driver of biodiversity loss.

4.6 Methodological contributions and advantages

This study demonstrates the value of combining deep learning with interpretable AI techniques such as SHAP for ecological modeling. Traditional statistical approaches (e.g., GLMs, GAMs) often struggle to capture nonlinear interactions and complex dependencies among variables, particularly in climate systems (Williams et al., 2007). While transformer and graph-based models are emerging as powerful alternatives, they typically require large, highly structured spatiotemporal datasets and explicit relational connectivity, which are not available in this study.

Compared to GLM and GAM baselines, the FCNN provides both higher predictive accuracy and the ability to identify nonlinear and synergistic climate effects, which are critical for mechanistic ecological interpretation. A limitation of the present study is the absence of direct comparison with state-of-the-art (SOTA) architectures; however, the heterogeneous, partially sparse, and species-aggregated nature of the dataset limits the feasibility and interpretability of such models in this context. As a result, the FCNN framework represents a pragmatic balance between predictive performance, data availability, and interpretability.

This approach addresses a key limitation in ecological modeling by bridging the gap between predictive performance and mechanistic understanding, a growing concern in ecological and climate modeling research.

4.7 Limitations and sources of uncertainty

Despite its strengths, this study has several limitations that should be considered when interpreting the results. The integration of heterogeneous datasets from multiple sources may introduce inconsistencies, even with normalization procedures. In addition, spatial interpolation of climate variables may obscure fine-scale microclimatic variation, which can play a critical role in overwinter survival (Johnston et al., 2021). Detection challenges for cryptic species such as Clonophis kirtlandii may also influence the accuracy of population estimates (Stewart et al., 2023).

Although cross-validation indicates strong model performance, caution is warranted when extrapolating beyond the studied region, species, or environmental conditions. Furthermore, the study does not include benchmarking against emerging state-of-the-art architectures, such as transformer or graph-based models. While these approaches may offer advantages in highly structured spatiotemporal contexts, their substantial data requirements and reduced interpretability limit their applicability in the present system.

Future research incorporating higher-resolution climate data, improved detection methods, experimental validation, and broader species representation would help reduce uncertainty and strengthen the generality of these findings.

4.8 Interpretability vs. causality: limitations of SHAP

While SHAP values provide a robust and interpretable framework for quantifying feature contributions to model predictions, it is critical to distinguish between statistical attribution and biological causation. SHAP estimates how individual variables influence model output relative to a baseline prediction, but these contributions reflect patterns learned from the data rather than definitive causal relationships (Cahill et al., 2013; Urban, 2015).

A key limitation lies in the distinction between correlation and causation. The model identifies snowpack loss and temperature variability as highly influential predictors of extirpation risk; however, SHAP attribution alone cannot confirm that these variables directly cause population decline without independent experimental or longitudinal validation. In some cases, high attribution values may reflect proxy relationships, where climatic variables are correlated with unmeasured environmental or anthropogenic factors.

Additionally, SHAP values can be sensitive to feature dependency. Climatic variables such as temperature, precipitation, and snowpack are inherently interrelated, and multicollinearity may influence how importance is distributed among predictors. As a result, the relative contribution assigned to each variable may reflect the statistical structure of the dataset rather than fully independent biological effects.

More broadly, SHAP explains the behavior of the model rather than the underlying ecological system. If the deep learning model captures spurious or dataset-specific relationships, SHAP will accurately attribute importance to those patterns, potentially leading to overinterpretation of non-causal drivers. Therefore, while SHAP-based results provide strong quantitative support for hypothesized mechanisms, they should be interpreted as hypothesis-generating rather than definitive evidence of causality.

4.9 Conservation and management implications

The identification of snowpack loss and increasing winter climate variability as primary predictors of extirpation risk has direct and actionable implications for conservation and resource management. These findings underscore the importance of winter processes, which remain underrepresented in most climate adaptation frameworks that traditionally emphasize growing-season temperature and precipitation extremes.

From a management perspective, habitat conservation strategies should prioritize the protection and restoration of microhabitats that buffer winter thermal variability, including deep soil substrates, wetlands, forested cover, and landscape features that promote snow accumulation and retention. Maintaining conditions that preserve at least 30–40% of historical snowpack buffering capacity is particularly critical, as model-derived thresholds indicate a sharp increase in extirpation risk below this range. These subnivean buffering environments are essential for stabilizing overwinter physiological conditions and reducing exposure to repeated freeze–thaw stress.

Future conservation planning should also explicitly incorporate winter climate variables into both predictive modeling and monitoring frameworks. Rather than relying solely on annual or growing-season climate summaries, ecological studies should integrate high-resolution winter-specific predictors such as snow water equivalent, snow cover duration, freeze–thaw frequency, and winter temperature variability. These variables can be incorporated into species distribution models (SDMs) and occupancy models by treating winter conditions as seasonally stratified covariates rather than aggregated annual metrics. In spatial modeling applications, coupling gridded climate products (e.g., snowpack and temperature reanalysis datasets) with species occurrence data at fine temporal scales would allow for more realistic representation of overwinter habitat suitability and mortality risk.

In addition, the integration of remote sensing products (e.g., MODIS-derived snow cover and reanalysis-based snow water equivalent datasets) into spatially explicit ecological models offers a practical pathway for scaling winter climate effects across landscapes. These datasets can be directly embedded into GIS-based habitat suitability models or machine learning frameworks to capture spatial heterogeneity in snow persistence and winter thermal stability. Coupling these layers with hydrological and land-cover data would further improve predictive capacity for identifying winter refugia and high-risk extinction zones.

Monitoring programs should also be redesigned to better capture winter-driven population dynamics. Specifically, increased emphasis on overwinter survival rates, timing of spring emergence, and interannual variability in recruitment would provide earlier and more sensitive indicators of population decline. Incorporating winter sampling protocols into long-term ecological monitoring would help close a critical observational gap in current biodiversity assessments.

Overall, these results suggest that conservation strategies focused exclusively on growing-season conditions may substantially underestimate species vulnerability. Explicit integration of winter climate dynamics into ecological modeling, spatial forecasting, and management planning will be essential for improving predictive accuracy and conservation outcomes under ongoing climate change.

5 Conclusion

This study provides a comprehensive, data-driven demonstration that winter climate dynamics specifically snowpack loss and increasing temperature variability are dominant and quantifiable predictors of amphibian and reptile extirpation risk across the Midwestern United States. By integrating heterogeneous long-term population datasets with high-resolution climatic variables and applying an interpretable deep learning framework, this research moves beyond simple correlation by providing quantitative, model-based attribution that is consistent with hypothesized ecological mechanisms.

The developed neural network achieved strong predictive performance across both continuous and classification tasks, with RMSE values ranging from 0.08 to 0.15 and R² between 0.71 and 0.87, indicating high accuracy in modeling population persistence and survival probability. For extirpation classification, the model demonstrated robust discrimination, with AUC-ROC values of 0.84–0.93 and F1-scores of 0.79–0.88, confirming reliable identification of high-risk conditions across species and regions.

SHAP-based attribution analysis revealed that snowpack (42–57%) is the most influential predictor of extirpation risk, exceeding the contributions of temperature anomalies (28–39%) and winter precipitation (12–21%). Critically, the results highlight strong nonlinear and synergistic interactions: combined snowpack reduction and temperature increases amplified extirpation probability by 35–62%, with interaction effects contributing an additional 18–26% to model output variance. These findings confirm that extirpation risk is not driven by single variables in isolation but by compound winter stressors.

The model further identified clear ecological thresholds, beyond which population collapse accelerates rapidly. Specifically, snowpack declines below 30–40% of historical mean and temperature anomalies exceeding +2 °C triggered nonlinear increases in extirpation probability, with risk escalating 2–3× faster than linear expectations. These threshold dynamics were consistent across species, although their magnitude varied according to ecological specialization and population structure.

Species-specific analyses demonstrated differentiated vulnerability. For example, Acris blanchardi exhibited a 6–11% increase in mortality probability per 10% snowpack reduction, while Sistrurus catenatus showed 22–38% higher mortality risk under increased temperature variability. Clonophis kirtlandii emerged as the most sensitive species, with extirpation probability increasing by up to 68% under combined stressors, reflecting compounding effects of fragmentation and low population density. In contrast, Regina septemvittata showed stronger sensitivity to hydrological interactions, with population persistence declining by 30–47% under combined snowpack loss and precipitation shifts.

At broader scales, the integrated dataset revealed substantial long-term ecological change, including 35–65% declines in amphibian populations, persistent low reptile densities (<150 individuals per site), and 40–70% reductions in occupied habitat patches. Additionally, post-2000 variability increased by 20–35%, reflecting intensifying climate instability and reinforcing the link between winter climate dynamics and population fluctuations. These findings confirm that extirpation risk is driven by interacting pathways of thermal stress (20–40% effect size), hydrological alteration (15–30%), and compound interaction amplification (30–60%). These coupled processes disrupt overwinter survival, breeding phenology, and habitat stability, ultimately leading to population collapse.

Collectively, these findings demonstrate that winter climate change is not a secondary component of global warming, but a primary driver of ecological instability and species loss in temperate ecosystems.

References

• 1

  AlkhayuonH.AshwinP.JacksonL. C.QuinnA. D.WoodR. A. (2022). Stochastic resonance in climate reddening increases the risk of ecosystem extinction via phase-tipping. Proc. R. Soc A.478, 20220273. doi: 10.1098/rspa.2022.0273. PMID:

  • Pubmed Abstract
  • CrossRef
  • Google Scholar
  • View reference in article

• 2

  BakerS. J.MullowneyJ. M.AgostaS. J.HromadaR. J.SchigerM. J. (2018). Temporal patterns of genetic diversity in an imperiled population of eastern massasauga rattlesnakes. Ichthyol. Herpetol.106, 682–688. doi: 10.1643/CG-17-682

  • CrossRef
  • Google Scholar
  • View reference in article

• 3

  BrodieJ. F.PostE.BergerJ. (2014). Trophic interactions and dynamic herbivore responses to snowpack. Clim. Change Responses1, 4. doi: 10.1186/s40665-014-0004-2. PMID:

  • Pubmed Abstract
  • CrossRef
  • Google Scholar
  • View reference in article

• 4

  BrookB. W.SodhiN. S.BradshawC. J. A. (2008). Synergies among extinction drivers under global change. Trends Ecol. Evol.23, 453–460. doi: 10.1016/j.tree.2008.03.011. PMID:

  • Pubmed Abstract
  • CrossRef
  • Google Scholar
  • View reference in article

• 5

  BukaitaW. (2024). Global Warming’s Influence on Temperature Increase. In: AraiK.. (eds) Proceedings of the Future Technologies Conference (FTC) 2024, Volume 3. FTC 2024. Lecture Notes in Networks and Systems, vol 1156. Springer, Cham. doi: 10.1007/978-3-031-73125-9_18

  • CrossRef
  • Google Scholar
  • View reference in article

• 6

  BukaitaW.AnyaiweO.NelsonP. (2024). An analysis of temperature variability using an index model. Adv. Inf. Commun. 192–212. doi: 10.1007/978-3-031-54053-0_15

  • CrossRef
  • Google Scholar
  • View reference in article

• 7

  BukaitaW.GhiurauA. (2025). Multivariate climatic drivers of local extirpation. Int. J. Environ. Monit. Anal.13, 328–346. doi: 10.11648/j.ijema.20251306.14

  • CrossRef
  • Google Scholar
  • View reference in article

• 8

  CahillA. E.Aiello-LammensM. E.Fisher-ReidM. C.HuaX.KaranewskyC. J.RyuH. Y.et al. (2013). How does climate change cause extinction? Proc. R. Soc B.280, 20121890. doi: 10.1098/rspb.2012.1890. PMID:

  • Pubmed Abstract
  • CrossRef
  • Google Scholar
  • View reference in article

• 9

  CrawfordJ. A. (2024). Conservation of Kirtland’s snake – a wet prairie species specialist. Outdoor Illinois J. (2024). Available online at: https://outdoor.wildlifeillinois.org/articles/conservation-of-kirtlands-snake-a-wet-prairie-species-specialist?utm_source=chatgpt.com (Accessed March 10, 2026).

  • Google Scholar

• 10

  DuffyK. A.GouhierT. C.GangulyA. R. (2022). Climate-mediated shifts in temperature fluctuations promote extinction risk. Proc. R. Soc B.289, 20220231. doi: 10.1038/s41558-022-01490-7. PMID:

  • Pubmed Abstract
  • CrossRef
  • Google Scholar
  • View reference in article

• 11

  Government of CanadaSpecies at risk consultation: eastern massasauga. Available online at: https://www.Canada.ca (Accessed March 10, 2026).

  • Google Scholar
  • View reference in article

• 12

  HilemanE. T.KingR. B.FaustL. J. (2018). Demography and extinction risk of eastern massasauga rattlesnakes under different management scenarios. J. Wildl. Manage.82, 965–977. doi: 10.1002/jwmg.21457. PMID:

  • Pubmed Abstract
  • CrossRef
  • Google Scholar
  • View reference in article

• 13

  HoskinsT. D.BooneM. D. (2017). Variation in larval survival of Blanchard’s cricket frog (Acris blanchardi) under pesticide exposure. Environ. Toxicol. Chem.36, 1613–1620. doi: 10.1002/etc.3715. PMID:

  • Pubmed Abstract
  • CrossRef
  • Google Scholar
  • View reference in article

• 14

  HoskinsT. D.BooneM. D. (2018). Atrazine exposure alters sex ratios in Blanchard’s cricket frog (Acris blanchardi). Environ. Toxicol. Chem.37, 471–478. doi: 10.1002/etc.3962. PMID:

  • Pubmed Abstract
  • CrossRef
  • Google Scholar
  • View reference in article

• 15

  Ibach (2022). Effects of snake fungal disease on the survival and growth of the queensnake. doi: 10.13023/etd.2022.348.

  • CrossRef
  • Google Scholar
  • View reference in article

• 16

  Illinois Natural Heritage Program (2024). Species guidance: Kirtland’s snake (Springfield, IL: Illinois Department of Natural Resources).

  • Google Scholar
  • View reference in article

• 17

  IPCC (2021). Climate change 2021: the physical science basis (Cambridge: Cambridge University Press). doi: 10.1017/9781009157896

  • CrossRef
  • Google Scholar
  • View reference in article

• 18

  JohnstonA. N.ChristophersenR. G.BeeverE. A.RansomJ. I. (2021). Freezing in a warming climate: marked declines of a subnivean hibernator after a snow drought. Ecol. Evol.11, 1264–1279. doi: 10.1002/ece3.7126. PMID:

  • Pubmed Abstract
  • CrossRef
  • Google Scholar
  • View reference in article

• 19

  JonesP. C.KingR. B.SuttonS. (2017). Demographic analysis of imperiled eastern massasaugas (Sistrurus catenatus). J. Herpetol.51, 264–270. doi: 10.1670/15-058. PMID:

  • Pubmed Abstract
  • CrossRef
  • Google Scholar
  • View reference in article

• 20

  KingR. B.et al. (2012). Range-wide analysis of eastern massasauga survivorship and demography. J. Wildl. Manage.76, 434–444. doi: 10.1002/jwmg.418. PMID:

  • Pubmed Abstract
  • CrossRef
  • Google Scholar
  • View reference in article

• 21

  LehtinenR. M.KrynakK. L.Lipps JrG. J.McCallJ. C.YoungquistM. B. (2025). Hotspots, population turnover, and long-term data reveal the dynamic nature of Blanchard’s cricket frog populations. Ecol. Evol.15, e72569. doi: 10.1002/ece3.72569. PMID:

  • Pubmed Abstract
  • CrossRef
  • Google Scholar
  • View reference in article

• 22

  LeuenbergerW.DavisA. G.McKenzieJ. M.DrayerA. N.PriceS. J. (2019). Evaluating snake density using PIT telemetry and spatial capture–recapture analyses. J. Herpetol.53, 272–281. doi: 10.1670/18-070. PMID:

  • Pubmed Abstract
  • CrossRef
  • Google Scholar
  • View reference in article

• 23

  LiuH.YiS.DingY.KangS.WangJ. (2024). Winter snowpack loss increases warm-season compound hot-dry extremes. Commun. Earth Environ.5, 567. doi: 10.1038/s43247-024-01734-8. PMID:

  • Pubmed Abstract
  • CrossRef
  • Google Scholar
  • View reference in article

• 24

  McCallumM. L.TrauthS. E. (2021). Habitat factors and restoration success in Blanchard’s cricket frog (Acris blanchardi). J. North. Am. Herpetol.2021, 1–12. doi: 10.17161/jnah.vi.14739

  • CrossRef
  • Google Scholar
  • View reference in article

• 25

  McKenzieJ. M.PriceS. J.ConnetteG. M.BonnerS. J.LorchJ. M. (2021). Effects of snake fungal disease on survival and movement. Ecol. Appl.31, e2251. doi: 10.1002/eap.2251. PMID:

  • Pubmed Abstract
  • CrossRef
  • Google Scholar
  • View reference in article

• 26

  MitchellM. A.TullyT. N. (2005). Salmonella in free-ranging reptiles from Pennsylvania. J. Wildl. Dis.41, 617–622. doi: 10.1053/saep.2001.19798

  • CrossRef
  • Google Scholar
  • View reference in article

• 27

  MuthsE. L.HossackB. R.GrantE. H.PilliodD. S.MosherB. A. (2020). Effects of snowpack, temperature, and disease on amphibian demography. Herpetologica76, 132–143. doi: 10.1655/0018-0831-76.2.132. PMID:

  • Pubmed Abstract
  • CrossRef
  • Google Scholar
  • View reference in article

• 28

  National Centers for Environmental Information (NCEI)Global summary of the month (GSOM). Available online at: https://www.ncei.noaa.gov(Accessed March 10, 2026).

  • Google Scholar
  • View reference in article

• 29

  ParmesanC. (2006). Ecological and evolutionary responses to recent climate change. Annu. Rev. Ecol. Evol. Syst.37, 637–669. doi: 10.1146/annurev.ecolsys.37.091305.110100. PMID:

  • Pubmed Abstract
  • CrossRef
  • Google Scholar
  • View reference in article

• 30

  PeclG. T.AraújoM. B.BellJ. D.BlanchardJ.BonebrakeT. C.ChenI. C.et al. (2017). Biodiversity redistribution under climate change: impacts on ecosystems and human well-being. Science355, eaai9214. doi: 10.1126/science.aai9214. PMID:

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 31

  PostE.StensethN. C. (1999). Climatic variability, plant phenology, and northern ungulates. Ecology80, 1322–1339. doi: 10.1890/0012-9658(1999)080[1322:CVPPAN]2.0.CO;2

  • CrossRef
  • Google Scholar

• 32

  RatschR.KingsburyB. A.JordanM. A. (2020). Environmental DNA detection of Kirtland’s snake (Clonophis kirtlandii). Animals10, 1057. doi: 10.3390/ani10061057. PMID:

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 33

  RenZ.ZhaoQ.ZhangS.KanJ. (2017). Winter is changing: trophic interactions under altered snow regimes. Food Webs13, 1–7. doi: 10.1016/j.fooweb.2017.02.006. PMID:

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 34

  SchefferM.CarpenterS. R.FoleyJ. A.FolkeC.WalkerB. (2001). Catastrophic shifts in ecosystems. Nature413, 591–596. doi: 10.1038/35098000. PMID:

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 35

  Stanley (2007). Distribution of the Queen Snake (Regina septemvittata). J. Arkansas Acad. Sci. 61, 103–108. doi: 10.54119/jaas.2007.6113.

  • CrossRef
  • Google Scholar

• 36

  StewartT. M.DreslikM. J.PhillipsC. A.KuhnsA. R.CzarneckiJ.KleopferJ. D.et al. (2023). Estimating the effort required to detect Kirtland’s snakes (Clonophis kirtlandii). Wildl. Soc Bull.51, e1498. doi: 10.1002/wsb.1498. PMID:

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 37

  United States Geological Survey (USGS)Status and trends of Blanchard’s cricket frog in the Upper Midwest. Available online at: https://www.usgs.gov (Accessed March 10, 2026).

  • Google Scholar

• 38

  United States Geological Survey (USGS) (2018). Kirtland’s snake range dataset. ScienceBase. doi: 10.5066/F7WS8SBN

  • CrossRef
  • Google Scholar

• 39

  United States Geological Survey (USGS) (2023). Genotype data for eastern massasauga from Bois Blanc Island, Michigan. ScienceBase. doi: 10.5066/P9HJW59U

  • CrossRef
  • Google Scholar

• 40

  UrbanM. C. (2015). Accelerating extinction risk from climate change. Science348, 571–573. doi: 10.1126/science.aaa4984. PMID:

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 41

  VanselowA.WieczorekS.FeudelU. (2019). When very slow is too fast: collapse of a predator–prey system. J. Theor. Biol.479, 1–12. doi: 10.1016/j.jtbi.2019.07.022. PMID:

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 42

  VelascoJ. A.EstradaF.Calderón-BustamanteO.SwingedouwD.UretaC.GayC.et al. (2021). Synergistic impacts of global warming and thermohaline circulation collapse on amphibians. Commun. Biol.4, 141. doi: 10.1038/s42003-021-01665-6. PMID:

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 43

  WaltherG. R.PostE.ConveyP.MenzelA.ParmesanC.BeebeeT. J. C.et al. (2002). Ecological responses to recent climate change. Nature416, 389–395. doi: 10.1038/416389a. PMID:

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 44

  WeiskopfS. R.ShiklomanovA. N.ThompsonL.WheedletonS.GrantE. H. C. (2022). Winter severity affects occupancy of anurans. Divers. Distrib.28, 2187–2199. doi: 10.1111/ddi.13620. PMID:

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 45

  WilliamsJ. W.JacksonS. T.KutzbachJ. E. (2007). Projected distributions of novel and disappearing climates by 2100 AD. Proc. Natl. Acad. Sci. U.S.A.104, 5738–5742. doi: 10.1073/pnas.0606292104. PMID:

  • Pubmed Abstract
  • CrossRef
  • Google Scholar