This study was conducted across the entire mainland of Bangladesh (Fig. 1B), situated between 23°41'N and 90°20'E, with a total land area of approximately 147,570 km² (Rashid 2019). Positioned at the confluence of the Indo-Himalayan and Indo-Chinese sub-regions within the Oriental biogeographic realm, Bangladesh serves as a critical corridor for species migration between the Indian subcontinent and Southeast Asia (Khan 2018). Its inclusion in the Indo-Burma biodiversity hotspot underscores its global ecological significance (Myers et al. 2000).

Figure 1. 

A. Distribution of Daboia russelii in Asia (Suraj et al. 2021); B. Occurrence points of D. russelii used for species distribution modelling; C. 12 bio-ecological zones (Nishat et al. 2002); D. Seven climatic sub-regions (Rashid 2019); E. Elevation (United States Geological Survey 2019); F. Flood plain (Bangladesh Agricultural Research Council 2019); G. Forest zones (Henry et al. 2021); H. Land use types (Buchhorn et al. 2019); I. River zones (Sarker et al. 2016). The details of the seven spatial variables are provided as Supplementary material 1.

Bangladesh experiences a subtropical monsoon climate, marked by three distinct seasons: pre-monsoon (March-May), monsoon (June to October), and winter (November to February), with heavy rainfall dominating the monsoon period (Rashid 2019). Its landscape is predominantly shaped by the dynamic river systems of the Brahmaputra, Ganges (Padma), and Meghna, forming an extensive riverine and deltaic terrain (Sarker et al. 2016).

The ecological diversity of Bangladesh is reflected in its 12 bio-ecological zones (Nishat et al. 2002), seven climatic sub-regions (Rashid 2019; Chowdhury et al. 2022), and five distinct forest types, encompassing tropical and subtropical ecosystems (Henry et al. 2021). This combination of climatic variability and topographic complexity makes Bangladesh an ideal region for studying ecological processes and species distribution patterns.

A trained snake rescuer team has been actively rescuing D. russelii from the western regions of the country from 2019 to 2023 through a continuous field survey and responding to the rescue call. The specimens of D. russelii were collected randomly from various locations along with their habitat notes, highlighting its occurrence in previously unreported areas. We compiled a total of 231 occurrence points of D. russelii from three sources: field collections, research-grade database, and reliable print or electronic media (Fig. 1B). Between 2019 and 2022, 69 specimens were collected through active field work across various localities in Bangladesh under a permit from the Bangladesh Forest Department (permission letter no: 22.01.0000.101.23.2019.3173). These specimens are currently housed at the Venom Research Centre, Bangladesh, and were documented with detailed information, including collection date, geographical coordinates, sex, and morphological characteristics. To supplement this dataset, area of occupancy (AOO) polygons were obtained from the IUCN Red List assessment (Rahman 2015), and 100 random occurrence points were collected from these polygons. Additional verifiable data were obtained from global biodiversity repositories such as GBIF (GBIF Secretariat 2023) and iNaturalist (iNaturalist 2024). After rigorous validation and cross-referencing, we incorporated 62 occurrence points derived from literature, preserved specimens, and authenticated sources in print, electronic, and social media.

We obtained 19 bioclimatic variables from WorldClim 2.0 at a spatial resolution of 30 arc-seconds (approximately 1 km² grid cells) for both current (1970–2000) and future (2061–2080) climate scenarios (Hijmans et al. 2005). To ensure consistency and minimise biases due to varying measurement scales, all variables were standardised using z-score transformation. Highly correlated variables (correlation coefficient, r > 0.9) were excluded to address multicollinearity concerns (De Marco and Nóbrega 2018), resulting in a final selection of 12 bioclimatic variables (Table 1). Future climate suitability for D. russelii was modelled using 23 climate projections for 2080 based on the shared socioeconomic pathway 370 (ssp370) scenario (Masson-Delmotte et al. 2021). All vector and raster datasets were prepared within the WGS84 coordinate system and clipped to match Bangladesh’s geographical boundaries.

