https://orcid.org/0000-0003-1089-2076
Figures
Abstract
The lack of information regarding biodiversity status hampers designing and implementing conservation strategies and achieving future targets. Northern Pakistan consists of a unique ecoregion mosaic which supports a myriad of environmental niches for anuran diversity in comparison to the deserts and xeric shrublands throughout the rest of the country. In order to study the niche suitability, species overlap and distribution patterns in Pakistan, we collected observational data for nine anuran species across several distinct ecoregions by surveying 87 randomly selected locations from 2016 to 2018 in Rawalpindi District and Islamabad Capital Territory. Our model showed that the precipitation of the warmest and coldest quarter, distance to rivers and vegetation were the greatest drivers of anuran distribution, expectedly indicating that the presence of humid forests and proximity to waterways greatly influences the habitable range of anurans in Pakistan. Sympatric overlap between species occurred at significantly higher density in tropical and subtropical coniferous forests than in other ecoregion types. We found species such as Minervarya spp., Hoplobatrachus tigerinus and Euphlyctis spp. preferred the lowlands in proximal, central and southern parts of the study area proximal to urban settlements, with little vegetation and higher average temperatures. Duttaphrynus bengalensis and D. stomaticus had scattered distributions throughout the study area with no clear preference for elevation. Sphaerotheca pashchima was patchily distributed in the midwestern extent of the study area as well as the foothills to the north. Microhyla nilphamariensis was widely distributed throughout the study area with a preference for both lowlands and montane terrain. Endemic frogs (Nanorana vicina and Allopaa hazarensis) were observed only in locations with higher elevations, higher density of streams and lower average temperatures as compared to the other seven species sampled. It is recommended to provide legal protection to amphibians of Pakistan, especially endemic species, through revision in the existing wildlife laws. We suggest studying the effectiveness of existing amphibian tunnels and corridors or designing new ones tailored to the needs of our species to prevent their local extinction due to ongoing or proposed urban development which might affect their dispersal and colonization.
Citation: Rais M, Nawaz MA, Gray RJ, Qadir W, Ali SM, Saeed M, et al. (2023) Niche suitability and spatial distribution patterns of anurans in a unique Ecoregion mosaic of Northern Pakistan. PLoS ONE 18(6): e0285867. https://doi.org/10.1371/journal.pone.0285867
Editor: Karen Root, Bowling Green State University, UNITED STATES
Received: October 28, 2021; Accepted: May 3, 2023; Published: June 15, 2023
Copyright: © 2023 Rais et al. This is an open access article distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original author and source are credited.
Data Availability: The data underpinning the analysis, R codes and all other files used in the Maxent reported in this paper are deposited in the Dryad Data Repository at Dryad (https://doi.org/10.5061/dryad.cz8w9gj74).
Funding: MR received the funds from International Foundation for Science, Sweden (Grant No. D-5459/1) and Higher Education Commission, Pakistan (Project No. 3172) and designed the study design and prepared the manuscript. The funders had no role in study design, data collection and analysis, decision to publish, or preparation of the manuscript.
Competing interests: The authors have declared that no competing interests exist.
Introduction
The Wallacean shortfall, which suggests our understanding of geographical distribution patterns of the species at large scales is generally poor, can be minimized by sampling ecosystem gradients at smaller scales and expanding our knowledge outward [1]. Despite some progress in achieving global strategic goals of the United Nations and the Aichi Targets, biodiversity has declined, particularly in developing countries [2]. The lack of information regarding biodiversity status precludes recognition of impacts of anthropogenic activities on biodiversity and implementation of conservation strategies and targets [3]. Therefore, to gain a large-scale understanding of any taxonomic group and to determine threats and conservation potential within their spatial distribution in accordance with environmental variability, there must first be fine-scale (local/regional) baseline data.
Amphibian occurrence and abundance are greatly influenced by localized variation in geomorphic, geologic, and environmental characteristics [4]. Species Distribution Modelling (SDM) approaches based on Geographical Information System (GIS) have been widely used to predict distribution of species of amphibians and other vertebrate groups such as rodents and passerine birds [5–11] for conservation purposes. However, with presence-only records often being the only available data, and potential for imperfect detection within a study location, one must account for distribution and absence given contributing environmental covariates [12,13].
Amphibians in Pakistan have long been ignored in research, conservation, management, and policy and legislation. Pakistan’s National Climate Change Policy (2012), National Biodiversity Strategy and Action Plan (2015) and Pakistan Wetland Action Plan (2000) proposed guidelines for the conservation of natural resources, including fauna and flora, and mitigation of threats [14,15]. Currently, 21 species of amphibian (all anurans) have been documented in Pakistan [16] of which nine are believed to be endemic [17]. However, no progress on integrating anurans in wildlife conservation, and policy development or legislation has been made. There is no national assessment of conservation status of anurans of Pakistan which cautions the use of global conservation status. Only a few published studies report the richness and abundance of various Pakistani anuran species [18–24]. For example, the common Skittering Frog (Euphlyctis cyanophlyctis) [18] was reported as abundant in the rice fields of Gujranwala, Punjab Province [19], and nine anuran species were recorded from Margalla Hills National Park, Islamabad [20]. Six anuran species were reported from Rawalpindi and Islamabad areas [21,22] including high abundance of Indus Valley Bull Frog (Hoplobatrachus tigerinus) [23] and Skittering Frog from Rawal Lake, Islamabad [24].
In the current study, we aimed to estimate niche suitability, niche overlap and model distribution of anuran species within the Rawalpindi District and Islamabad Capital Territory, which encompasses a unique mosaic of deserts and xeric shrublands, montane grasslands, temperate broadleaf and mixed forests, temperate conifer forests, and tropical/subtropical coniferous forests [25]. Our findings were expected to generate new information on factors affecting niche suitability, spatial distribution size, and important ecoregions of the studied species and for anuran diversity to inform future conservation plans in Pakistan.
Materials and methods
Study area
The study was carried out in seven administrative units (Gujar Khan, Kahuta, Kallar Sayedan, Kotli Sattian, Murree, Rawalpindi, and Taxila) of the Rawalpindi District and Islamabad Capital Territory, Pakistan. The Rawalpindi District (33.4620° N, 73.3709° E) is located in the north-west of Punjab Province and covers an area of 5312 km2. Islamabad Capital Territory (ICT) (33.7205° N, 73.0405° E) is in the north-east of the country and covers an area of 906.50 km2. The elevation ranges from 457–2286 m and 457–610 m in Rawalpindi District and Islamabad Capital Territory, respectively. The climate of the study area is humid subtropical (Koppen climate classification). The summers produce more rain than the winter due to the monsoon season (July-August). The average rainfall in Islamabad Capital Territory is about 940 mm, in areas of Rawalpindi, Gujar Khan and Taxila ranges from 970–990 mm and in Murree, Kotli Sattian, and Kahuta is 1249 mm [26–28].
The study area is dominated by tropical and subtropical coniferous forests in the north; broad-leaf, mixed forest and montane grasslands and shrublands in the mid and midwest; and arid shrublands in the east and south (Fig 1). The tropical and subtropical coniferous forests are dominated by Pinus wallichiana and Pinus roxburgii and have relatively fewer human settlements. The proximal, central and southern regions feature urban and semi-urban areas with vegetation species such as Acacia modesta, along with Olea cuspidata, and Dodonea viscosa [29].
Map showing eleven major ecoregions in Pakistan (left) and ecoregions present in the extent of the present study area (right). The study area is shown with black boundary and grey shaded area while dotted lines represent river and watershed distributions throughout the areas. The study area is dominated by tropical and subtropical coniferous forests in the north; broad-leaf and mixed forests, montane grasslands, and montane shrublands in the mid and midwest; and arid shrublands in the east and south. The capital city Islamabad (shown as a blue box) is located in the midwestern portion of the study area. Ecoregion boundaries obtained from RESOLVE Ecoregions and Biomes database (Bioscience, An Ecoregions-Based Approach to Protecting Half the Terrestrial Realm DOI: https://doi.org/10.1093/biosci/bix014), available in ArcGIS Online under a CC by 4.0 license.
Of the 9 anuran species that are endemic to Pakistan [17], and therefore are of increased conservation interest, our study area encompasses most of the range of two of these endemics within Pakistan. Hazara Torrent Frog (Allopaa hazarensis) [30] is endemic to the springs and streams of Northern Pakistan [30]. Murree Hills Frog (Nanorana vicina) [31] is endemic to Pakistan and India [32].
Data collection
We surveyed 87 randomly selected sites in the seven administrative units of Rawalpindi District and Islamabad Capital Territory and gathered anuran presence data from 2016 to 2018 (March–September) on a monthly basis; we visited each site at least twice from during the study period. The duration and number of observers varied during the surveys; however, on average each visit consisted of two days of non-standard field observations from early morning till late night carried out on a weekly basis, with 4–6 observers. Each visit had at least one observer familiar with the identification of anurans of the area. Anuran species richness was recorded using time-constrained visual encounter survey (VES) technique [33]. Observers actively and thoroughly searched the sampling locations for a predefined time (~60–120 min.) in order to record the species. Adult specimens and tadpoles were collected with dip nets or simply picked up by hand, were examined, and identified following descriptive identification by Khan [29] and released back.
Data analysis and species distribution modelling
All analyses, models, maps and plots were generated using R statistical software version 3.6.3 [34]. A preliminary check on observed environmental values for each species was conducted by using the extract function in the ‘raster’ package [35]. Environmental ranges of species were checked before running models in order to identify potential flaws in the output (e.g. high spatial distribution probability in montane ranges when the species observations have initial low-elevation values). Bioclimatic variables (BIO 1–19) were retrieved from the WorldClim repository [36] via the ‘raster’ package at 0.5 arcmin resolution, and were tested for collinearity using Variation Inflated Factors (VIF). The following bioclimatic variables were retained: BIO7: temperature (C°) annual range; BIO8: mean temperature (C°) of wettest quarter (3 consecutive months); BIO9: mean temperature (C°) of driest quarter; BIO15: precipitation seasonality; BIO18: precipitation (mm) of the average warmest quarter; and BIO19: precipitation (mm) of the average coldest quarter. To map vegetation, we used MODIS (https://developers.google.com/earthengine/datasets/catalog/MODIS_006_MOD13Q1#bands; https://code.earthengine.google.com/e3a10b1ec6086c3ee7c598cfaca7dd98; resolution: 250m; scale 0.0001). MODIS has a spatial resolution of 250 m and provides a Vegetation Index (VI) value at a per pixel basis. We used two primary vegetation layers: Normalized Difference Vegetation Index (NDVI), which is the continuity index to the existing National Oceanic and Atmospheric Administration-Advanced Very High-Resolution Radiometer (NOAA-AVHRR) derived NDVI; and Enhanced Vegetation Index (EVI), which minimizes canopy background variations and maintains sensitivity over dense vegetation conditions. In Google Earth Engine (GEE), both MODIS layers The MODIS NDVI and EVI products were computed from atmospherically corrected bi-directional surface reflectance that had been masked for water, clouds, heavy aerosols, and cloud shadows. A 30 m Digital Elevation Model (DEM) was retrieved from the USGS earth explorer database (https://earthexplorer.usgs.gov/). The DEM was used to generate slope, aspect, and terrain roughness via the terrain function in the ‘raster’ package [37]. Distance to rivers was calculated using a local shape file of all waterways in Pakistan and the fasterVectToRastDistance function of the ‘fasterRaster’ package [38]. River distance values (in meters) were log normalized to prevent outlier values from overwhelming the model contribution estimates. All environmental variables were resampled to 0.5 arcmin resolutions thereafter using the bilinear method, and masked to a shapefile of Pakistan regional boundaries. The explanatory variables were then stacked and tested for multicollinearity using the vifstep and vifcor functions of the ‘usdm’ package [39] with a threshold VIF of 10 and a threshold correlation coefficient with an absolute value of 0.7; whereby we excluded predictors that were strongly correlated with ⸺ but considered to be less important for anurans ⸺ than other predictors.
The Maximum Entropy or “Maxent” modelling method has been shown to perform as well as, or better than ensemble models when modelling species distributions, with additional benefits of lower computational power requirements and increased simplicity of use [40]. Species Distribution Models for this study were created following guidelines on Maxent parameterization [41,42]. Species occurrence data was thinned to one point per raster cell to omit spatial biases. Maxent models were created using the Maxent software wrapper through the ‘dismo’ package [37]. A kernel density bias mask was created by querying all available anuran occurrences from the Global Biodiversity Information Facility (GBIF) database for Pakistan (S1 File) as a measure of restrictive effort using the ‘MASS’ package [43] using the reference bandwidth smoothing factor. We used the occ_search function in the ‘rgbif’ package [44] in R to access the GBIF database, using the taxon “Anura”, the country Pakistan (“PK”), and “coordinates = TRUE” as filters to create a text file containing 345 dataset keys and number of observations per dataset) (S2 File). We generated 10,000 geographically randomized background points within the bias mask estimate. Models were tuned using the ENMeval function in the ‘ENMeval’ package [45] to identify feature selection variables and regularization multiplier (beta-multiplier) values selected from the lowest delta Akaike Information Criterion (AICc) using five random k-folds. Maxent datasets were partitioned into training (3/4) and testing (1/4) data using a 4-way partitioned k-fold. Features selected in the final models were linear and quadratic (lq) with a beta-multiplier set to between 1–2 for sensitivity testing, and 1 for the final model. We selected these features for the final Maxent models for each species based on the statistics we used to validate our model: area-under-the-curve (AUC) values within receiver-operating characteristic (ROC) curves and the true skill statistic (TSS [= TP + TN—1, where TP = proportion of true positive predictions and TN = proportion of true negative predictions]). Models were replicated with replacement using the bootstrap method.
Final models were visualized using the raw model output using the ‘ggplot2’ package [46] and each model prediction was condensed using the highest true skill statistics (TSS) threshold value as a filter (hereafter called “threshold model”). The true skill statistics is defined based on the components of the standard confusion matrix representing matches and mismatches between observations and predictions [47]. Model outputs were defined as niche suitability on a 0–1 probability scale (1 = highest niche suitability; 0 = lowest niche suitability), and threshold models were considered to be the niche distribution of a given species in terms of overlap with others. Each threshold model was converted to points and extracted from the explanatory variables to derive summary statistics of environmental values for each species within their niche distribution. All species threshold predictions were then combined to create niche overlap raster to determine which eco-regions were most suitable for anuran species in Pakistan. High-density niche overlap areas were identified by filtering 50% (4.5 species per raster cell) and 100% (9 species per raster cell). We then converted the two final niche spaces to a spatial polygon to identify percentage overlap between each ecoregion and high-density niche distributions.
Results
Niche suitability and overlap
Final predictions for niche suitability of all nine anuran species using Maxent models indicated good fit in terms of AUC (Fig 2). Environmental variables varied among species (Fig 2); however, distance to rivers and precipitation in the warmest and coolest quarters were shown to have significant contribution (% contribution >10; higher permutation importance and variable importance in jacknife-testing) (Fig 2) for Duttaphrynus bengalensis and M. nilphamariensis; distance to rivers and precipitation of the warmest quarter for Duttaphrynus stomaticus; distance to rivers, precipitation of the warmest quarter, and NDVI for A. hazarensis and N. vicina; and precipitation of the warmest quarter and NDVI were most influential for H. tigerinus, Euphlyctis spp., Minervarya spp., and Sphaerotheca pashchima (S2–S10 Files). Habitat suitability for most species was higher in locations with greater precipitation in the warmest quarter and less precipitation in the coolest quarter. Habitat suitability for A. hazarensis and N. vicina was higher in locations with greater ecosystem productivity (higher NDVI) closer to rivers (S3–S11 Files).
Numbers above each box plot indicate measures of average model run performance (AUC = area-under-the-curve; TSS = true skill statistic). Numbers below each box plot indicate the percent contribution (top row) and permutation importance (bottom row) of each variable in predicting habitat suitability for each species, averaged over ten model runs per species. BIO7: Temperature (C°) annual range; BIO8: Mean temperature (C°) of wettest quarter; BIO9: Mean temperature (C°) of driest quarter; BIO15: Precipitation seasonality (%); BIO18: Precipitation (mm) of the warmest quarter, BIO19 Precipitation (mm) of the coldest quarter; NDVI: Normalized Difference Vegetation Index; RivDist: (log) Distance to Rivers (m).
Spatial distribution
Nine anuran species showed a divergent pattern of spatial distribution (Fig 3). We found that Euphlyctis spp., Minervarya spp., and H. tigerinus showed preference for the lowlands in proximal, central and southern parts of the study with high urban settlements, little vegetation and higher average temperatures. The toads (D. bengalensis and D. stomaticus) had scattered distributions throughout the study area with no clear preference for elevation. S. pashchima showed a patchy distribution in the midwestern extent of the study area as well as the foothills to the north. Microhyla nilphamariensis also showed a wide distribution throughout the study area with a preference for both lowlands and montane terrain. Endemic species such as N. vicina and A. hazarensis were observed in locations with higher elevations, proximal to streams and lower average temperatures as compared to the other seven species sampled (Fig 3).
The predictions were unconstrained which allowed the model to predict outward throughout Pakistan, but predictions mostly stayed in that habitat mosaic in the north.
Niche overlap at 50% (4.5 species per area) within the studied ecoregions showed 0.2% overlap in montane grasslands and shrublands, 4.8% in temperate broadleaf and mixed forests, 23.6% in temperate coniferous forests and 33.9% in tropical and subtropical coniferous forests. Whereas, the niche overlap (9 species per area) revealed 0% overlap in montane grasslands and shrublands, 3.2% overlap in deserts and xeric shrublands, 9.4% overlap in temperate broadleaf and mixed forests, 17% overlap in temperate coniferous forests, and 70.2% overlap in tropical and subtropical coniferous forests (Fig 4), providing evidence that tropical and subtropical coniferous forests support the highest diversity of anuran species in the region.
Lighter areas have more species (richness) and darker areas have less based on suitability predictions above the TSS threshold values for each species. A species was predicted to be present (= 1) if a location’s niche suitability for that species exceeded the TSS threshold for that species; otherwise, the species was predicted to be absent (= 0). Species richness was then summed for each location. Of locations with maximal amphibian species overlap (9 species), 70.2% occurred within tropical/subtropical coniferous forests (pink polygons); 17.0% occurred within temperate coniferous forests (dark green); 9.4% occurred within temperate broadleaf and mixedwood forests (medium green); 3.2% occurred within xeric desert and shrublands (pale blue); and 0% occurred in montane grasslands and shrublands (lightest green). Ecoregion boundaries obtained from RESOLVE Ecoregions and Biomes database (Bioscience, An Ecoregions-Based Approach to Protecting Half the Terrestrial Realm DOI: https://doi.org/10.1093/biosci/bix014), available in ArcGIS Online under a CC by 4.0 license.
Discussion
Development of analytical habitat distribution models has rapidly increased in ecology due to the invention of new GIS tools and statistical techniques. Such models statistically relate the geographical distribution of species or communities to their present environment. Neither any scientific study on the distribution patterns nor any species occurrence database of anurans of Pakistan exists in the country. In forested mountainous regions outside of but similar to those of Pakistan, amphibians, particularly Ascaphus truei and Dicamptodon tenebrosus, were more abundant in stream habitats within older coniferous forests [48]. Precipitation and soil temperature influences probability of occurrence for Polish amphibian species [10] while vegetation contributed significantly in the prediction for salamander species of central Portugal [11]. Our results suggest similar trends, with precipitation (of the warmest and coldest quarter), distance to rivers and vegetation being highly deterministic factors of suitability for anurans in Pakistan.
While Maxent models can accurately predict grid-based habitat suitability and presence of species when observation data of those species is limited, predictions from such models could also be used to identify the likeliest locations of species for further monitoring and obtain actual presence/absence or abundance data. Actual presence/absence or abundance data at precise locations enables biologists to use finer-scaled environmental variables that influence habitat suitability and abundance of species. For anurans and other aquatic wildlife, the type of wetland habitat or changes in wetland variables over time at repeatedly monitored sites can increasingly be quantified or classified at fine spatial scales over large regions away from ground-truthed locations [49–54]. Increased availability of wetland data is due to: 1) the accumulation of decades of remotely sensed environmental data by satellites; 2) the development of free, open-source, online platforms (e.g., Google Earth Engine) for processing remotely sensed data and extracting this data to survey locations; and 3) development and sharing of open-source machine learning techniques for predicting and classifying wetlands [49–54]. To improve the chances of detecting anurans or other vocalising species when they are present at sites, monitoring could involve obtaining multiple visits per site by using passive acoustic monitoring with pre-programmed acoustic recorders [55,56]. Acoustic data could then be transcribed to obtain detections per visit and hierarchical models can be used to account for varying detection probability of species among sites when estimating effects of environmental variables on occupancy or abundance [57].
The toads (Family Bufonidae) of the study area: D. bengalensis and D. stomaticus, had somewhat scattered distribution throughout the region with no clear preference for elevations; however, D. stomaticus showed a clear preference to the northwestern extent of the study area, mostly concentrating within and proximal to the lowlands. D. bengalensis and D. stomaticus have been recorded as widespread species found up to 1800 m [58] and 4500 m [59] elevations, respectively. D. bengalensis is adapted to various types of habitats, even degraded ones, and around human habitations [60–62]. The geographic range of D. bengalensis has now been extended due to its introduction and D. bengalensis has attained a status of invasive species in various parts of the world [58]. The two toads are found in the plains, lowlands, sub mountain areas as well as hilly areas in Pakistan which experience monsoon season (July-August) during which the species breed [29]. We recorded D. bengalensis at elevations higher than previously reported. One reason for their widespread distribution in our study area is their adaptation to a wide range of habitats. M. nilphamariensis, a diminutive frog species (Family Microhylidae), showed a widespread distribution throughout the study area with a preference for both lowlands and montane terrain. The species has been recorded from areas up to 2000 m elevation [63] as well as from lowlands, sub mountain areas and foothills in Pakistan [29].
The dicroglossid frogs of the study area showed varied distribution patterns. S. pashchima, a burrowing frog species, showed a patchy distribution in the midwestern extent of the study area as well as the foothills to the north. S. pashchima has been recorded as a widespread species from lowlands and forested areas up to 1500 m [64]. This species remains under soft soil for most parts of the year and emerges during summers to breed during the monsoon [29] and avoids high altitude areas possibly due to their burrowing habit, since hard substrate makes it difficult for them to dig. Further, mountains in the north are under less influence of monsoon, which likely provides less suitable breeding conditions. Other dicroglossid frogs such as Euphlyctis spp. showed scattered distribution probability through the lowlands. Minervarya spp. was shown to prefer proximal, central, southern, and western lowland areas. H. tigerinus also showed variability in preference between lowland and elevated areas, with a concentrated distribution toward the middle to northern extents of the study area. H. tigerinus has previously been recorded from areas up to 2000 m [65]. These areas are amongst most built up parts of the region in addition to encompassing other human modified habitats such as croplands [44,66]. Euphlyctis spp. has been recorded from areas up to 2500 m [67] while Fejervarya spp. from areas up to 2000 m (Dijk 2004) [58], but these dicroglossid frogs are also widespread in lowlands and forested areas. Euphlyctis cyanophlyctis and Zakerana syhadrensis occur along stream banks and water pools between forest edges, agricultural areas, and residential gardens [62]. Our findings are consistent with the available information. However, we have thus forth provided empirical data on the response of these species to the studied environmental factors for the first time.
The endemic frogs, A. hazarensis and N. vicina, showed restricted occurrence within the northern and north-eastern mountain ranges in areas of high elevation compared to the proximal lowlands to the southeast. Most of these areas feature subtropical pine forest (900–1500 m) dominated by Pinus roxburghii trees, while the northernmost areas possess Himalayan moist temperate forest (1500–3000 m) dominated by Pinus wallichiana and Pinus roxburghii. The wetlands throughout this range exist in the form of freshwater streams [29]. A. hazarensis is endemic to Pakistan while N. vicina is known from Pakistan and India. A. hazarensis and N. vicina are known from streams and pools in forested mountainous areas as high as 1500 m [68] and 3000 m, respectively [64]. During our study, we found that A. hazarensis could occur at a higher elevation (>1500) than the previously reported range.
Conservation implications and suggestions
Urban developments have been shown to have negative impacts on amphibians [69,70]. Several types of culverts, tunnels and corridors have been developed and their effectiveness assessed. Many amphibian species in North America and Europe have used these structures, resulting in reduced mortality from vehicular collisions [71]. Since anurans in Pakistan enjoy no legal protection, no such consideration is given during urban planning. With recent urban expansion in Rawalpindi District and the Islamabad Capital Territory [44,66], there has been a noticeable reduction in forests, open spaces and watersheds [72]. As shown by our results, some of the anuran study species are tolerant to habitat degradation while others are not. The creation of human modified habitats may further facilitate the spread of native species such as Euphlyctis spp., while also accommodating invasive species such as D. bengalensis and D. stomaticus which may pose the threat of resource competition against native species. We recommend studying the effectiveness of such existing amphibian tunnels and corridors or designing new ones tailored to the needs of our species. These could then be incorporated in future urban development programs.
Conservationists put more emphasis on conserving threatened species. Studies, however, have shown that even common species are subjected to population decline and local extinction especially if these species are associated with a particular type of habitat or set of environmental conditions [73,74]. A 29% decrease in the population of the Moor Frog (Rana arvalis), which is found in heathlands and moorlands in parts of Europe, was reported between 1950–2006 [73]. This decline was attributed to cultivation of heathlands and moorlands, lowering of ground water levels and intensification of agricultural practices. In a forest landscape study, Brown Frogs (Rana arvalis, R. temporaria) bred more in various wetland habitats (e. g., naturally flooded areas, beaver ponds, mitigation pools with shallow littoral zones, cleaned ditches) than in ditches overgrown with forest vegetation [74].
Climate change during the past two decades has affected several species of plants and animals in Pakistan [75]. The tadpoles of A. hazarensis and N. vicina responded (under laboratory conditions) to higher temperature (>26°C) through faster metamorphosis, reduction in the body size, more frequent developmental complications or deformities such as edema and tail kinks, lower fitness and higher mortality [76]. Being associated with a particular set of environmental conditions in the north and northeast of the study area, it is feared that these two species endemic frogs, which are currently evaluated as least concern in the IUCN Red List of threatened species [77], may experience local extinction in the future. Land use simulators can be used to project changes over space and time in environmental variables used both in Maxent models (e.g., climate) and in models based on surveys whose locations were informed by Maxent model predictions. Thus, Maxent models can be used to project changes in distribution of species over time under different climate scenarios, and raster layers based on these distributions can be used to identify potential refugia for anurans either under current or future climate conditions [78]. These raster layers may also be used as inputs in raster overlay-based conservation planning tools (e.g., Marxan, Zonation) to prioritise locations for protection or management of threatened species [79].
Supporting information
S3 File. Maxent Output for Allopaa hazarensis.
https://doi.org/10.1371/journal.pone.0285867.s003
(PDF)
S4 File. Maxent Output for Duttaphrynus bengalensis.
https://doi.org/10.1371/journal.pone.0285867.s004
(PDF)
S5 File. Maxent Output for Duttaphrynus stomaticus.
https://doi.org/10.1371/journal.pone.0285867.s005
(PDF)
S7 File. Maxent Output for Hoplobatrachus tigerinus.
https://doi.org/10.1371/journal.pone.0285867.s007
(PDF)
S8 File. Maxent Output for Microhyla nilphamariensis.
https://doi.org/10.1371/journal.pone.0285867.s008
(PDF)
S11 File. Maxent Output for Sphaerotheca pashchima.
https://doi.org/10.1371/journal.pone.0285867.s011
(PDF)
Acknowledgments
We would like to especially thank Mr. Wajehuddin and Mr. A. Bhutta, Punjab Forest Department, Government of Punjab, Pakistan for providing assistance during the field work. We are thankful to Snow Leopard Foundation, Pakistan, for providing logistics. We are also thankful to Dr. Hussain Ali, Snow Leopard Foundation, for his help.
References
- 1. Whittaker RJ, Araújo MB, Jepson P, Ladle RJ, Watson JE, Willis KJ. Conservation biogeography: assessment and prospect. Diversity and distributions. 2005; 11(1):3–23.
- 2. Zhang Y, Tariq A, Hughes AC, Hong D, Wei F, Sun H, et al. 2023. Challenges and solutions to biodiversity conservation in arid lands. Science of the Total Environment. 857; (3). 159695. pmid:36302433
- 3. Mihoub JB, Henle K, Titeux N, Brotons L, Brummitt NA, Schmeller DS. Setting temporal baselines for biodiversity: the limits of available monitoring data for capturing the full impact of anthropogenic pressures. Scientific Report. 2017; 7(1):1–13. https://doi.org/10.1038/srep41591.
- 4. Wilkins RN, Peterson NP. Factors related to amphibian occurrence and abundance in headwater streams draining second-growth Douglas-fir forests in southwestern Washington. Forest Ecology and Management. 2000; 139:79–91. https://doi.org/10.1016/S0378-1127(99)00336-9.
- 5. Weiers S, Bock M, Wissen M, Rossner G. Mapping and indicator approaches for the assessment of habitats at different scales using remote sensing and GIS methods. Landscape Urban Planning. 2004; 67:43–65. https://doi.org/10.1016/S0169-2046(03)00028-8.
- 6. Anderson RP, Lew D, Peterson AT. Evaluating predictive models of species’ distributions: Criteria for selecting optimal models. Ecological Modelling. 2003; 162: 211–232. https://doi.org/10.1016/S0304-3800(02)00349-6.
- 7. Arntzen JW. From descriptive to predictive distribution models: A working example with Iberian amphibians and reptiles. Frontiers in Zoology. 2006; 3:8. pmid:16674803
- 8. Dayton GH, Fitzgerald LA. Habitat suitability models for desert amphibians. Biological Conservation. 2006; 132: 40–49. https://doi.org/10.1016/j.biocon.2006.03.012.
- 9. Kaky E. Potential habitat suitability of Iraqi amphibians under climate change. Biodiversitas Journal of Biological Diversity. 2020; 21: 731–742. https://doi.org/10.13057/biodiv/d210240.
- 10.
Lai J. Amphibian Species Distribution Modelling in Poland.2009; (Master’s thesis, University of Twente).
- 11.
Negga H. Predictive modelling of amphibian distribution using ecological survey data: a case study of central Portugal. 2007; Mphil thesis. ITC, Netherlands.
- 12. Rodríguez-Rodríguez EJ, Beltrán JF, Tejedo M, Nicieza AG, Llusia D, Márquez R, et al. Niche models at inter-and intraspecific levels reveal hierarchical niche differentiation in midwife toads. Scientific Reports. 2020;10(1):1–1. https://doi.org/10.1038/s41598-020-67992-6.
- 13. Lahoz-Monfort JJ, Guillera-Arroita G, Wintle BA. Imperfect detection impacts the performance of species distribution models. Global Ecology and Biogeography. 2014; 23: 504–515. https://doi. 0.1111/geb.12138.
- 14. GOP. Government of Pakistan. The Pakistan National Conservation Strategy. Environment and Urban Affairs Division. 1992; Government of Pakistan and IUCN-Pakistan.
- 15. Khurshid SN. Pakistan’s Wetlands Action Plan. World Wild Fund and National Council for Conservation of Wildlife. 2000; Government of Pakistan.
- 16. Rais M, Ahmed W, Sajjad A, Akram A, Saeed M, Hamid HN, et al. Amphibian fauna of Pakistan with notes on future prospects of research and conservation. ZooKeys. 2021; 1062:157. pmid:34720620
- 17. Ali W, Javid A, Hussain A, Bukhari SM. Diversity and habitat preferences of amphibians and reptiles in Pakistan: a review. Journal of Asia-Pacific Biodiversity. 2018; 11(2):173–87. https://doi.org/10.1016/j.japb.2018.01.009.
- 18. Schneider JG. Historia Amphibiorum Naturalis et Literarariae. Fasciculus Primus. Continens Ranas, Calamitas, Bufones, Salamandras et Hydros in Genera et Species Descriptos Notisque suis Distinctos. 1799; Jena: Friederici Frommanni.
- 19. Yousaf S, Mahmood T, Rais M, Qureshi IZ. Population variation and food habits of ranid frogs in the rice-based cropping system in Gujranwala region, Pakistan. Asian Herpetological Research. 2010; 1(2):123–30.
- 20. Masroor R. An annotated checklist of amphibians and reptiles of Margalla Hills National Park, Pakistan. Pakistan Journal of Zoology. 2011; 43:1041–1048.
- 21. Rais M, Baloch S, Rehman J, Anwar M, Hussain I, Mahmood T. Diversity and conservation of amphibians and reptiles in North Punjab, Pakistan. Herpetology Bulletin 2012; 16–25.
- 22. Akram A, Rais M, Asadi MA, Jilani MJ, Balouch S, Anwar M, et al. Do habitat variables correlate anuran abundance in arid terrain of Rawalpindi–Islamabad Areas, Pakistan?. Journal of King Saud University-Science. 2015; 27(3):278–83. https://doi.org/10.1016/j.jksus.2015.02.001.
- 23. Daudin FM. "An. XI". Histoire Naturelle des Rainettes, des Grenouilleset des Crapauds. 1802; Quarto version. Paris: Levrault.
- 24. Tabassum F, Rais M, Anwar M, Mehmood T, Hussain I, Khan SA. Abundance and breeding of the common skittering frog (Euphlyctis cyanophlyctis) and bull frog (Hoplobatrachus tigerinus) at Rawal Lake, Islamabad, Pakistan. Asian Herpetological Resource. 2011; 2:245–50.
- 25. Dinerstein E, Olson D, Joshi A, Vynne C, Burgess ND, Wikramanayake E, et al. An ecoregion-based approach to protecting half the terrestrial realm. BioScience. 2017; 67(6):534–45. pmid:28608869
- 26. GOP. Government of Pakistan. Summary environmental impact assessment Rawalpindi environmental improvement project in the Islamic republic of Pakistan. 2005; Government of Pakistan, Islamabad.
- 27. GOP. Government of Pakistan [Internet]. Meteorological data of Rawalpindi from 1931 to 2006. 2006 [cited 2019 June 05]. Available from: http://www.met.gov.pk/cdpc/islamabad.html. Accessed:
- 28. Malik RN, Husain SZ. Classification and ordination of vegetation communities of the Lohibehr reserve forest and its surrounding areas, Rawalpindi, Pakistan. Pakistan Journal of Botany. 2006; 38(3):543.
- 29. Khan MS. Amphibians and reptiles of Pakistan. Malabar, Florida, USA: Krieger Publishing Company. 2006; 1–311
- 30. Dubois A, Khan MS. A new species of frog (genus Rana, subgenus Paa) from northern Pakistan (Amphibia, Anura). Journal of Herpetology. 1979; 403–10.
- 31. Stoliczka F. Notes on some new species of Reptilia and Amphibia, collected by Dr. W. Waagen in North-Western Punjab. In Proceedings of the Asiatic Society of Bengal 1872 (Vol. 1872, pp. 124–131).
- 32. Rais M, Abbassi S, Batool T, Jilani MJ, Assadi MA, Mubarak H, et al. A note on recapture of Nanorana vicina (Anura: Amphibia) from Murree, Pakistan. Journal of Animal and Plant Sciences. 2014; 24:455–8.
- 33. Graeter GJ, Buhlmann KA, Wilkinson LR, Gibbons JW. Inventory and monitoring: Recommended techniques for reptiles and amphibians. Partners in Amphibian and Reptile Conservation Technical Publication IM-1, Birmingham, Alabama. 2013.
- 34.
R Core Team. [Internet]. R A language and environment for statistical computing. R Foundation for Statistical Computing, Vienna, Austria. Version 3.6.3. 2020 [cited 2019 June 05] Available from: https://www.R-project.org/.
- 35. Hijmans RJ, Van Etten J, Sumner M, Cheng J, Baston D, Bevan A, et al. [Internet]. Package ‘raster’. R package. Version 3.3–13. 2015 [cited 2019 June 05]. Available from https://cran.r project.org/web/packages/raster/index.html.
- 36. Hijmans RJ, Cameron SE, Parra JL, Jones PG, Jarvis A. Very high-resolution interpolated climate surfaces for global land areas. International Journal of Climatology: Journal of the Royal Meteorological Society. 2005; 25: 1965–1978. https://doi.org/10.1002/joc.1276.
- 37. Hijmans RJ, Phillips S, Leathwick J, Elith J. [Internet]. Package ‘dismo’. 2011[cited 2019 June 05]. Available from: http://cran.r-project.org/web/packages/dismo/index.html.
- 38. Smith BA. [Internet]. Faster Raster: Faster raster processing in R using GRASS GIS. R package version 0.6.0. 2020 [cited 2019 June 05]. Available from: http://www.earthSkySea.org.
- 39. Naimi B. Uncertainty analysis for species distribution models. R package version. 2015; 1:1–12.
- 40. Kaky E, Nolan V, Alatawi A, Gilbert F. A comparison between Ensemble and MaxEnt species distribution modelling approaches for conservation: A case study with Egyptian medicinal plants. Ecological Informatics. 2020; 60:101150. https://doi.org/10.1016/j.ecoinf.2020.101150.
- 41. Feng X, Gebresenbet F, Walker C. Shifting from closed-source graphical-interface to open-source programming environment: a brief tutorial on running Maxent in R (No. e3346v1). Peer J Preprints. 2017.
- 42. VanDerWal J, Shoo LP, Graham C, Williams SE. Selecting pseudo-absence data for presence-only distribution modeling: how far should you stray from what you know?. Ecological modelling. 2009; 220(4):589–94. https://doi.org/10.1016/j.ecolmodel.2008.11.010.
- 43.
Venables WN, Ripley BD. [Internet]. Modern Applied Statistics with S, Fourth edition. Springer, New York. ISBN 0-387-95457-0, 2002 [cited 2019 June 05]. Available from: https://www.stats.ox.ac.uk/pub/MASS4/.
- 44. Zahra H, Rabia S, Sheikh SA, Amir HM, Neelam A, Amna B, Summra Eet al. Dynamics of land use and land cover change (LULCC) using geospatial techniques: A case study of Islamabad, Pakistan. Springer Plus. 2016; 5:8–12. pmid:27390652
- 45. Muscarella R, Galante PJ, Soley‐Guardia M, Boria RA, Kass JM, Uriarte M, Anderson RPet al. ENM eval: An R package for conducting spatially independent evaluations and estimating optimal model complexity for Maxent ecological niche models. Methods in Ecology and Evolution. 2014; 5(11):1198–205. https://doi.org/10.1111/2041-210X.12261.
- 46.
Wickham H. [Internet]. 2016.ggplot2: Elegant Graphics for Data Analysis. Springer-Verlag New York. ISBN 978-3-319-24277-4.2016 [cited 2019 June 05]. Available from: https://ggplot2.tidyverse.org.
- 47. Fielding AH, Bell JF. A review of methods for the assessment of prediction errors in conservation presence/absence models. Environmental conservation. 1997; (1):38–49.
- 48. Welsh H, Lind A. Habitat Correlates of the Southern Torrent Salamander Rhyacotriton variegatus (Caudata: Rhyacotritonidae) in Northwestern California. Journal of Herpetology. 1996;30:385–398. https://doi.org/10.2307/1565176.
- 49. Alonso A, Muñoz-Carpena R, Kennedy RE, Murcia C. Wetland landscape spatio-temporal degradation dynamics using the new Google Earth Engine cloud-based platform: Opportunities for non-specialists in remote sensing. Transactions of the ASABE. 2016; 59(5):1331–42.
- 50. Dong J, Xiao X, Menarguez MA, Zhang G, Qin Y, Thau D, Biradar C, Moore III Bet al. Mapping paddy rice planting area in northeastern Asia with Landsat 8 images, phenology-based algorithm and Google Earth Engine. Remote Sensing of Environment. 2016; 185:142–54 https://doi.org/10.1016/j.rse.2016.02.016. pmid:28025586
- 51. Chen B, Xiao X, Li X, Pan L, Doughty R, Ma J, Dong J, Qin Y, Zhao B, Wu Z, Sun Ret al. A mangrove forest map of China in 2015: Analysis of time series Landsat 7/8 and Sentinel-1A imagery in Google Earth Engine cloud computing platform. ISPRS Journal of Photogrammetry and Remote Sensing. 2017; 131:104–20. https://doi.org/10.1016/j.isprsjprs.2017.07.011.
- 52. Hird JN, DeLancey ER, McDermid GJ, Kariyeva J. Google Earth Engine, open-access satellite data, and machine learning in support of large-area probabilistic wetland mapping. Remote Sensing. 2017; 9(12):1315.https://doi.org/10.3390/rs9121315.
- 53. Fekri E, Latifi H, Amani M, Zobeidinezhad A. A Training Sample Migration Method for Wetland Mapping and Monitoring Using Sentinel Data in Google Earth Engine. Remote Sensing. 2021; 13(20):4169.https://doi.org/10.3390/rs13204169.
- 54. Long X, Li X, Lin H, Zhang M. Mapping the vegetation distribution and dynamics of a wetland using adaptive-stacking and Google Earth Engine based on multi-source remote sensing data. International Journal of Applied Earth Observation and Geoinformation. 2021;102:102453. https://doi.org/10.1016/j.jag.2021.102453.
- 55. Ribeiro JW, Sugai LS, Campos-Cerqueira M. Passive acoustic monitoring as a complementary strategy to assess biodiversity in the Brazilian Amazonia. Biodiversity and Conservation. 2017; 26(12):2999–3002. https://doi.org/10.1007/s10531-017-1390-0.
- 56. Desjonquères C, Gifford T, Linke S. Passive acoustic monitoring as a potential tool to survey animal and ecosystem processes in freshwater environments. Freshwater Biology. 2020; 65(1):7–19. https://doi.org/10.1111/fwb.13356.
- 57. Mazerolle MJ, Desrochers A, Rochefort L. Landscape characteristics influence pond occupancy by frogs after accounting for detectability. Ecological Applications. 2005; 15(3):824–34. https://doi.org/10.1890/04-0502.
- 58. van Dijk PP, Iskandar D, Inger R, Lau MWN, Ermi Z, Baorong G, Dutta S, De-Silva KMA, Bordoloi S, Kaneko Y, Matsui M, Khan MSet al. [Internet]. Fejervarya limnocharis. The IUCN Red list of Threatened species 2004:e.T58275A86154107. 2004 [cited 2020 May 03]. Available from: https://dx.doi.org/10.2305/IUCN.UK.2004.RLTS.T58275A11747569.en.
- 59. Stöck M, Khan MS, Papenfuss T, Anderson S, Kuzmin S, Rastegar-Pouyani N, Dutta S, Ohler A, Sengupta S, Manamendra-Arachchi K, de Silva A, Anderson S, Sharifi Met al. [Internet]. Duttaphrynus stomaticus. The IUCN Red list of Threatened species 2009: e.T54768A11201081. 2009. [cited 2020 May 03]. Available from: en.
- 60.
Saidapur SK, Girish S. Growth and metamorphosis of Bufo melanostictus tadpoles: effects of kinship and density. Journal of Herpetology. 2001; 249–54. https://doi.org/10.2307/1566115.
- 61. Narzary J, Bordoloi S. A Study on certain common pond breeding Anurans and their tadpoles in a pond of Western Assam, India. International Journal of Advanced Biological Research. 2012; 2(2):342–8.
- 62. Botejue WMS, Wattavidanage J. Herpetofaunal diversity and distribution in Kalugala proposed forest reserve, Western province of Sri Lanka. Amphibian and Reptile Conservation 2012; 5: 65–80.
- 63. Dutta S, Shrestha TK, Arachchi KM, Khan MS, Roy D. 2008 [Internet]. Microhyla ornata. The IUCN Red List of Threatened Species 2008: e.T57886A11686884. 2008 [cited 2020 May 03] Availablefrom:https://dx.doi.org/10.2305/IUCN.UK.2008.RLTS.T57886A11686884.en.
- 64. Ohler A, Khan MS, Van Dijk PP, Wogan G, Dutta S, Inger RF, Kumar Shrestha T, M Manamendra-Arachchi K, de Silva Aet al. [Internet]. Sphaerotheca breviceps. The IUCN red list of threatened species 2004: e.T58755A86145885. 2004 [cited 2020 May 03] Availablefrom:https://dx.doi.org/10.2305/IUCN.UK.2004.RLTS.T58755A11837725.en.
- 65. Padhye A, Manamendra-Arachchi K, deSilva A, Dutta S, Kumar Shrestha T, Bordoloi S, Papenfuss T, Anderson S, Kuzmin S, Khan MS, Nussbaum Ret al. [Internet]. Hoplobatrachus tigerinus. The IUCN red list of threatened species 2008. 2008 [cited 2020 May 03] Available from: e.T58301A11760496. https://dx.doi.org/10.2305/IUCN.UK.2008.RLTS.T58301A11760496.en.
- 66. Mahboob MA, Atif I, Riaz A. Spatio temporal mapping and monitoring of land cover dynamics of Islamabad using multi-temporal remote sensing data. Pakistan Journal of Science. 2016; 68:146.
- 67. Khan MS, Papenfuss T, Anderson S, Rastegar-Pouyani N, Kuzmin S, Dutta S, Manamendra-Arachchi K, Sharifi Met al. [Internet]. Euphlyctis cyanophlyctis. The IUCN red list of threatened species 2009: e.T58260A86626211. 2009 [cited 2020 May 03]. Available from https://dx.doi.org/10.2305/IUCN.UK.2009.RLTS.T58260A11745753.en.
- 68. Khan MS, Dutta S, Ohler A. [Internet]. Allopaa hazarensis. The IUCN Red List of Threatened Species 2008:e.T58426A11779666. 2008 [cited 2020 May 03] Available from:https://dx.doi.org/10.2305/IUCN.UK.2008.RLTS.T58426A11779666.en.
- 69.
Andrews K,Gibbons JW,Jochimsen DM. Ecological effects of roads on amphibians and reptiles: a literature review. In Mitchell J.C., Jung Brown R.E., Bartholomew B. (Eds), Urban Herpetology, Society for the Study of Amphibians and Reptiles. Salt Lake City, Utah, USA. 2008; 121–143.
- 70. Becker , Fonseca CR, Haddad CF, Prado PI. Habitat split as a cause of local population declines of amphibians with aquatic larvae. Conservation Biology. 2010; 24:287–294. pmid:19758391
- 71.
Smith RK, Meredith H, Sutherland WJ. Amphibian Conservation. In Sutherland W.J., Dicks L.V., Ockendon N., Petrovan S.O., Smith R.K. (Eds), What works in conservation. 2019; Open Book Publishers, Cambridge, UK. 2019; 9–65.
- 72. Saeed MA, Ashraf A, Ahmed B, Shahid M. Monitoring deforestation and urbanization growth in Rawal watershed area using remote sensing and GIS techniques. A scientific journal of COMSATS. 2011; 17:19–104.
- 73. Delft JV, Creemers R. Distribution, status, and conservation of the in the Netherlands. Zeitschrift für Feldherpetologie. Supplement. 2008; 13:255–268.
- 74. Remm L, Vaikre M, Rannap R, Kohv M. Amphibians in drained forest landscapes: Conservation opportunities for commercial forests and protected sites. Forest Ecology and Management. 2018; 428:87–92. https://doi.org/10.1016/j.foreco.2018.06.038.
- 75. Saadat HB, Nawaz CM, Manzoor F, Nasim G. Effect of climate change on butterfly population of selected coniferous forests of Murree Hills and adjacent areas, Pakistan. Pakistan Journal of Zoology. 2016; 48:1963–1969.
- 76. Saeed M, Rais M, Gray RJ, Ahmed W, Akram A, Gill S, et al. Rise in temperature causes decreased fitness and higher extinction risks in endemic frogs at high altitude forested wetlands in northern Pakistan. Journal of Thermal Biology. 2021; 95:102809. pmid:33454039
- 77. Ohler A, Dutta S. 2004. Nanorana vicina: IUCN Red List of Threatened Species. International Union for Conservation of Nature. IUCN Red List of Threatened Species.
- 78. Barrett K, Nibbelink NP, Maerz JC. Identifying priority species and conservation opportunities under future climate scenarios: amphibians in a biodiversity hotspot. Journal of Fish and Wildlife Management. 2014; 5(2): 282–297. https://doi.org/10.3996/022014-JFWM-015.
- 79. Carvalho SB, Velo-Anton G, Tarroso P, Portela AP, Barata M, Carranza S, et al. Spatial conservation prioritization of biodiversity spanning the evolutionary continuum. Nature ecology & evolution. 2017; 1(6):1–8. pmid:28812637