Introduction
Environmental factors play a major role in determining the distribution of species (Dias-Terceiro et al. Reference Dias-Terceiro, Kaefer, de Fraga, de Araújo, Simões and Lima2015, Vellend Reference Vellend2010). Being poikilothermic animals, snakes are especially sensitive to environmental gradients, which also makes them highly susceptible to climate change (Needleman et al. Reference Needleman, Neylan and Erickson2018). In this regard, Bangladesh is expected to be particularly affected in the near future, not only via changes in temperature and precipitation but also by sea-level rise (Shahid Reference Shahid2011). This might threaten the status of snakes in this tropical country. Therefore, studying the effects of climate change on the distribution of snakes in Bangladesh is essential for designing effective conservation measures.
At large scales, the amount of favourable habitat of tropical snakes may be reduced due to climate change (Lourenço-de-Moraes et al. Reference Lourenço-de-Moraes, Lansac-Toha, Schwind, Arrieira, Rosa, Terribile, Lemes, Fernando Rangel, Diniz-Filho, Bastos and Bailly2019). Studies on local scale have shown that snakes are highly impacted by any alteration in their ecosystems and by changes in forest structure and land use (Breininger et al. Reference Breininger, Mazerolle, Bolt, Legare, Drese and Hines2012, dos Santos et al. Reference dos Santos, de Oliveira, Gonçalves, de Almeida, Loebmann and Tozetti2013, Rahman et al. Reference Rahman, Rashid, Das and Luiselli2013b). The distribution of the ecosystem shapes the occurrence range of tropical squamates (Reed Reference Reed2003). The spatial inaccessibility of favourable habitats due to wide rivers and regular floods also impacts the range expansion of snakes living at low altitudes (Moraes et al. Reference Moraes, Pavan, Barros and Ribas2016). Moreover, elevation gradients have a direct influence on snake diversity (Colwell et al. Reference Colwell, Brehm, Cardelús, Gilman, Longino, Cardelus, Gilman and Longino2008, Fu et al. Reference Fu, Wang, Pu, Zhang, Chen, Zhao, Chen and Wu2007). All these factors are active within the territory of Bangladesh and are assumed to highly affect the distribution of terrestrial snakes. Being situated in the Indo-Himalayan region, Bangladesh has seven different climatic sub zones, which are characterised by seasonal variations of temperature, rainfall, humidity, fog-mist-dew, winds, cyclones, tornadoes and atmospheric pressure (Rashid Reference Rashid1991). Twelve different bio-ecological sub zones have been officially identified for Bangladesh according to physiography, soil composition, temperature, rainfall and floral and faunal diversity (Nishat et al. Reference Nishat, Huq, Barua, Reza and Khan2002). The country also has five distinct types of forest. For instance, the largest mangrove ecosystem in the world, the Sundarbans, located in the south-western section (Iftekhar Reference Iftekhar2006). This forest hosts a unique diversity of plants and animals (Aziz & Paul Reference Aziz and Paul2015, Iftekhar Reference Iftekhar2006, Sarker et al. Reference Sarker, Reeve, Paul and Matthiopoulos2019). Tropical hill forest occurs in hilly areas of the eastern part of the country. The top canopy trees of this forest reach a height of 25–62 m and about 800 species of flowering plants are found in this forest (Iftekhar Reference Iftekhar2006). This forest host the highest herpetofaunal diversity among forest types (Payo et al. Reference Payo, Mukhopadhyay, Hazra, Ghosh, Ghosh, Brown, Nicholls, Bricheno, Wolf, Kay, Lázár and Haque2016). Tropical deciduous forest, commonly known as ‘Sal’ (Sohrea robusta) forests, are distributed in the central part of the country. The tropical evergreen forest occurs in the north-western part of the country and the village forest made up of a variety of plantations throughout the country supporting diverse wildlife (Mukul et al. Reference Mukul, Biswas, Manzoor Rashid and Pullaiah2018) (Supplementary Figure 1e). All these forest are facing high rate of degradation due to change in climate and socio-economic condition (Mukul et al. Reference Mukul, Biswas, Manzoor Rashid and Pullaiah2018) that may cause decrease in suitable habitats for snakes. The diversity of forest in Bangladesh, but also other particularities of this country, creates a complex and patchy habitats for snakes. For instance, around 400 rivers originating mostly from the Himalayas cross the country and spill into the Bay of Bengal in the south of Bangladesh (Banglapedia 2021). Three of these rivers are more than a kilometre wide, with a strong current, creating natural barriers for terrestrial animals all the way from their origin upstream down to the estuary (Banglapedia 2021). Heavy rainfall during the Monsoon causes seasonal floods in wide floodplains alongside the rivers (Deo et al. Reference Deo, Adamowski, Begum, Salcedo-Sanz, Kim, Dayal and Byun2019).
Bangladesh has a rich snake diversity including 83 terrestrial snake species belonging to 11 families (IUCN Bangladesh 2015). In the latest Red List assessment, eight of these are assessed as threatened and 32 are either not evaluated or data deficient (IUCN Bangladesh 2015). Regular and updated editions of field-based publications on the occurrence and distribution of these species have provided details of their occurrence within the country (Khan Reference Khan2015, Reference Khan2018). However, an ecological understanding of the distribution of snakes in Bangladesh is limited, hampering conservation efforts. Like in every tropical country, snakebite is a health burden, with about a million incidents in rural Bangladesh and 6,000 deaths per year (Rahman et al. Reference Rahman, Faiz, Selim, Rahman, Basher, Jones, d’Este, Hossain, Islam, Ahmed and Milton2010). A total of 33 terrestrial snake species in the country use venom to subdue their prey. Thirteen of these venomous species (Bungarus caeruleus, Bungarus fasciatus, Bungarus lividus, Bungarus niger, Bungarus walli, Naja kaouthia, Naja, Ophiophagus hannah, Daboia russelii, Trimeresurus albolabris, Trimeresurus erythrurus, Trimeresurus gramineus and Trimeresurus popeiorum) have a large impact on public health because of their venom toxicity, whereas the remaining ones are considered non-lethal (e.g. mild or no toxicity symptoms) to humans (Rahman et al. Reference Rahman, Faiz, Selim, Rahman, Basher, Jones, d’Este, Hossain, Islam, Ahmed and Milton2010). Designing an appropriate antivenom, considering interspecific venom variations and assuring its sufficient supply in local health services is crucial for effective management of snakebites (Gutiérrez Reference Gutiérrez2014), thereby also considering an intraspecific variation of venom (Laxme et al. Reference Laxme, Attarde, Khochare, Suranse, Martin, Casewell, Whitaker, Sunagar, Senji Laxme, Attarde, Khochare, Suranse, Martin, Casewell, Whitaker and Sunagar2021). Transportation of a patient to a hospital is often time consuming, and the duration between bite and hospital reach is found to be inversely related to a patient’s survival probability (Mondal et al. Reference Mondal, Chowdhury, Rani, Mohammad, Islam, Haque and Faiz2012). Thus, an occurrence map of venomous snakes along with accurate predictions of suitable habitats may guide the distribution of resources in affected areas (Yañez-Arenas et al. Reference Yañez-Arenas, Peterson, Mokondoko, Rojas-Soto and Martínez-Meyer2014), but either approach relies on an accurate prediction of the distribution of venomous snakes and their preferred habitats.
This study aimed to describe terrestrial venomous snake occurrences in Bangladesh using ecological niche modelling. We analysed the present distribution of the terrestrial venomous snakes using geographical, climatic and habitat parameters as well as by projecting their future distributions. The study was also to find out the limiting factors for the distribution of the venomous snake in this tropical country. Based on our results, we proposed guidelines for the conservation of venomous snake diversity and snakebite management in Bangladesh.
Methods
Study area
Our analyses were focused on Bangladesh, a tropical country situated at 23°41' N, 90°20' E and with an area of 147,570 km2. It is bordered by India to the east, west and north and by the Bay of Bengal to the south. Bangladesh exhibits highly heterogeneous environmental conditions, provided by a wide range of climatic sub zones and diverse vegetation types (Rashid Reference Rashid1991) (Figure 1 and Supplementary Figure 1).
Geoposition of Bangladesh in Indo-Himalayan and Indo-Chinese realm (above). Occupancy of venomous snakes in Bangladesh in respect to seven spatial variables (below), (a) 12 different ecosystems, (b) seven climatic sub-zones, (c) five forest systems, (d) various elevation levels, (e) 11 types of land use, (f) inside and outside the regular flood zone and (g) three zones separated by wide rivers. Categories of these variables are given in Table 1.
Data collection on species distribution
We compiled occurrence information for 29 of the 33 terrestrial venomous snakes occurring within Bangladesh. The occurrence map of each snake was constructed by compiling the occurrence records of each species from the published literature on snake occurrence and distribution in Bangladesh (Hasan et al. Reference Hasan, Khan and Feeroz2014, IUCN Bangladesh 2015, Khan Reference Khan2015, Reference Khan2018) and we digitised the polygons using ArcGIS v10.4.1. Therefore, occurrence vector file of each species was consisted of small polygons that presented the patches of small habitat across the occurrence range of the particular species. These spatial polygon files were used as the current range (CR) of the studied species. Four species were excluded from the analysis because of insufficient occurrence information to generate their CR maps.
Variable selection
Seven environmental variables (Table 1) were collected from different national and international data repositories. The bio-ecological sub zones (hereafter denoted as ecosystem) and the climatic sub zones of Bangladesh were digitised from published maps (Nishat et al. Reference Nishat, Huq, Barua, Reza and Khan2002, Rashid Reference Rashid1991) and then rasterised. The regular flood zone (flood) and parts of the country which are segmented by wide river (river) spatial data were collected from the Bangladesh Agricultural Research Council and Water Development Board then rasterised (Bangladesh Agricultural Research Council 2019, Bangladesh Water Development Board 2010). These four raster files were in one km2 pixel resolution. The SRTM 3 arc-second (90 × 90 m2) land surface elevation (elevation) data were extracted from the Earth Resources Observation and Science Center (United States Geological Survey 2019), which was further resampled to 1 km2 spatial resolution using the nearest neighbour resampling technique in ArcGIS. The latest land use/land cover (landuse) raster was in 100 × 100 m2 pixel resolution and the updated forest cover raster was in 30 × 30 m2 pixel resolution, both were collected from the Copernicus global land service, which is widely used for sustainable development goals such as conservation (Buchhorn et al. Reference Buchhorn, Smets, Bertels, Lesiv, Tsendbazar, Herold and Fritz2019), and the Bangladesh Forest Department (Bangladesh Forest Department 2019), respectively. These last two rasters were used with comparatively higher spatial resolutions because there were several small size land use and forest patches having areas of less than one km2. The values of the raster cell were categorical. For example, there are only two categorical cells in raster, outside flood zone (cell value = 1) and inside flood zone (cell value = 2) (Table 1). The other six environmental variables, i.e. ecosystem, climatic sub zone, elevation, forest, land use and river, have twelve, seven, ten, five, eleven and three categories, respectively, and they were assigned to different values in each raster file (Table 1 and Supplementary Figure 1).
Seven independent variables and their categories used for the current analysis (described in methods)
Nineteen bioclimatic variable raster files (1 km2 resolution) were downloaded from Worldclim for the present and the future (prediction for the year 2070) climate (www.worldclim.org accessed on 27.09.2019) (Hijmans et al. Reference Hijmans, Cameron, Parra, Jones and Jarvis2005). The cell values in bioclimatic raster files were derived or calculated through the different measuring units and so they were standardised (z-scored). To minimise the bias and collinearity of variables, between highly correlated variables (r > 0.9) the one which is important to species was selected (Brun et al. Reference Brun, Thuiller, Chauvier, Pellissier, Wüest, Wang and Zimmermann2020, De Marco & Nóbrega Reference De Marco and Nóbrega2018) and finally 11 bioclimatic variables (bio2, bio3, bio4, bio7, bio10, bio11, bio14, bio15, bio16, bio17 and bio18) were selected for analysis. For future (the year 2070) climate prediction, we chose a moderate approach, i.e. the Representative Concentration Pathway 4.5 (RCP 4.5), for a future global warming scenario. All vectors and raster files were developed in the WGS84 coordination system and masked with the Bangladesh border.
Occupancy calculation
Using the spatial polygon (CR) as a boundary of patchy habitats, the frequency of pixels in the different categories of a raster file (independent variable) was obtained by using the extract() function of the ‘raster’ package (Hijmans et al. Reference Hijmans, Van Etten, Sumner, Cheng, Bevan, Bivand, Busetto, Canty, Forrest, Ghosh, Golicher, Gray, Greenberg, Hiemstra, Karney, Mattiuzzi, Mosher, Nowosad, Pebesma, Lamigueiro, Racine, Rowlingson, Shortridge and Wueest2019). It was done separately on each environmental raster because of the aforementioned variation in the resolution of the raster files. The number of pixels was then converted to square kilometre unit to get the size of the segments of a habitat by multiplying pixel values and resolution of the particular raster. Then, the percentage proportion of CR area over different categories of a variable was calculated. The presumption is that the bigger the portion of CR is spread over a category, the more preferable the category is to the species. An average value of this proportional distribution of 29 species over a category was calculated and set as an occupancy score for that category. For example, we averaged the size of 29 CRs over village forest and assigned that value to all pixels of village forest in the forest raster file by using the reclassify() function of the ‘raster’ package. In the same way, average scores were calculated and assigned to all 50 categories in the seven raster files.
Species distribution models
Fifty occurrence coordination points were randomly selected from the CR polygons of each snake by using the bas.polygon() function of the ‘Balanced Acceptance Sampling package’ (Mcdonald Reference Mcdonald2016). This function prevents the repetition of points from the same pixel. It is to note that each of 29 CR is consisted of multiple small polygons presenting occurrence of each of 29 venomous snakes. These coordination points were then used to extract values of 11 bioclimatic and seven aforementioned environmental raster files. The bioclim() function of the ‘Dismo’ package (Hijmans et al. Reference Hijmans, Phillips, Leathwick and Elith2017) was trained using those values and the predict() function was used to produce potential niches for each venomous snake in the study area. We designed a climatic niche by using only bioclimatic variables and an environmental niche by using seven spatial variables and a combined niche model by using all 18 variables.
Ecological niche models (ENMs) are often used to predict the effects of climatic change on the distribution of a species (Rangel & Loyola Reference Rangel and Loyola2012). Bioclim algorithms applied on the climatic variables of the presence-only data of a species to generate a projection of the climatically suitable areas (Beaumont et al. Reference Beaumont, Hughes and Poulsen2005). This algorithm use percentile distribution of the data in the n-dimensional climatic space and thus the areas inside the climatic tolerance sampled for the species (Hijmans et al. Reference Hijmans, Phillips, Leathwick and Elith2017).
In each CR, the pixels of the small polygons were set to one to show the favourable pixel for the respective snake species and zero for unfavourable. The clustering of such 29 CRs produced a species richness (SR) map that is based on field knowledge and publications from Bangladesh (see data collection on species distribution section). The cell values (from 0 to 29) of the SR map presented the number of snake species using that pixel as their potential habitat. Similarly, the pixel values of potential climatic niche of each study species were set to either 0 or 1 to present unfavourable or favourable pixel for that species, respectively. Again, these 29 potential climate niches are clustered into one raster where the pixel values ranged from 0 to 29. For example, if a pixel value is 10, it meant that pixel was favourable to 10 study species. Similarly, we clustered the current environmental niches, the combined niche (combination of current climatic niche and current environmental niche) and the future climatic niches of the 29 study species.
Model analysis
The resemblance between the SR map derived from the CRs and the three projected niches of the studied species were tested in two ways. First, the Pearson correlation was employed to test the significance of relationships between frequencies of pixels supporting the same number of species in two raster files. Second, the values of pixel from identical locations in the current SR maps and in any of the three projected niche maps were compared. This pixel-level comparison is expected to reveal the degree of accuracy of the prediction comparing the current occurrences. All these statistical operations were conducted on the R platform using the packages ‘raster’ (Hijmans et al. Reference Hijmans, Van Etten, Sumner, Cheng, Bevan, Bivand, Busetto, Canty, Forrest, Ghosh, Golicher, Gray, Greenberg, Hiemstra, Karney, Mattiuzzi, Mosher, Nowosad, Pebesma, Lamigueiro, Racine, Rowlingson, Shortridge and Wueest2019), ‘maptool’ (Bivand et al. Reference Bivand, Lewin-Koh, Pebesma, Archer, Baddeley, Bearman, Bibiko, Brey, Callahan, Carrillo, Dray, Forrest, Friendly, Giraudoux, Golicher, Rubio, Hausmann, Hufthammer, Jagger, Johnson, Luque, MacQueen, Niccolai, Pebesma, Lamigueiro, Plunkett, Short, Snow, Stabler, Stokely and Turner2019b), ‘rgdal’ (Bivand et al. Reference Bivand, Keitt, Rowlingson, Pebesma, Sumner, Hijmans, Rouault, Warmerdam, Ooms and Rundel2019a), ‘sp’ (Pebesma et al. Reference Pebesma, Bivand, Rowlingson, Gomez-Rubio, Hijmans, Sumner, MacQueen, Lemon, O’Brien and O’Rourke2019) and ‘sf’ (Pebesma et al. Reference Pebesma, Bivand, Racine, Sumner, Cook, Keitt, Lovelace, Wickham, Ooms, Müller, Pedersen and Baston2020).
Results
The extent of the CR of the five countrywide distributed species is almost as large as the area of the country, while five further species have a less than 1,000 km2 area of occupancy within Bangladesh. Except for these 10 species, the CR size of the other studied species varies from 5,722 to 47,672 km2.
Characterising habitat
Out of the 50 categorical states of the seven variables (details in Table 1), 35 of them are either absent or present in a negligible portion in a CR while 15 other categories are found in a wide part of CR of the studied species. Of these 15 contributory categories, 11 can partially characterise the habitat of multiple studied species while each of the remaining four parameters partially contributes to the habitat of a single species (Supplementary Table 1). An average of 88% (58-100%) area of the CR of the study species stretched outside the regular flood zone of the country, while the portion of CR distributed in central Bangladesh was mostly situated inside the regular flood zone (Figures 1f, 2f) and in the flood plain ecosystem (Figure 2a). Of the 12 types of ecosystem distributed in Bangladesh, the hilly ecosystem in the eastern side of the country constitutes the major portion (average ∼80%) of the CR of a group of 16 species (Figure 1a; category 9 in Figure 2a). They are distributed over two slightly different climatic sub zones, i.e. the south-eastern climatic sub zone characterises the habitat of (>90% area of CR) 12 venomous snakes of this group while the north-eastern sub zone is favoured (>83% area of CR) by the four other snake species. Altogether, the habitat of this group of snakes is shaped by a set of seven parameters (Supplementary Table 1).
Area of current ranges (in %) of 29 venomous snakes distributed in different categories of seven environmental variables. (a) 12 types of ecosystems, (b) seven climatic zones, (c) five forest types, (d) 10 categoric elevations, (e) 11 types of land use, (f) regular flood areas and (g) three zones separated by wide rivers. “N” is the number of species. Categories of these variables are given in Table 1.
Village forest is the next preferred habitat component (average 86% area of CR) to another group of 11 venomous snake species (category 1 in Figure 2c). Except for Ahaetulla nasuta, all these snakes are surface dwellers, nocturnal and pose a significant threat to public health, including two species of cobras, four kraits and one viper. Sixty-six percent of CR of these snakes is also associated with the regular cropland, a category of land use variable. The other 10 types of land use occupied a small portion of the CR (∼7%) of the study species. On the other hand, the occurrence of venomous snakes in Bangladesh gradually decreases in higher elevations (Figure 2d). In fact, of the 29 studied species, 27 have their habitat (average 81% area CR) within a 100-meter vertical elevation. The remaining two species, the Pope’s Pit Viper (Trimeresurus popeiorum) and Mountain Pit Viper (Ovophis monticola), have a large portion of their habitat at higher elevations.
The terrestrial habitats of Bangladesh are separated by the Padma-Meghna-Jamuna river flow, which plays a major role in restricting the spread of 22 of the studied snake species. Eighteen of them have their major portion (88.84%) of habitat in the eastern part of the country, two occur mostly in the south-western part and two snakes in the north-eastern region (Figure 2g). However, some species show a wide distribution, having a range spanning over multiple categories of variables. For example, seven species, including five countrywide distributed taxa, have a wide CR in all three segments of the country separated by the river. Similarly, the habitats of 12 of the studied snake species extend over multiple ecosystems and climatic sub zones in Bangladesh. Two species are found in all types of forest and 16 occurred in each category of the land-use variable. However, species belonging to the same family did not necessarily belong to similar environmental and climatic variables (Supplementary Figure 2).
Climatic niche
The geographical area of favourable climatic niches and the CR of each species are not identical, and neither is the whole CR inside the potential climate surface nor is the whole climatic niche inside the CR. On average, 85% (75.49–100%) of the area of CR are situated inside their potential climatic niche and the rest is outside of the favourable climate surface. Although the climatic niches are always larger than the CR for each species, the size of the 29 CRs and climatic niches have a highly significant association (r = 0.88, p = 0.000). The predicted present climatic niches of five countrywide distributed venomous snakes are as large as the area of the country (Table 2), while the climatic niche of other studied species ranges from 59 to 125,026 km2. The cumulative climatic niche of the studied species shows a greater assemblage of snakes in the southwest, southeast and northeast regions of the country, which is similar to the CRs of the venomous snake community (Figure 3a–d).
The size of the current range (CR) and predicted climatic niches of 29 venomous snake species of Bangladesh
▴ = Distributed throughout Bangladesh, DD = Data Deficient, EN = Endangered, LC = Least Concern, NT = Near Threatened, VU = Vulnerable.
Potential habitat (frequency of pixels) supporting various numbers of terrestrial venomous snakes (cell value = number snakes may share the 1 km2 cell) in Bangladesh, (a) species richness based on current occurrence range, (b) present potential climatic niche, (c) present environmental niche, (d) present potential niche (a combined model) and (e) future potential climatic niche.
Future climatic niche
We estimated that the potential climatic space of each species will decrease dramatically and shift to the northern part of the country in the next half-century (Figure 3e). Output of the bioclim model suggested that around 80% area of the country will have no potential climatic surface in the near future for any venomous snake currently occurring in the country. This future climatic niche may support less than five species, and a negligible area in the eastern part of the country may support as many as 12 of the study species (Figure 3e). Except for the countrywide distributed snakes, each venomous snake species will lose an average of 99% of its present climatic niche individually within the next 50 years (Table 2). Eleven of the study species, i.e. Amphiesma platyceps, Boiga siamensis, Chrysopelea ornata, Pseudoxenodon macrops, Rhabdophis himalayanus, Rhabdophis subminiatus, Bungarus lividus, Ophiophagus hannah, Daboia russelii, Ovophis monticola and Trimeresurus popeiorum will have almost no potential climatic niche (<10 km2) within the country, while 12 further species will have less than 1000 km2 of favourable climate space (Table 2).
Model analysis
The pattern of the SR maps generated from the three models (climate, environment and combined models) and CRs are apparently similar (Figure 4). Despite these similarities, the frequency of pixels supporting a different number of venomous snake species in the climatic and environmental niches is not correlated to that of the CR map (Figure 3a). Niches predicted by using all 18 studied variables are correlated (r = 0.44, p < 0.05) to the CR of the study species which was sourced from the publication of the field experts (IUCN Bangladesh 2015, Khan Reference Khan2015, Reference Khan2018) (Figures 3, 4). Around two percent of the pixels located in the identical position in the combined niche and the CR support the same number of snake species, whereas 73% of pixels in the combined niche may support more species than an identically positioned pixel in the CR (Figure 4).
(a) Frequency of pixels (size 1 km2 each) supporting a different number of snake species and (b) identically located pixels in the combined niche support different snake species richness than those of current occurrence map.
Discussion
Our results indicate that the regular flood zone located in the central part of the country might be the reason for having a large assemblage of venomous snakes in the southeast, southwest and northeast corners of the country. Floods force burrowers and surface dwellers to higher land for survival, which regulates seasonal prey dynamics (Madsen & Shine Reference Madsen and Shine1996, Piatti et al. Reference Piatti, Rosauer, Nogueira, Strussmann, Lúcia and Martins2019) and increases accidental encounters between snakes and humans (Tauzer et al. Reference Tauzer, Borbor-Cordova, Mendoza, De La Cuadra, Cunalata and Stewart-Ibarra2019) during the Monsoon in rural Bangladesh (Rahman et al. Reference Rahman, Faiz, Selim, Rahman, Basher, Jones, d’Este, Hossain, Islam, Ahmed and Milton2010). Also, the regular occurrence of floods is one of the drivers in the variation in reptiles species diversity in the tropics (Moraes et al. Reference Moraes, Pavan, Barros and Ribas2016). In Bangladesh, the catchment area alongside the major rivers inundates up to 25% to 70% of land every year (Dewan et al. Reference Dewan, Nishigaki and Komatsu2003, World Meteorological Organization 2003). In addition, all major rivers of Bangladesh are shared with other countries situated upstream to those rivers, and the floodwater in those countries is carried down to Bangladesh gradually. The extensive flood plains at sea level cause slow movement of water, which means floods linger from weeks to months during the Monsoon season in central Bangladesh, causing unsuitable habitats for terrestrial snakes. In addition, there is an increasing trend in annual and pre-monsoon rainfall over the last 50 years, which may increase the degree of flooding in the future (Khalequzzaman Reference Khalequzzaman1994, Shahid Reference Shahid2011). Under this predicted worsening of floods, habitat disruption by flooding will be an increasing threat to the reptile communities in Bangladesh. As most of the studied snake species breed during the rainy season (IUCN Bangladesh 2015), the co-occurrence of the breeding season with the flood during the Monsoon might disrupt the species’ reproductive success and may result in poor population structure (Rahman et al. Reference Rahman, Rashid, Das, Jenkins and Luiselli2013a) (Table 3).
Habits and habitats of 29 studied venomous snakes of Bangladesh
Ar = Arboreal, S = Surface dweller, and F = Fossorial, M = Small mammals, B = Small birds, R = Reptiles, A = Amphibians, F = Fishes, E = Eggs, and I = Insects, - = No information.
Moreover, a negligible portion of the CR (around one-tenth) of 27 species is located within the regular flood zone. This is because of the occurrence of favourable habitat in central Bangladesh (Figure 2f), which in turn becomes inaccessible for terrestrial species during floods (Figure 1f). For example, N. naja and N. kaouthia are repeatedly observed taking shelter in human habitations during the flood though they are still well distributed in the flood plain (IUCN Bangladesh 2015). To counteract flooding, several plans were executed in Bangladesh (Dewan et al. Reference Dewan, Nishigaki and Komatsu2003, Khalequzzaman Reference Khalequzzaman1994), mostly focused on minimising the impact of flooding on the economy and humans, but not on biodiversity (Deo et al. Reference Deo, Adamowski, Begum, Salcedo-Sanz, Kim, Dayal and Byun2019). However, the recently initiated delta plan has adopted sustainable development goals, which are expected to address this impact (Planning Commission Bangladesh 2018). Targeted actions on flood management may reduce the threat of habitat loss due to the worsening of flooding in the future (Sexton et al. Reference Sexton, Drda, Sexton and Bramble2007) and might positively assist terrestrial snake populations in expanding the limits of their favourable niche.
The wide and strong flows of three major rivers – the Padma, the Meghna and the Jamuna – divide the country into three distinct segments (Figure 1g). Our study found that a large portion of the study species is widely dispersed across the eastern part of the country, whereas their distribution to the west is limited by the Padma-Meghna-Jamuna flow. In Bangladesh, the extent of Daboia russelii is restricted by these rivers to the west part of the country (Wüster et al. Reference Wüster, Otsuka, Malhotra and Thorpe1992). Strong flows of wide rivers have been documented to restrict the dispersion of terrestrial snakes and also impede gene flow over evolutionary time scales (Gibbs et al. Reference Gibbs, Sovic, Amazonas, Chalkidis, Salazar-Valenzuela and Moura-Da-Silva2018, Pla et al. Reference Pla, Sanz, Quesada-Bernat, Villalta, Baal, Chowdhury, León, Gutiérrez, Kuch and Calvete2019).
The valleys and plains elevated above flood level in this hilly area in the eastern part of the country provide favourable habitat for a group consisting of 16 venomous snake species, and most of them are arboreal while some are surface dwellers (Table 3). This hilly terrain in the eastern part of the country is made up of hilly forests and ecosystems, diverse climatic sub zones and situated outside the flood plain. Two major microclimatic sub zones dominate this area (Mia et al. Reference Mia, Ali and Roy2016, Rashid Reference Rashid1991); in the north-eastern corner, the mean maximum temperature is rarely above 32°C and the mean minimum is <10°C, while the temperature in the south-eastern corner rarely exceeds a mean of 32°C and goes below a mean of 13°C (Rashid Reference Rashid1991). The average humidity is higher in the north than in the south, and the northeast is the cloudiest part in Bangladesh while the south-eastern climate is characterised by heavy rainfall usually over 2,540 mm (Banglapedia 2012). To conserve biodiversity, this tropical and sub-tropical hilly forest dominated with tropical evergreen plants, dense canopy, diverse ecosystems and the largest watershed in the eastern part of the country are currently managed by more than a dozen protected areas (Ahsan et al. Reference Ahsan, Aziz and Morshed2016). However, as per our future climatic niche projection, the suitability of those protected areas currently located in the south-eastern and south-western part of Bangladesh might be reduced (Figure 3e). Studies have suggested that a higher altitude with lower temperatures in abundantly shaded environments with low evapotranspiration might provide favourable conditions and serve as thermal refuges or retreat sites for snakes in future global warming (Colwell et al. Reference Colwell, Brehm, Cardelús, Gilman, Longino, Cardelus, Gilman and Longino2008, Fu et al. Reference Fu, Wang, Pu, Zhang, Chen, Zhao, Chen and Wu2007, Kearney et al. Reference Kearney, Shine and Porter2009). So to face climate change, a proper management plan to improve the habitat of this area might conserve the high richness of venomous snakes (Rashid Reference Rashid1991, Reed Reference Reed2003, Moraes et al. Reference Moraes, Pavan, Barros and Ribas2016) and other terrestrial vertebrates of Bangladesh (Colwell et al. Reference Colwell, Brehm, Cardelús, Gilman, Longino, Cardelus, Gilman and Longino2008, Lourenço-de-Moraes et al. Reference Lourenço-de-Moraes, Lansac-Toha, Schwind, Arrieira, Rosa, Terribile, Lemes, Fernando Rangel, Diniz-Filho, Bastos and Bailly2019). In addition, the world’s largest mangrove forest, Sundarbans, located in the south-western corner of the country will shrink due to climate change driven sea level rise (Payo et al. Reference Payo, Mukhopadhyay, Hazra, Ghosh, Ghosh, Brown, Nicholls, Bricheno, Wolf, Kay, Lázár and Haque2016). This demands a scientifically sound conservation strategies to protect suitable habitats of these venomous snakes from the impact of climate change (Raha et al. Reference Raha, Das, Banerjee and Mitra2012).
We also found that a group of 11 venomous snake species is distributed in village forests and cropland. Snakes of this group are also found in household patches of plantation and even inside the human residence that might further impact public health. These venomous snakes also have similar diets consisting of small mammals, birds and amphibians, which are abundant in cropland and village forests (Table 3). Both these habitat components are subject to rapid change and alteration; for example, there is a rotation of crops is a recent agricultural practice in Bangladesh that might impact prey availability (John et al. Reference John, Zaitsev and Wolters2021). As the change in land use is associated with the survival of species (Breininger et al. Reference Breininger, Mazerolle, Bolt, Legare, Drese and Hines2012), attention must be paid to sustainable land-use change for the conservation of biodiversity in Bangladesh. Moreover, the forest cover, including the village forests, in Bangladesh is less than 16% of the country area and is decreasing at a rate of 2.1% per year (Choudhury Reference Choudhury, Enters, Qiang and Leslie2003). It is predicted that this current forest cover may be reduced to half by 2050 (Choudhury Reference Choudhury, Enters, Qiang and Leslie2003). In addition, the recent influx of Rohinga refugees due to internal conflict in northern Myanmar has transformed 13 km2 of the forest into a human settlement, which is a small part of the continuous destruction of natural habitation and forest (Hassan et al. Reference Hassan, Smith, Walker, Rahman and Southworth2018). Thus, a sustainable multi-stakeholder-based conservation and forest extension process should be emphasised through social forestation and promoting public awareness (Iftekhar Reference Iftekhar2006). The detailed study of snake habitat in the seasonal cropland and village forests in Bangladesh may provide insights for conservation (Balouch et al. Reference Balouch, Rais, Hussain and Akram2016).
In a moderate climate scenario (RCP 4.5), the bioclim model predicted that the favourable climatic niche will be reduced to about one-fifth of the country size in the next 50 years (Figure 3b, e). Except for the north-central portion of the country, the rest of the country may not have any preferable climate surface for any snakes (Figure 3e). On an individual level, almost 99% of the potential present climatic niche of each studied species may be lost in a moderate global warming scenario. Our result showed that 11 snake species including the venomous Ophiophagus hannah and Daboia russelii will have no favourable climate surface (0–10 km2) in Bangladesh by 2070 (Table 2). These species are already on the threatened list for conservation due to habitat degradation and unstable population size (IUCN Bangladesh 2015). The favourable future climatic surface for the remaining snake species, including all five country-wide distributed species, may be reduced and restricted to the north-central part of the country. Such geographically shrunk and shifted potential niches were also predicted for terrestrial snakes, e.g., Bothrops spp., Agkistrodon contortrix, etc. (La Sorte & Jetz Reference La Sorte and Jetz2012, Needleman et al. Reference Needleman, Neylan and Erickson2018, Sahlean et al. Reference Sahlean, Gherghel, Papeş, Strugariu and Zamfirescu2014). Moreover, snakes have a low dispersal ability and are sensitive to environmental change, so the shrinkage of their potential niches may overwhelm their phenotypic plasticity, resulting in the extinction of those species in Bangladesh (Chevin et al. Reference Chevin, Lande and Mace2010, Maclean & Wilson Reference Maclean and Wilson2011, Segura et al. Reference Segura, Feriche, Pleguezuelos and Santos2007, Walther et al. Reference Walther, Post, Convey, Menzel, Parmesan, Beebee, Fromentin, Hoegh-Guldberg and Bairlein2002). The quick change in habitat superseding the adaptability of the snake may explain why snakes are declining in different areas of the world (Gibbons et al. Reference Gibbons, Scott, Ryan, Buhlmann, Tuberville, Metts, Greene, Mills, Leiden, Poppy and Winne2000, Reading et al. Reference Reading, Luiselli, Akani, Bonnet, Amori, Ballouard, Filippi, Naulleau, Pearson and Rugiero2010).
As per our analysis, five species have less than 1,000 km2 CR. Two of these are already assessed as threatened for having small and fragmented habitats (IUCN Bangladesh 2015). For the same reason, a further three species may meet the criteria for the threatened category and hence a revised assessment is demanded (IUCN Standards and Petitions Committee 2019). Our study shows that the absence of observed occurrences of some taxa in the potential climatic niche might be due to unfavourable flood zones, forests and ecosystems and even climatic sub zones which may not be accessible. Similarly, the geographically shifted future climatic space over such aforementioned unfavourable habitats may limit the occurrence of the snakes greatly (Dias-Terceiro et al. Reference Dias-Terceiro, Kaefer, de Fraga, de Araújo, Simões and Lima2015, Sahlean et al. Reference Sahlean, Gherghel, Papeş, Strugariu and Zamfirescu2014). On the contrary, snake species have also occurred in a favourable habitat located outside the suitable climate space (Morelli et al. Reference Morelli, Daly, Dobrowski, Dulen, Ebersole, Jackson, Lundquist and Millar2016). This is an adaptation of species through behavioural tolerance to the climatic factors outside the favourable climatic niche (Putman & Clark Reference Putman and Clark2017). In a climate change scenario, this may help snakes to survive even in an unfavourable climate. For instance, some studied snake species, including five country-wide distributed taxa, extend over an area of various types of land use, forests and ecosystems. Such wide adaptability might withstand the impact of future environmental change (Chevin et al. Reference Chevin, Lande and Mace2010, Phillimore et al. Reference Phillimore, Had, Jones and Smithers2010). Moreover, the persistence of species in their original ranges will be dependent on the degree of physiological and phenotypic plasticity or evolutionary adaptation of each species under changing conditions, or a combination of both processes (Holt Reference Holt1990, Visser Reference Visser2008).
The model analysis of our study indicates that combining environmental variables with climatic variables produced a realistic prediction for snake distribution (Rödder et al. Reference Rödder, Lötters, Öz, Bogaerts, Eleftherakos and Veith2011), and it may effectively update conservation plans. Adding more environmental variables may improve these models, and further precision might be achieved by including the low dispersal ability of snakes (Sahlean et al. Reference Sahlean, Gherghel, Papeş, Strugariu and Zamfirescu2014). The spatial distribution of categories of any variable over the territory of Bangladesh differs from the pattern of others (Figure 1). In the case of venomous snake species, such diversified habitats may also cause intraspecific venom variation in different populations of a species (Strickland et al. Reference Strickland, Smith, Mason, Schield, Borja, Castañeda-Gaytán, Spencer, Smith, Trápaga, Bouzid, Campillo-García, Flores-Villela, Antonio-Rangel, Mackessy, Castoe, Rokyta and Parkinson2018) and maybe differences in diet (Barlow et al. Reference Barlow, Pook, Harrison and Wüster2009). This demands designing antivenom by considering the intraspecific venom variations. For example, a population of pit vipers found in the hilly forest on the eastern side of the country and another in the mangrove forest in the south-western corner of the country are not connected by a continuation of suitable habitat and due to wide river flow (IUCN Bangladesh 2015), and thus may differ in their venom composition. Analysis of venoms of the medically important species distributed in different habitats will provide further details of intraspecific venom variation, allowing antivenom design with better efficacy (Chang et al. Reference Chang, Olenzek and Duda2015).
Regarding public health, 13 species including four country-wide distributed species contribute to human envenomation in the rural part of the country. Four species of kraits (Bungarus spp.), two cobras (Naja spp.), Russell’s viper (Duboia russelii) and three Green Pit Vipers (Trimeresusus spp.) are included in this group. These venomous snakes are the most threatening issue to public health in villages where human households are poorly managed (Harrison et al. Reference Harrison, Hargreaves, Wagstaff, Faragher and Lalloo2009). Besides the lack of public awareness, insufficient treatment facilities cause more human mortalities due to snakebite in village areas than in other parts of a tropical country (Harrison et al. Reference Harrison, Hargreaves, Wagstaff, Faragher and Lalloo2009), and this was also demonstrated to be the case in Bangladesh (Rahman et al. Reference Rahman, Faiz, Selim, Rahman, Basher, Jones, d’Este, Hossain, Islam, Ahmed and Milton2010). Another group of 17 mildly venomous snake species have no or very rare bite incidence to humans or their venom potency cannot kill humans but is still hazardous. None of these venomous snakes has medical records of deadly bites, as they are distributed in the less populated hilly areas of the country. Minimising impingement on these hotspots zones will reduce incidental bites. Particularly, assuring snakebite treatment in villages and regular flood areas may succeed in reducing the human death toll due to snake envenomation (Ghose & Faiz Reference Ghose, Faiz, Gopalakrishnakone, Faiz, Fernando, Gnanathasan, Habib and Yang2015, Ochoa et al. Reference Ochoa, Bolon, Durso, Ruiz de Castañeda, Alcoba, Babo Martins, Chappuis and Ray2020, Rahman et al. Reference Rahman, Faiz, Selim, Rahman, Basher, Jones, d’Este, Hossain, Islam, Ahmed and Milton2010).
Conclusion
Venomous snake habitat in Bangladesh is highly modulated by environmental and climatic factors. The major portion of the habitat of the study species is located at a low elevation outside the flood area. By analysing the spatial features of their current range, we found two groups of venomous snakes consisting of 16 and 11 species, respectively, whose habitat can each be described by a single set of spatial features. The first group prefers the hilly tropical forest and its associated ecosystems. Twelve of these species prefer the southeastern climatic sub zone, while four are distributed in the northeastern climatic sub zone. The distribution of the first 16 species in the eastern side of the country is constrained by the combined flow of the Jamuna-Padma-Meghna rivers, while the second group is distributed in the village forests and cropland habitats but have no preference for any particular climatic sub zone and riverside. Including country-wide distributed snakes, the second group consists of deadly venomous snake species that mostly occur near human habitation. We thus suggest that investing resources in snakebite management in village areas may reduce snakebite mortality among low-income villagers. Our study also found that the same species is occurring in variety of habitats thus increase the chance of intraspecific venom variations. This should be considered for medically relevant snake species when designing antivenom. From a conservation point of view, reducing flood areas, protecting the hilly and village forests and managing seasonal cropland land use may benefit the terrestrial snake species – and wider terrestrial biodiversity – of Bangladesh. For this, the hilly ecosystem on the eastern side of the country may provide a climatic refuge for terrestrial snakes species. Conserving this area by declaring more protected areas may be a partial remedy for the negative impact of climate change on snake biodiversity.
Supplementary material
To view supplementary material for this article, please visit https://doi.org/10.1017/S0266467422000359
Acknowledgements
We thank to Deutscher Akademisher Austauschdienst (DAAD) for the scholarship granted to the first author (Personal reference no 91609101). We acknowledge Ulrich Kuch and Faysal Bibi for their valuable comments on the earlier versions of the manuscript, and we appreciate anonymous reviewers for their suggestions to improvise the English.
Financial support
None.
Conflict of interest
None.
Ethical statement
None.