In addition, seven spatial variables (e.g., bio-ecological zones, climatic sub-regions, elevation, flood plains, forest types, land use, and river zones) were included in the analysis (Table 1). These datasets were converted into raster format following the methodologies by Chowdhury et al. (2022) and Dutta et al. (2024). Autocorrelation tests confirmed that the spatial variables were independent. The 12 bio-ecological sub-zones and seven climatic sub-regions were digitised from published literature (Nishat et al. 2002; Rashid 2019). Floodplain and river zone data were sourced from the Bangladesh Agricultural Research Council (Bangladesh Agricultural Research Council 2019) and the Water Development Board, as cited by Sarker et al. (2016).

Elevation data were collected from the Earth Resources Observation and Science Centre using SRTM 3 arc-second (90 × 90 m²) datasets (United States Geological Survey 2019). Land use and land cover data, with a spatial resolution of 100 × 100 m², were obtained from the Copernicus Global Land Service (Buchhorn et al. 2019). Forest cover data, at a resolution of 30 × 30 m², were acquired from the Bangladesh Forest Department repository (Henry et al. 2021). All raster datasets were resampled to a standard resolution of 1 km² using the nearest neighbour resampling technique in the ArcGIS environment. These datasets featured categorical pixel values to represent distinct environmental attributes. For example, the flood raster categorised pixels as either outside the flood zone (value = 1) or within it (value = 2). Similarly, five other variables—bio-ecological zones, climatic sub-regions, forest zones, land use, and river zones—had complex categorical structures, with varying numbers of categories per raster file. This categorisation provided a detailed representation of the environmental characteristics across the study area, supporting comprehensive ecological analyses.

The spatial and temporal dispersion patterns of D. russelii were analysed using 69 field-collected occurrence records. To evaluate yearly range expansion, the occurrence data were organised chronologically and visualised using ArcGIS. Each record was sorted by year and mapped at the district level, with newly recorded districts added incrementally for each subsequent year. This stepwise approach enabled a progressive mapping of the viper’s distribution across Bangladesh over time. By examining the assumed yearwise distributional changes, potential routes of dispersion were inferred. This provided insights into possible geographic pathways and expansion trends of D. russelii, highlighting areas of recent colonisation and facilitating the identification of ecological factors driving this expansion.

To predict the suitable climatic niches of D. russelii in Bangladesh, we utilised Species Distribution Modelling (SDM) following the methodology described by Rangel and Loyola (2012). Initially, the six distinct spatial variables with categorical values—bio-ecological zones, climatic sub-regions, flood plain, forest types, land use, and river zones—were converted to raster with continuous values. This was done by calculating the proportion of the area of each category under the species’ Area of Occupancy (AOO). First, each categorical raster was masked using the AOO of this viper, retaining the number of categorical cells within the species’ AOO. The proportion of the AOO area within each category was then calculated to produce a category-specific percent occupancy. The percent occupancy values were then used to reclassify the categories by replacing each category value with the corresponding percent occupancy. For example, if the AOO of this viper species distributed 40% over the cells of category A and 60% over the cells of category B of a spatial raster, we assigned 0.4 to all cells of category A and 0.6 to all cells of category B. The rest of the cell of that raster was assigned as 0. Finally, raster stack of 19 variables (Table 1) were used to perform the species distribution model.

We employed the Bioclim () and the Maximum Entropy Maxent () algorithms to determine the suitable climate space for D. russelii in Bangladesh. The Bioclim () is a simple and non-parametric approach suitable for predicting climate-driven niches (Beaumont et al. 2005; Hallgren et al. 2019), while the Maxent () model is a machine-learning algorithm to determine the suitable climate space for D. russelii in Bangladesh (Phillips 2021). We used the ‘DISMO’ (Hijmans et al. 2017) and ‘RJAVA’ (Urbanek 2024), ‘MAXNET’ (Phillips 2021) packages on R platform (version 4.4.1) with the integration of several supporting packages, including ‘RASTER’ (Hijmans and Van Etten 2019), ‘MAPTOOL’ (Bivand et al. 2016), ‘RGDAL’ (Bivand et al. 2019), ‘SP’ (Pebesma and Bivand 2019) and ‘SF’ (Pebesma 2018). This approach provided a robust framework for identifying potential climate-driven habitats of D. russelii, facilitating further analysis of its distribution dynamics across Bangladesh. To ensure transparency, reproducibility, and comprehensive documentation of the species distribution modelling process, we followed the ODMAP (Overview, Data, Model, Assessment, and Prediction) protocol (Zurell et al. 2020). ODMAP provides a standardised framework for documenting each step of the modelling workflow, from data collection and pre-processing to model development, evaluation, and prediction. By adopting this structured approach, we aimed to enhance the clarity of our methodological decisions and facilitate future replication or extension of our study (Suppl. material 2). The models were evaluated using the area under curve (AUC), threshold-dependent sensitivity, specificity, and Cohen’s Kappa.

For this analysis, we first reclassified the predicted raster outputs from both the Bioclim and Maxent models using their respective threshold values. Cells with values above the threshold were assigned a value of 1 (indicating suitable areas), while cells with values below the threshold were assigned 0 (indicating unsuitable areas). These reclassified raster files were then ensembled. Since Maxent provided more realistic predictions with higher and more consistent AUCs for both current and future climate scenarios, it was assigned a higher weight in the ensemble operation. The ensemble formula was as follows:

Ensembled _{(current/future)} = (0.3 × reclassified Bioclim _{(current/future)}) + (0.7 × reclassified Maxent _{(current/future)})

The resulting ensembled raster was then reclassified into three distinct categories to indicate climatic suitability. To determine the threshold, we multiplied the threshold values from each model by their respective weights. For example, the threshold probability for the Bioclim prediction in the current climate is 0.0022, and for Maxent, it is 0.1356. Thus, the weighted threshold for the current climate is: (0.3 × 0.0022) + (0.7 × 0.1356) = 0.09558.

This weighted threshold was used to reclassify the ensembled predictions into three categories:

• Cells with no occurrence probability were assigned a value of 0.
• Cells with probability values from −∞ to 0.09558 were classified as 1 (unsuitable areas).
• Cells with probability values greater than 0.09558 were classified as 2 (suitable areas).

A similar process was applied to future climate model predictions, where the following categories were used:

• Cells with no occurrence probability were assigned a value of 0.
• Cells with probability values from −∞ to 0.03984 were classified as 3 (unsuitable areas).
• Cells with probability values greater than 0.03984 were classified as 6 (suitable areas).

This reclassification enabled a clear categorization of each cell within the raster for both the current and future climate conditions. To facilitate direct comparison between current and future suitability, we then overlapped the raster files from both time periods and summed the corresponding cell values to generate distinct integer values, as follows:

To identify the key factors influencing the distribution of D. russelii in Bangladesh, we performed a Principal Component Analysis (PCA) using SPSS, following the methodology described by McGarigal et al. (2000). PCA is a robust statistical tool that reduces dimensionality by identifying the most important variables that account for the maximum variance in a dataset. This allows us to isolate the climatic factors that most significantly drive the species’ distribution (Bueno de Mesquita et al. 2021; Joswig et al. 2022). Data were collected for 13 predictor variables, including 12 bioclimatic variables and elevation, at the occurrence points of D. russelii. Additionally, six spatial variables were converted from categorical to continuous format by calculating the percentage occupancy of each category within the species’ known range in Bangladesh. This transformation ensured consistency across variables, enabling a more accurate PCA. The results of the PCA highlighted the relative importance of climatic and environmental variables, providing valuable insights into the ecological drivers of the distribution of D. russelii and guiding subsequent habitat modelling efforts.

From 2019 to 2022, 61 vipers were collected during fieldwork and through rescue operations. Their sexes were determined using sex determination probes following established methodologies (Marais 1984; Mader 2006). This analysis aimed to evaluate the population’s sex ratio and its demographic structure over multiple years. However, data from 2022 and 2023 were excluded from the analysis due to incomplete records. By determining the sex ratio, we gained insights into the population dynamics of D. russelii, which could have implications for understanding its reproductive ecology.

Animals known to prey on snakes were identified as potential predators of D. russelii. Information on the conservation status and population trends of these predators (Suppl. material 3) was sourced from published literature, including the Red List of Bangladesh Volumes 2, 3, and 4 (IUCN Bangladesh 2015 a, b, c). Similarly, potential competitor species were identified based on their dietary overlap and shared habitat preferences with D. russelii (Suppl. material 4). Data on these possible competitors were gathered by analysing their feeding habits, behaviours, and geographic distributions as described in the same literature sources.

To assess habitat conditions, data on rice production and paddy fields, which serve as primary habitats for D. russelii, were collected from the Yearbook of Agricultural Statistics (Bangladesh Bureau of Statistics 2023) spanning the years 2012 to 2022. Additionally, information on annual floods, a key factor influencing snake dispersal and habitat dynamics, was extracted from the Annual Flood Report 2021 published by the Bangladesh Water Development Board (Barua et al. 2021). This multi-faceted approach enabled a comprehensive understanding of the ecological factors influencing the expansion and population trends of D. russelii, providing insights into its interactions with predators, competitors, and its flood-prone habitat.

The data underpinning the analysis reported in this paper are deposited at GBIF, the Global Biodiversity Information Facility, and are available at https://ipt.pensoft.net/resource?r=russells_viper_bangladesh.

Our field observation found that D. russelii might be dispersed through two possible ways—one exploiting aquatic channels and the other relying on terrestrial plains—enabling the species to expand its range across previously unoccupied habitats. As of 2019, the species’ distribution was primarily concentrated in the north-western regions of Bangladesh along the Padma and upper Meghna Rivers (Fig. 2A). From 2019 to 2021, dispersion was observed through riverine channels, where individuals moved from the upper to the lower Meghna River (Fig. 2B, C). This pathway might enable the viper to colonise new riverside localities, including tidal floodplains and coastal islands within the Bay of Bengal. On the other hand, sightings recorded in previously unoccupied localities such as Meherpur, Jhenaidah, and Chuadanga by 2022 can be deduced as terrestrial dispersion (Fig. 2D). By 2023, D. russelii had reached Satkhira, a coastal district in southwestern Bangladesh, through terrestrial corridors that included crossings via Jashore district (Fig. 2E).

Figure 2. 

Occurrence records of Daboia russelii in new localities in Bangladesh, illustrating its dispersal pattern: A. In 2019; B. 2020; and in C. 2021, this species was observed along the riverine channels of the Padma and Meghna Rivers; D. Since 2022, this species has occurred in new localities independent of river channels; E. In 2023, the species had reached coastal districts in the southwestern region of the country.

Climate model analysis revealed that D. russelii currently has a broad suitable climatic area of approximately 76,716 km² in western Bangladesh (Fig. 3A, B). However, projections for 2080 indicate a slight reduction in this suitable space for this viper. Across both present and future scenarios, approximately 75% of the study area remains unsuitable for D. russelii (Fig. 4).

Figure 3. 

Suitable climate map of D. russelii in Bangladesh. A. Weighted present suitable climate space (0.7*Maxent current) + (0.3* Bioclim current). B. Weighted future (2080) suitable climate area (0.7*Maxent future) + (0.3* Bioclim future).

Figure 4. 

Suitability index based on weighted average ensemble of Bioclim and Maxent for present and future climate areas (0.7 * Maxent) + (0.3 * Bioclim) shows that currently 23% of the area are suitable and the remaining 77% unsuitable for Daboia russelii in Bangladesh. In the future, 79% of the country area is suspected to be climatically unsuitable.

Currently, the estimated suitable climate space (Fig. 3A, B) is nearly five times larger than the species’ observed area of occurrence (AOO: 14,344 km²). These areas are predominantly located along major river systems, including the Padma, Meghna, and Jamuna, and encompass regular floodplains, low tidal plains, and diverse forest types such as tropical deciduous, mangrove, and village forests.

The current suitable climate space of D. russelii in Bangladesh is confined to the western half of the country. According to the weighted ensemble of predicted models (Maxent and Bioclim), the major portion of the suitable space is situated at the banks and char lands of the Padma, Meghna, and the Jamuna River. Most of these areas are riverbanks, and lower flat plains, and cover the major portion of the village forest of the southern, central, north-west and northern part of the country.

The river bank area of the north-west and central parts is the most suitable place for this viper (p > 0.6) under the current climate scenario (Fig. 3A), but its suitability is slightly reduced in 2080 (p < 0.4) (Fig. 3B). The southern coastal area and the mangrove forest also provide current suitable space (0.4 > p < 0.6), but reduced in the future (p < 0.2) (Fig. 3A, B).

Both Bioclim and Maxent models strongly support that the riverbank, char lands, and lower floodplains of Padma, Meghna and Jamuna rivers of the country provide the current suitable climatic condition for D. russelii in Bangladesh (p > 0.8 for Maxent, 0.5 for Bioclim). In the future (2080), this suitable climate situation slightly decreases (p ≤ 0.005 for Bioclim and p ≤ 0.5 for Maxent) (see details in Suppl. material 5).

The Bioclim model proved to be highly effective in the current climatic condition, with high discriminatory power (AUC = 0.9851), in addition to high sensitivity (0.9821) and specificity (0.9739), thus indicating its ability to correctly differentiate between suitable and unsuitable climate space. However, its performance significantly declined when tested against future climate scenarios (AUC = 0.6512), as indicated by the low specificity (0.3087) and low Kappa value (0.3054) (Table 2).

In contrast, the Maxent model showed consistently strong performance in both the current and future contexts, as evidenced by consistent AUC values (0.9730 for the present context and 0.9718 for the future context), maximum sensitivity (1.0), and high specificity (0.9522) in both cases. This highlights the ability of Maxent to clarify complex environmental interactions and its resilience to changing climatic conditions, thus making it a more reliable tool for predicting the distribution of D. russelii in future contexts (Table 2).

The SDM analysis for D. russelii in Bangladesh shows that 23% of the country’s land is climatically suitable, while 77% is unsuitable. Future projections reveal that 21% of the land will be suitable, compared to 79% that will remain unsuitable. The areas found to be suitable in the current and future evaluations include 19%, while 76% of the country is unsuitable in both temporal analyses. Additionally, about 4% of the currently suitable areas are expected to shift to an unsuitable category, while about 1% of the areas found to be unsuitable will be suitable in the future (Fig. 4).

Despite these changes, the regular floodplain remains a consistent feature within both present and future suitable climatic spaces, underscoring its importance as a stable habitat component for D. russelii in Bangladesh.

Principal Component Analysis (PCA) underscored the complex interplay of climatic and spatial factors driving the distribution and habitat suitability of D. russelii in Bangladesh. It identified three components (PC1, PC2, and PC3) that collectively explained over 79% of the variance in the occurrence of D. russelii. These components encompassed 15 significant variables (loading value > 0.7), which either positively or negatively influence the species distribution (Table 3).

PC1 accounted for the largest variance and highlighted key variables influencing the species distribution. Positive contributors included minimum temperature of the coldest month (0.937), annual precipitation (0.805), precipitation of the driest month (0.942), and precipitation of the driest quarter (0.905). Conversely, variables with negative influences included mean diurnal range (-0.941), temperature seasonality (-0.920), temperature annual range (-0.868), and elevation (-0.724). PC2 identified additional variables affecting the distribution of D. russelii. Positively associated factors were annual mean temperature (0.725) and maximum temperature of the warmest month (0.722). In contrast, negative influences came from isothermality (-0.737) and precipitation of the warmest quarter (-0.935).

Over a four-year study period, the collected individuals of D. russelii exhibited a distinct female-biased trend (Fig. 5). The trend lines in Fig. 5 further support this observation, showing an upward trajectory for females and a corresponding decline in males. The high coefficient of determination (R² = 0.98) for both male and female trends suggest a strong, consistent pattern in the sex ratio dynamics over the analysed period. In 2019, the sex ratio of collected specimens was balanced at 1:1, indicating equal male and female representation. By 2020, females began to outnumber males with a ratio of 1:1.3, reflecting an increase in the female proportion to approximately 56.5% of the sample. This female dominance became more pronounced in 2021, with females accounting for 66.7% of the population and a sex ratio of 1:2 favouring females.

Figure 5. 

Sex-based population structure of Daboia russelii collected in Bangladesh from 2019 to 2021 suggests a gradual increase in females in natural habitats.

The apparent population expansion of D. russelii in Bangladesh may also be influenced by a significant decline in both its predators and competitors. Specifically, the population of predators has declined by 29%, while competitors have seen a more substantial decline of 46% (Fig. 6A, B). Notably, the status of many species remains uncertain, with 36% of predators and 32% of competitors categorized as “unknown” according to IUCN (2015).

Figure 6. 

Decreasing population trends of potential predators (A) and competitors (B) might facilitate population expansion of Daboia russelii. A total of 32% of potential predators (C) and 30% of potential competitors (D) were assessed under threatened or near threatened categories regionally by IUCN Bangladesh in 2015. Lists of potential predators and competitors of D. russelii with their population trends and status are provided in Supplementary material 3 and 4, respectively. E. Increasing trends of land use for rice cultivation and gross production of rice over the years in Bangladesh (2001-02 to 2021-22) provided additional paddy field habitat for D. russelii and more food for their rodent prey. F. Flood affected areas in Bangladesh over the last 20 years (2001-2020) might provide dispersal pathways for D. russelii towards low plains.

Conversely, only a small proportion of species—4% of predators and 1% of competitors—have shown an increasing population trend. Threat levels are also concerning, with approximately 14% of predators and 11% of competitors regionally classified as threatened or near-threatened (CR, EN, VU) by IUCN Bangladesh (2015) (Fig. 6C, D). Additionally, 18% of predators and 19% of competitors are classified as Near Threatened (NT) due to declining population trends.

The analysis of habitat selection based on land use indicates a pronounced trend toward increasing rice cultivation over the past two decades. Specifically, gross rice production has shown a rapid increase with a strong correlation (R² = 0.93) between time and production trends (Fig. 6E). Between the fiscal years 2001-02 and 2021-22, the cultivated area dedicated to rice expanded modestly by 0.49%, while total rice production surged by 2.32%. This reflects an intensification of agricultural practices, likely driven by higher yields per year rather than a significant expansion in cultivated area. Greater and more frequent availability of rice, for example, due to multiple planting and harvest cycles per year, is expected to result in higher rodent populations which in turn could sustain higher populations of rodent-eating snakes like D. russelii.

Over the last two decades, Bangladesh has shown an increasing trend in flood-affected lowland areas (Fig. 6F). Notably, major flood events in 2004, 2007, 2017, and 2020 impacted 38%, 42%, 42%, and 40% of the country’s total area, respectively. These data indicate a consistent rise in the extent of flood-affected regions. Floods can facilitate the terrestrial species over the inundated landscape.