Ecological correlates of extinction risk and persistence of direct-developing stream-dwelling frogs in Mesoamerica

dendrobatidis was the main driver of DSFS declines and of the current lack of recent observations of fifteen DSFS. We used our results to propose specific actions for all DSFS included in this study. This study highlights the importance of considering multiple threats and spatial scales when assessing the status of declined and threatened species.

Understanding the threats to wildlife across space and time is essential for developing effective conservation strategies.In Mesoamerica (i.e., the region that extends from Central Mexico to the most southern point in Panama) at least 40% of amphibian species declined between the late 1970s and the early 2000s.Most of these declines have been linked with the destruction of suitable habitats for amphibians as well as the introduction and spread of the fungal pathogen Batrachochytrium dendrobatidis, the causal agent of the disease known as chytridiomycosis.In this study, we quantified geographic and elevational ranges for direct-developing, stream-dwelling frog species (DSFS) in Mesoamerica.Within the range of each DSFS, we estimated the extent of suitable area for the occurrence of Batrachochytrium dendrobatidis and the types of land use to quantify habitat deterioration.To date, 33% of DSFS remain undetected since 2005 or before.At the regional level (i.e., Mesoamerica), as in previous studies, we found that narrow geographic and elevational ranges increased vulnerability to extinction.Nevertheless, the ranges of 83% of DSFS were composed of 50% or more high-quality habitats between the 1980s and 2005, when most species declined.We also found that on average, 80% of the range of each species currently overlaps with predicted suitable areas for the occurrence of Batrachochytrium dendrobatidis.At the local and site level (i.e., focusing on two species ranges where extensive monitoring has been conducted even before decline occurred), we found that the present suitable habitat for Batrachochytrium dendrobatidis corresponded with the reduction in predicted habitat suitability and climatic niche of DSFS.We also found that the location of remnant populations can be predicted by environmental factors, which can help identify regions where remnant populations of these declined species and others with similar ecology may occur.Combined, results from our regional and local analyses support the idea that Batrachochytrium dendrobatidis was the main driver of DSFS declines and of the current lack of recent observations of fifteen DSFS.We used our results to propose specific actions for all DSFS included in this study.This study highlights the importance of considering multiple threats and spatial scales when assessing the status of declined and threatened species.

Introduction
Environmental threats can vary greatly across space, are scale-dependent, and can determine species ranges (Isbell et al., 2017;Johnson et al., 2019).However, only a few studies have evaluated threats for entire vulnerable clades across spatial scales (Brittain et al., 2010;Stohlgreen et al., 2002).Among vertebrates, amphibians represent the most endangered class (Monastersky, 2014).According to the International Union for Conservation of Nature Red List of Threatened Species (hereafter 'IUCN Red List') over 56% of the approximately 8460 described species (Frost, 2021) are likely threatened by extinction (IUCN, 2019).Although threats to amphibians are considered individually, many amphibian species are facing multiple threats at the same time, which can increase extinction risk (Hof et al., 2011; but see Greenville et al., 2021).Thus, one of the main challenges scientists have is to identify the various threats to amphibians across landscapes and spatial scales, so that appropriate conservation actions can be implemented in regions where they will have the largest impact (Catenazzi, 2015).
Amphibians represent an ideal group on which to conduct threat assessments because species tend to exhibit discontinuous ranges across the landscape (Semlitsch and Bodie, 2003), and the threats to amphibians vary widely across spatial scales (Bosch et al., 2004;Grant et al., 2016).For example, the effects of introduced pathogens (Scheele et al., 2017) can be exacerbated by climate change, pollutants, and habitat deterioration (Becker et al., 2007;Rollins-Smith, 2017), and the range of pathogens and host species can be driven by processes that occur at different spatial scales (Chase et al., 2019;Keesing et al., 2010).

Amphibian declines in Mesoamerica
Mesoamerica (here defined as the geographic region that extends from Central Mexico to the most southern point in Panama) is a global hotspot of amphibian biodiversity and endemism (Savage, 2002;Wilson and Johnson, 2010) and is home to approximately 850 described amphibian species (AmphibiaWeb, 2020).However, the destruction and alteration of habitat have jeopardized many amphibian species in the region (Frías-Alvarez et al., 2010;Mayani-Parás et al., 2019;McCranie and Wilson, 2002;Savage, 2002).Furthermore, at least 40% of amphibian species have experienced declines even in seemingly undisturbed environments (Crawford et al., 2010;Whitfield et al., 2016).Many of these declines have been associated with the introduction and spread of epizootic lineages of the pathogenic fungus Batrachochytrium dendrobatidis (hereafter 'Bd ' Longcore et al., 1999;Rosenblum et al., 2013) between the late 1970s and early 2000s (Lips et al., 2008;Puschendorf et al., 2006a).During this Bd epizootic period (here defined as the period between the introduction of Bd in Mesoamerica -sometime in the mid-1970s or early 1980s-and 2005 when Bd became enzootic; Cheng et al., 2011), Bd rapidly spread across Mesoamerica causing massive die-offs of multiple species and even the collapse of entire amphibian communities (Crawford et al., 2010;Lips et al., 2006).Currently, several studies at the local (country) level (Bolom-Huet et al., 2019;Kilburn et al., 2010;Zumbado-Ulate et al., 2021a) suggest that Bd is now widespread and has not caused high rates of mortality since the early 2000s.For example, in a well-studied system in Panama, Bd seems to be endemic and persisting species may have evolved a degree of resistance to the pathogen that allows them to coexist with Bd, but perhaps not recover rapidly from infection (Voyles et al., 2018).In this study, we considered 2005 or later as 'Bd enzootic period' following spatiotemporal models that suggest that Bd was enzootic across Mesoamerica by 2005 (Cheng et al., 2011;Lips et al., 2008).Although Bd has been detected in all Mesoamerican countries, political, economic, and social reasons have historically prevented systematic sampling across the region and instead favored the collection of opportunistic data which ends up in searchable and shared databases in southern Mesoamerica (Costa Rica and Panama); this method of data collection is an important factor to be considered when using existing data to estimate threats across space.
In Mesoamerica, direct-developing, stream-dwelling frog species (hereafter 'DSFS') group within two phylogenetically related clades: the Craugastor punctariolus group (34 spp.) and the C. milesi group (12 spp.) (Hedges et al., 2008).All 46 DSFS are endemic to Mesoamerica and are regionally known as 'robber frogs' or 'stream frogs' (Campbell and Savage, 2000;Hedges et al., 2008).Once considered extremely common, DSFS catastrophically declined during the Bd epizootic period (Brem and Lips, 2008;Crawford et al., 2010;Lips et al., 2005Lips et al., , 2003;;Ryan et al., 2008;Whitfield et al., 2016).Currently, 90% of DSFS are considered threatened according to the IUCN Red List (IUCN, 2019).Similarly, the regional index of environmental vulnerability scores (EVS; Wilson and McCranie, 2004) shows that 98% of DSFS have the highest vulnerability to extinction (Wilson et al., 2013;Johnson et al., 2015).Each of the 46 DSFS has been monitored and studied since the 1960s by the authors of this study and other scientists, and this information has been included in the most recent IUCN Red List workshops conducted in Mesoamerica in 2018 and 2019 (Rodríguez et al., 2019).Results from this workshop included the absence of field records for fifteen of our study species since 2005 or before despite exhaustive sampling.Also, these workshops have served to highlight the rediscovery of several DSFS species, including some believed to be extinct (Chaves et al., 2014;Jiménez and Alvarado, 2017;Köhler et al., 2012;Kolby and McCranie, 2009;Kubicki, 2016).However, in most of these cases, rediscoveries are represented by the observation of one or few individuals, which greatly differs from the high abundance of DSFS observed during pre-decline times (Campbell and Savage, 2000;McCranie and Wilson, 2002;Savage, 2002).
In this study, we aimed to examine the effect of natural and human-induced threats (range size, climatic suitability for pathogens, and habitat deterioration) in the current range of DSFS.At the regional level (i.e., Mesoamerica), we first quantified the geographic and elevational ranges of each of the 46 DSFS contained within the C. punctariolus group and the C. milesi group.We also estimated the suitable habitat for the occurrence of Bd and quantified types of land use during the Bd epizootic and enzootic periods in Mesoamerica.Based on previous studies that have shown high susceptibility of DSFS to Bd, and a history of population decline in seemingly pristine environments (Lips et al., 2003;Puschendorf et al., 2006b;Ryan et al., 2008;Zumbado-Ulate et al., 2019b), we predicted that the occurrence of Bd was the main driver of the decline of many DSFS since 2005 or before.At the local level, we predicted that the habitat suitability and climatic niche of two specific focal DSFS (C.ranoides and C. taurus) contracted after the introduction of Bd.We also identified biotic and abiotic factors that could predict the occurrence of remnant populations.This work will help understand the causes of the assumed local absence of declined species.Furthermore, it will provide specific actions to inform which regions and species are most likely to experience future declines and extinctions.

Methods
Our spatial analyses were set to the current range of DSFS in Mesoamerica (N 5-40º, W 70-118º).We hold the resolution of our raster layers constant (30 arc-s) and modified the extent (i.e., scope ;Schneider, 2001).Our maps and geographic analyses were created with ArcMap 10.8 (ESRI® Redlands, CA, USA) and R 4.1.2(R Development Core Team, 2021).We used shapefiles from Atlas Digital Costa Rica 2014 (Ortiz-Malavasi, 2014) and the Database of Global Administrative Areas (GADM; https://gadm.org/data.html).In all our analyses we used the World Geodetic System datum (WGS84) as the reference coordinate system.

Regional-level analyses
For each of the 46 DSFS (Supplementary materials 1: Fig. S1 and S2; Table S1), we estimated range size and quantified land use and the predicted suitable area for Bd within each species range polygon (see section 3.1.1).

Range polygons
Updated species-range polygons for each of the 46 DSFS were provided by the IUCN Amphibian Specialist Group.The IUCN range polygons included the spatial information generated in the IUCN Red List workshops conducted in Mesoamerica in 2018 and 2019 (Rodríguez et al., 2019).Because IUCN has different classifications assigned to range polygons (Bland et al., 2017), we only considered the polygons corresponding to the extant and extinct range of our focal species and excluded those that represented regions of uncertain presence.

Bd occurrence dataset
To generate robust predictions of the climatic suitability of Bd across Mesoamerica, we exhaustively consulted peer-reviewed literature and compiled an initial presence-only dataset for Bd that included 181 occurrence points from Mesoamerica and 41 occurrence points from northern Colombia, southern United States, and Cuba (Supplementary materials 1: Fig. S3; Table S2).We cleaned our dataset by checking typos, removing unreferenced records, cross-checking geographic coordinates, and removing coordinates with a geographic inaccuracy > 1000 m and duplicates using the R packages 'scrubr' (Chamberlain, 2022) and 'Coor-dinateCleaner' (Zizka et al., 2019).To reduce the effect of spatial clustering caused by opportunistic sampling, especially in southern Mesoamerica (Costa Rica and Panama), we filtered and excluded the occurrence points that were separated by a distance > 15 km with the R package 'spThin' using 100 random repetitions (Aiello-Lammens et al., 2015).Our final dataset consisted of 133 filtered occurrence points.

Regional abiotic dataset
We downloaded the 19 BIOCLIM variables (Booth et al., 2014) from WorldClim v2.1 (Fick and Hijmans, 2017) and the 16 environmental raster layers for ecological modeling 'ENVIREM' (Title and Bemmels, 2018), both at a spatial resolution of 30 arc-s.The BIOCLIM variables have been widely used to predict the distribution of Bd, which facilitates the comparison of our predictions with other regional models (e.g., Bolom-Huet et al., 2019;Puschendorf et al., 2009;Zumbado-Ulate et al., 2021a).Similarly, the ENVIREM dataset has been shown to potentially improve the performance of habitat suitability maps that use BIOCLIM variables.Following Title and Bemmels (2018), we excluded the ENVIREM layers 'aridityIndexThornthwaite' (which is redundant with the 'climaticMoistur-eIndex'), and 'monthCountByTemp10' (which is categorical).For data extraction, we cropped the 33 environmental layers with the R package 'raster' (Hijmans et al., 2015).We used the R package 'usdm' (Naimi, 2015) to detect collinearity among predictors by quantifying and selecting the predictors which had a variance inflation factor < 3 and pair-wise correlations < 0.6.We retained eight environmental predictors in our abiotic dataset: 'isothermality', 'maximum temperature of warmest month', 'precipitation seasonality', 'precipitation of warmest quarter', 'precipitation of coldest quarter', 'mean monthly potential evapotranspiration of driest quarter', 'mean monthly potential evapotranspiration of warmest quarter', and 'sum of mean monthly temperature for months with mean temperature greater than 5 ℃ multiplied by number of days' (Supplementary materials 1: Table S3).

Species range size
Narrow-ranged species have been associated with an increased risk of extinction compared with wide-ranged species (Gaston and Fuller, 2009).Therefore, we were interested in determining whether narrow ranges (geographic and elevational ranges) correlated with those species lacking records in the field since 2005 or before (i.e., species that may be extinct).We used the 'calculate geometry' GIS tool to estimate the area in km 2 contained within the minimum convex polygon of each species-specific IUCN range polygon (Nilsen et al., 2008) as a proxy of the geographic range.The elevational range inhabited by each species was defined as the difference between the upper and lower elevation limits (e.g., White and Bennett, 2015;Yu et al., 2017).The elevation limits for each of the 46 DSFS species were obtained from the IUCN Red List website (https://www.iucnredlist.org).

Habitat suitability of Bd in Mesoamerica
We were interested in identifying if declined species' ranges substantially overlapped with the predicted suitable habitat for Bd in Mesoamerica compared to species that persisted.For this, we built a habitat suitability map of Bd in Mesoamerica during the enzootic period (2005 and after) and estimated the extent of suitable habitat (hereafter 'ESH') (Brooks et al., 2019) for Bd within the range polygons of each of the 46 DSFS.We were not able to predict the range of Bd during the epizootic period (mid-1970s to 2005) due to the lack of uniform sampling in most of Mesoamerica and global methodical limitations to detect low levels of Bd before the 2000s.However, evidence shows that Bd spread southwards in Mesoamerica during the Bd epizootic period, especially across mountain ecosystems (Cheng et al., 2011;Lips et al., 2008;Phillips and Puschendorf, 2013).Although other habitat suitability maps have predicted the climatic suitability of Bd in Mesoamerica (e.g., James et al., 2015;Ron, 2005), we have included a more comprehensive and updated collection of occurrence points than previous models, used more robust methods to improve the accuracy of our predictions, and followed strict criteria to control spatial clustering of species records caused by unbalanced sampling, spatial autocorrelation, and multicollinearity (see below).
Our predictions of suitability were generated with the MaxEnt algorithm (Phillips et al., 2006) implemented through the R package 'ENMeval' (Muscarella et al., 2014).We obtained 36 candidate models using the Bd occurrence dataset (section 3.2.2), the regional abiotic dataset (section 3.2.3),and the following settings: partition method = 'block'; random points = 10000; algorithm = maxent.jar;regularization of multiplier values = 0.5-3 with increments of 0.5; feature classes = L, LQ, H, LQH, LQHP, LQHPT.We selected the model with the highest value of average test of the area under the receiver operating characteristic curve (AUC mean).If two or more candidate models had the same AUC mean value, we selected the model with the lowest omission rate at minimum training presence (orMTP; Radosavljevic and Anderson, 2014).
Our model predictions were transformed into raster objects with the R package 'dismo' (Hijmans et al., 2021).With the logistic output of our habitat suitability maps, we generated absence-presence maps of the potential range of Bd in Mesoamerica using the equal training sensitivity and specificity threshold (Liu et al., 2005).We transformed our predictive absence-presence maps of Bd into shapefiles and calculated the ESH of Bd in Mesoamerica in square kilometers (km 2 ) using a projected coordinate system (WGS84).

Suitable habitat for DSFS
Habitat loss has been associated with increased extinction risk in amphibians (Becker et al., 2007).We aimed to determine if suitable habitats for declined DSFS (i.e., old-growth and secondary forest ecosystems) substantially decreased from the Bd epizootic to enzootic periods.We estimated the land use within the range polygon of each of the 46 DSFS during both the Bd epizootic and enzootic periods.For this, we downloaded 'The Human Footprint 2018 Release' (Sanderson et al., 2002;Venter et al., 2016) at a spatial resolution of 30 arc-s.This dataset consists of two raster layers that quantify the cumulative human pressure on terrestrial ecosystems in 1993 (within the Bd epizootic period) and 2009 (within the Bd enzootic period).The human pressure is measured using eight variables: built-up environments, population density, electric power infrastructure, crop lands, pasture lands, roads, railways, and navigable waterways (Venter et al., 2016).In each raster, scores ranging from 0 to 50 are assigned to each cell, with the highest values representing cells with the highest human footprint.To build an accurate categorical representation of the land-use in Mesoamerica, we superimposed land-use layers of Central America (Ortiz-Malavasi, 2014) on both human footprint layers and reclassified the raster values into four categories according to each cell score: suitable habitat for DSFS (old-growth and secondary forest ecosystems, 0-10), intermediate disturbance (small forest patches, croplands, and pastures; 11-20), rural development (21− 30), and urban development (31− 50).We transformed both reclassified (i.e., categorized) 1993 and 2009 human footprint layers into vector data (i.e., shapefiles).Then, we overlaid both shapefiles to quantify changes in land use during these 16 years (Supplementary materials 1: Fig. S4).

Assessing threat risk
We conducted multiple logistic regression using stepwise, backward selection (Hosmer et al., 2013) to examine the current known occurrence of focal DSFS using binomial response variables (detected or undetected).'Undetected' species were all species lacking records in the field since 2005 or before.Here, we are not assuming undetected species as extinct, however, all focal DSFS have been extensively searched for in the last four decades (IUCN, 2019), so we are certain that those species have declined and remain undetectable throughout their known range.The year 2005 was set according to models that suggest that Bd likely became enzootic across Mesoamerica by 2005 (Cheng et al., 2011;Lips et al., 2008).To deal with quasi-perfect separation we followed a Bayesian analysis with non-informative prior assumptions (Gelman et al., 2008) using the R package 'arm' (Gelman et al., 2021).Output models were ranked according to Akaike's information criterion (Burnham and Anderson, 2004).
Most of our predictors consisted of continuous variables to improve model performance and avoid arbitrary categorization (Hosmer et al., 2013;Worster et al., 2007).We used 1) the area of minimum convex polygon, 2) the elevational range inhabited by a species, 3) the percent of overlap of Bd with the range of each DSFS, 4) the percent of suitable habitat for DSFS within the range of each focal species during the Bd epizootic period, and 5) clade (C.punctariolus or C. milesi group).Before running the logistic regression, we ran a correlation test between continuous predictors to identify high multicollinearity (pair-wise correlations > 0.6).Since all pair-wise correlations were below our threshold value we did not exclude any predictor.

Local, species-range level analyses
We wanted to quantify historical niche dynamics of Bd within the specific range of DSFS.Since these analyses require extensive long-term field data including detection of Bd across the species' ranges, we were only able to use data from two species that fully meet these criteria: the 'dry forest robber frog' (C.ranoides), and the 'Golfito robber frog' (C.taurus).Therefore, we used a bounding box that included Costa Rica, and neighboring areas in Nicaragua and Panama (N 7-14º, W 80-88º). Due to the considerably large amount of detail needed to describe our methods, we have created a separate document (Supplementary materials 2) to ensure the repeatability of our approaches.

Focal species
The 'dry forest robber frog' is a subspecies of C. ranoides (Supplementary materials 2: Fig. S1a) that was historically distributed across the Santa Elena Peninsula and the Guanacaste Volcanic range in Costa Rica; although its historical distribution could be wider (Puschendorf et al., 2019).Currently, this species is known only from scattered remnant populations across the Santa Elena Peninsula (Puschendorf et al., 2019(Puschendorf et al., , 2009;;Sasa and Solórzano, 1995;Zumbado-Ulate and Willink, 2011).The 'Golfito robber frog,' C. taurus (Supplementary materials 2: Fig. S1b), was historically distributed from Dominical in central Pacific Costa Rica to the Burica Peninsula, a region that includes extreme southeastern Costa Rica and the eastern part of the Chiriquí Province of Panama.Currently, it is only known in a small portion of the Burica Peninsula (Chaves et al., 2014).Extensive search efforts have been conducted across its range, especially in the localities of Rincón de Osa, Golfito, and along 47 streams and tributaries across Paso de la Danta Biological Corridor, where it remains undetected (Chaves et al., 2014;Zumbado-Ulate et al., 2021b).

DSFS dataset
We compiled an initial presence-only dataset for both focal species (historical dataset) using occurrence data from 1960 to 2020 obtained from the Museo de Zoología at Universidad de Costa Rica, the Global Biodiversity Information Facility (GBIF.org;https://doi.org/10.15468/dl.c7yupyfor C. ranoides and https://doi.org/10.15468/dl.dnmjrdfor C. taurus), and new occurrence points generated from expeditions conducted during this study.Then, we created a new dataset (present dataset) by including only the observations recorded during the Bd enzootic period.In total, we compiled 88 records for the dry forest robber frog (30 observations during the Bd enzootic period) and 286 records for the Golfito robber frog (19 observations during the Bd enzootic period).We used the same methods described in section 3.2.2 to clean both datasets independently; but this time we used a spatial filter of 1 km to reduce spatial clustering, matching the resolution of the climatic layers of the 'local abiotic dataset' (section 3.6.2.).Our final DSFS dataset consisted of 124 occurrence points: 31 for the dry forest robber frog (23 used to model the historical range and eight to model the present range during the Bd enzootic period) and 93 for the Golfito robber frog (79 to model the historical range and 14 to model the present range during the Bd enzootic period).

Local abiotic dataset
To select environmental predictors at the local level, we followed the same methods described in section 3.2.3.Our final local abiotic dataset consisted of six environmental predictors: isothermality, temperature annual range, precipitation of the wettest month, precipitation of driest month, minimum temperature of the warmest month, and mean monthly potential evapotranspiration of driest quarter (Supplementary materials 2: Table S1).

Local environment
Although our two focal species have declined in most of their range, they persist in remnant population in some stream networks of Costa Rica (Zumbado-Ulate et al., 2019a, 2014).Therefore, we characterized the biotic and abiotic environment in sites with and without remnant populations to identify environmental conditions that might predict the location of remnant populations of DSFS and other species with similar ecology.

Fieldwork
We conducted six comparative field surveys on six stream networks in six regions of Costa Rica: the Santa Elena Peninsula, Rincón de la Vieja Volcano, Uvita, Rincón de Osa, Golfito, and the Punta Burica Peninsula during the dry seasons of 2017 and 2018 (Supplementary materials 2: Fig. S2).Our field trips were scheduled during the dry season to maximize the collection of biotic and abiotic data.During the rainy season, tropical stream networks highly increase in volume and flow speed, reducing the detectability of DSFS and making sampling extremely dangerous (Savage, 2002).

Survey transects
In each study stream network, we selected one or two streams where the two focal species were common before the Bd epizootic period and surveyed 25 m linear transects to search for the focal species and to collect environmental data.The total number of transects varied according to the accessibility of each study stream (Supplementary materials 2: Table S2).Transects were surveyed twice a day, in the morning starting at 9 am and in the evening starting at 6 pm.The duration of each survey varied according to the amount of data collected.

Biotic environment
High community heterogeneity has been associated with suppression of pathogen transmission in susceptible species (i.e., dilution effect Schmidt and Ostfeld, 2001).Therefore, we aimed to determine if remnant populations of both focal species of DSFS occur in amphibian communities with high heterogeneity.For this, we counted individuals of all amphibian species that were visually or acoustically detected across the transects, in both morning and evening surveys (Supplementary materials 2: Table S3 and S4).For acoustic detection, we searched for individuals calling from rocks or vegetation within the transect.Once an individual was found, we confirmed the species and added the data to our counts.Using the data from our counts, we estimated the community diversity (H) using Hill numbers (q = 1; Norris and Pollock, 1998;Chao et al., 2014) and the community heterogeneity (J) using the Shannon evenness index (Krebs, 1989;Norris and Pollock, 1998).We also characterized abiotically each transect by measuring the following variables in its midpoint: stream width and depth (as a proxy of water volume in m 3 by multiplying width x depth x 25), flow speed (m/s), percent of canopy cover, air and water temperature ( • C), percent of relative humidity, water pH, and elevation (m) (Supplementary materials 2).

Reduced habitat suitability of DSFS
We aimed to identify if geographic range contraction occurred for both focal DSFS after the introduction of Bd.Therefore, we generated habitat suitability maps for both focal DSFS with two scenarios: 1) 'no decline' (using all collection of occurrence points) which assumes no local extinctions and generates the historical habitat suitability; and 2) 'decline' (i.e., known present range from 2005 which matches the Bd enzootic period).To generate our predictions, we used the DSFS dataset (section 3.6.1),the local abiotic dataset (section 3.6.2),and the same methods used to predict the habitat suitability of Bd in Mesoamerica (section 3.3.2).
However, we used the 'N-1 Jackknife' as the method of data partition in three of our sets of predictions, aiming to maximize the limited available information when species occurrences < 25 (Muscarella et al., 2014).The block method was used to predict the suitability of the Golfito robber frog with the 'no decline' scenario (occurrences = 79).For both focal species, we generated absence-presence maps and calculated the ESH for both scenarios as described in section 3.3.2.

Niche contraction
We aimed to determine if climatic niche contraction occurred for both focal DSFS after the introduction of Bd.We used the R package 'humboldt' to convert the geographic space of both focal DSFS into environmental space (Supplementary materials 2).Then, we estimated the density of occurrences as probability distributions defined over the multivariate climatic niche (hereafter 'Z d ') (Broennimann et al., 2012(Broennimann et al., , 2007;;Petitpierre et al., 2012) with the R package 'ecospat' (Di Cola et al., 2017).We estimated Z d for the historical and present distributions of the two focal DSFS and Bd during the enzootic period.We were not able to obtain Z d for epizootic Bd because the number of occurrence points was < 5 within the focal study species' ranges before 2005.
We estimated niche overlap with the Schoener's D overlap index (Schoener, 1968), which varies between 0 (no overlap) and 1 (identical niches).We predicted low similarity between the historical and present distributions of both focal DSFS.We used a niche similarity test to quantify niche divergence.The niche similarity test also reports a D value that varies between 0 (different niches) and 1 (identical niches).We predicted niche divergence after the introduction of Bd for the two focal DSFS.We also estimated niche unfilling (hereafter 'contraction'), niche expansion, and niche stability (Warren et al., 2008).We predicted that 1) the climatic niche of both focal species contracted after the introduction of Bd, 2) the climatic niche of enzootic Bd expanded within the environmental

Table 1
Range area and elevational distribution of 46 species of direct-developing, stream-dwelling frogs in Mesoamerica.For each species, we present the area of minimum polygon convex (A MCP ) as a proxy of the geographic area, the lower and upper elevation limits, and the percent of a species range that overlaps with 1) enzootic Batrachochytrium dendrobatidis (Bd), 2) suitable habitat (i.e., old-growth and secondary forests; SH), 3) intermediate disturbance (ID), 4) rural development (RD), and 5) urban development (UD) during the Bd epizootic period (Ep) and the Bd enzootic period (En).space of both focal species, and 3) the existence of a reduced environmental space where enzootic Bd and current distributions of both focal DSFS co-occur, suggesting the existence of environmental refugia where Bd and DSFS can coexist.All observed distributions from the niche equivalency test were compared with null distributions generated from 100 randomly simulated values to determine if the observed niche differences are larger than those expected when comparing random niche distributions.All calculations were implemented with 'ecospat'.

Local environmental conditions
To identify the climatic conditions that might predict the occurrence of remnant populations of focal DSFS, we generated buffers (radius = 10 km) around each of the six study stream networks and generated 500 random points with no replacement within each buffer.For each point, we extracted values for all six environmental predictors of the local abiotic dataset (section 3.6.2) to achieve a full climatic description of the six study stream networks.722 random points were projected on the ocean surface (pixels with missing values); therefore they were automatically excluded from the analysis.We compared networks using a principal component analysis (PCA) generated with the R package 'FactoMineR' (Lê et al., 2008).The results of the PCA were interpreted and visualized with the R package 'FactoExtra' (Kassambara and Mundt, 2017).Similarly, to identify microhabitat conditions that might favor the occurrence of remnant populations of focal DSFS, we first quantified collinearity among stream predictors as shown in section 3.2.2.We retained six stream predictors (elevation, water pH, volume, flow speed, percent of canopy cover, and community heterogeneity; Supplementary materials 2: Table S2).We also compared the study networks using a principal component analysis (PCA) as described earlier in this section.
Fig. 2. a) Habitat suitability map of two focal direct-developing, stream-dwelling species assuming no decline (using all historical collection of occurrences.The upper left polygon shows the predicted historical climatic suitability for the dry forest robber frog (subspecies of Craugastor ranoides).The bottom right polygon displays the predicted historical climatic suitability for the Golfito robber frog (C.taurus); b) the climatic suitability of both focal species decreased assuming the current extinction of undetected populations; c) under a no decline scenario, the extent of suitable (ESH) habitat is approximately 2600 Km 2 for the dry forest robber frog and 6200 Km 2 for the Golfito robber frog; d) assuming current extinction of declined populations, the ESH habitat has reduced in 85% for the dry forest robber frog and 95% for the Golfito robber frog.

Regional-level analyses
Overall, 46% of DSFS had geographic ranges smaller than 1000 km 2 , with 21% between 1000 and 5000 km 2 , and 33% larger than 5000 km 2 (max.C. vocalis 170,000 km 2 ).The elevational range inhabited by a species was greater than 1000 m for 17 DSFS, but only seven of those species had an elevational range greater than 1500 m (max.C. vocalis 2200 m).The habitat suitability maps for Bd in Mesoamerica derived from our most robust model (FC = LQ, RM = 3, AUC= 0.82; OrMTP = 0.05; Fig. 1) showed high suitability for Bd across the Sierra Madre Oriental and the Sierra Madre Occidental range in México, the Sierra de las Minas in Guatemala, and most of Costa Rica and Panamá.Intermediate suitability was predicted across the lowlands of Nicaragua, Honduras, and southwestern México (Fig. 1a).The environmental predictors with the highest permutation (Supplementary materials 1: Table S3) were 'max temperature of warmest month' (74.6%), and 'precipitation of warmest quarter' (14.4%), and suggested the most suitable areas for Bd occur in areas with high precipitation (> 2000 mm) and temperatures below 25 • C.
The ESH of Bd in Mesoamerica during the enzootic period was approximately 850,000 km 2 (41.6% of Mesoamerica; Fig. 1b).On average, 80% of the range of each species overlapped with the predicted suitable habitat for Bd.This overlap was greater than 90% in 31 species, 15 of which have been undetected since 2005 or before (Table 1).Furthermore, 72% of suitable area for Bd during the enzootic period occurred in old-growth and secondary forests, which is considered suitable habitat for DSFS, and 23% in intermediate disturbance habitats, which may be inhabited by some DSFS (e.g., C. amniscola).
In our combined threats assessment, our most robust model (Supplementary materials 1: Table S4) fit significantly better than an empty model (X 2 = 29.9,df = 3, p < 0.001) and showed that current species undetectability significantly increased with Bd overlap and clade (the C. milesi group is more vulnerable).

Local-level analyses
The habitat suitability maps (Fig. 2) derived from our most robust models (Supplementary materials 2: Table S5) showed that the current ESH of the dry forest robber frog reduced from 2597.5 km 2 to 397.4 km 2 (85% of the historical suitable range) after the introduction of Bd.The suitable habitat for this species was restricted to the tropical dry forest of the Santa Elena Peninsula in Costa Rica.Similarly, the ESH of the Golfito robber frog decreased from 6200 km 2 to 328.2 km 2 (95% of the historical suitable range) under a Bd-driven scenario.The suitable habitat for this species was limited to the Burica Peninsula in Costa Rica and Panama.The predictors with the highest permutation importance in our habitat suitability maps (Supplementary materials 2: Table S1) were 'precipitation of driest month' for the dry forest robber frog and 'mean monthly potential evapotranspiration of driest quarter' for the Golfito robber frog.These specific predictor contributions suggested that the most suitable habitats for both focal DSFS species are restricted to dry and warm regions, even during the rainy season when Bd prevalence might increase (Whitfield et al., 2017).
Our estimations of niche overlap between the historical and present distributions of the dry forest robber frog and the Golfito robber frog (D = 0.09 and D = 0.23, respectively) and niche similarity (D = 0.05 and D = 0.15, respectively) suggested the current niche of the dry forest robber frog has diverged from the historical niche and the centroid is now oriented towards dry conditions.The observed differences in niche dynamics were higher than expected if compared with random niches only for the dry forest robber frog (P = 0.04) (Supplementary materials 2: Table S6).We also found a strong contraction in the climatic niche of both focal DSFS after the spread of Bd.For the dry forest robber frog, the climatic niche contracted 91% when compared with the historical climatic niche (Fig. 3a).This contraction corresponded with 91% of niche expansion of Bd within the climatic niche of the dry forest robber frog.The current shared environmental space between the dry forest robber frog and Bd was estimated at 8.5% (Fig. 3b).Similarly, the Golfito robber frog experienced a reduction of 87% in the climatic niche (Fig. 3c), which also coincided with 76% of expansion of Bd within the climatic niche of the Golfito robber.The current shared environmental space between the dry forest robber frog and Bd was estimated at 12% (Fig. 3d).
The highest loadings of our local environmental predictors in our PCAs (Fig. 4; Supplementary materials 2: Table S7 and Table S8) suggested that the location of remnant populations of focal DSFS (temperature annual range = 0.82 in PC1; minimum temperature of the warmest month = 0.86 and the mean monthly potential evapotranspiration of the driest quarter = 0.75 in PC2) were associated with dry and warm climatic conditions.Regarding the microhabitat predictors, our second PCA suggested that remnant populations of focal DSFS were associated with high community heterogeneity (PC1, − 0.91) in the Santa Elena Peninsula, and water pH > 8 (PC2, 0.91) in the Punta Burica Peninsula.

Discussion
Amphibian conservation is a global priority, especially after the accelerated global decline of amphibian populations during the last 40 years (Mendelson et al., 2019).Therefore, it is essential to measure the impact of anthropogenic threats across the landscape (Cohen et al., 2016) to develop specific management plans for threatened species.Here, we focused on evaluating the threat risk and current known occurrence of 46 DSFS.For our analyses, we used a combination of published data and new data derived from fieldwork conducted across Mesoamerica.Our study is innovative and relevant to other amphibian clades, especially direct-developing and stream-associated amphibians, and can be replicated in different taxa.At the regional level (Mesoamerica), we found evidence that suggests that Bd was a primary cause of the decline of most DSFS during the Bd epizootic period (Supplementary materials 1: Table S4).We found that Bd is predicted to be widespread across the historical and present ranges of most DSFS (Fig. 2).Additionally, our results suggest that the species within the Craugastor milesi group are more vulnerable to extinction than the species within the C. punctariolus group, which might be linked to their narrow range size or genetic effects.Species with narrow range sizes (geographic and elevational) have been commonly associated with low competitive ability, low aptness to adapt to novel environments, and low capacity to tolerate environmental variability (Gaston and Fuller, 2009;Lawton, 1993;White and Bennett, 2015).Similarly, low genetic diversity has been linked to reduced evolutionary potential, and elevated extinction risk (Spielman et al., 2004).At the local level, our results suggest that the spread of Bd has caused contraction of the geographic range and climatic niche of two focal DSFS, and that variation in biotic and abiotic conditions within their original ranges has allowed remnant populations to coexist with this pathogen in some habitats (Fig. 3 and Fig. 4).Our findings suggest that using different spatial scales when conducting threat assessments could allow scientists to obtain a more accurate picture of the threats and status of endangered taxa.
We found that the geographic ranges of 81% of the species within the C. punctariolus group and 91% within the C. milesi group contained at least 50% of suitable habitat for DSFS (old-growth and secondary forest cover; Table 1) during the Bd epizootic period.These findings suggest that habitat loss was unlikely to be the primary driver of declines in DSFS.However, the increase of 65% of urban development in Mesoamerica between the Bd epizootic and enzootic periods (approximately 20,000 km 2 , Fig. 2) might also have affected populations that historically inhabited areas close to human settlements (Campbell and Savage, 2000;McCranie and Wilson, 2002), and some changes in habitat quality may be responsible for range contractions.
Our analyses at the local level allowed us to assess the persistence of two focal DSFS using methods that are less accurate at a larger scale.We found that the spread of Bd during the enzootic period coincided spatially and temporally with the climatic niche contraction and the extinction of most populations of the dry forest robber frog and the Golfito robber frog across their historical ranges (Fig. 3, Supplementary materials 2: Table S6).Our suitability and climatic niche predictions suggest that seasonal dry and semidry conditions have the largest impact on the present suitable habitat of both DSFS (Fig. 2 and Fig. 3).This finding coincides with other studies that have suggested that the persistence of remnant populations of C. ranoides and C. taurus is associated with dry and warm regions that constrain the growth of Bd (Chaves et al., 2014;Granados-Martínez et al., 2021;Puschendorf et al., 2009;Zumbado-Ulate et al., 2014).Here we provide new evidence that high community heterogeneity (e.g., Schmidt and Ostfeld, 2001;Searle et al., 2011) and abiotic conditions outside of the optimal range for Bd transmission (e.g., pH 6-8; Piotrowski et al., 2004) could also be suppressing pathogen transmission in locations where remnant populations exist (Fig. 4, Supplementary materials 2: Table S8).Future field and laboratory studies are needed to assess if these environmental factors alone or in combination with other factors are constraining the spread of Bd in environmental refugia from decline (Puschendorf et al., 2011).

Conservation implications
We consider the conservation of this neglected group of endangered amphibians a regional priority; therefore, we have created a separate document (Supplementary materials 3) that reports the species' rediscoveries and provides specific conservation actions for each of the 46 focal species.These actions have been reviewed by conservation experts, including those within the IUCN Amphibian Specialist Group.We have included general conservation actions (e.g., long-term monitoring of remnant populations in specific locations, threat quantification, and mitigation) as well as highly specific actions depending on the status of each species (e.g., Ex-situ reproduction, expeditions to specific targeted regions).The recent rediscoveries of several DSFS are a beacon of hope that other species are surviving and suggest partial recovery and host-pathogen coexistence (e.g., Whitfield et al., 2017).Thus, new studies are needed to determine if these rediscoveries are linked to environmental factors that reduce pathogen transmission (Puschendorf et al., 2011;Scheele et al., 2017) or eco-evolutionary rescue (i.e., rapidly evolved, reduced susceptibility to Bd) (Christie and Searle, 2018;DiRenzo et al., 2018).Here we report the rediscovery of a remnant population of C. amniscola at the state park of La Pera, in Chiapas, México (Supplementary materials 1: Fig. S2).

Conclusions
Results of this work can be used to make field inventories more effective by identifying priority areas for locating remnant, or even undiscovered, populations of DSFS.However, because a remnant population cannot be assumed to represent a recovering population (Mendelson et al., 2019), a long-term commitment to monitoring and conservation interventions is required to ensure persistence.For example, management actions that create or preserve habitat refugia from chytridiomycosis and target other threatening processes such as habitat loss, have great potential for supporting species' recovery (Skerratt et al., 2016).These conservation efforts will become increasingly important as future shifts in environmental conditions or other emerging threats may trigger a re-emergence of Bd, resulting in further declines or extinctions (Scheele et al., 2017).Furthermore, this study can be used as a reference for the evaluation of other direct-developing, stream-dwelling species, or even other Neotropical amphibian groups that were also decimated by Bd (e.g., Atelopus spp., Isthmohyla spp.).
Studies that assess threat risk across different spatial scales are needed to develop more successful conservation strategies and target species that need urgent attention.Our fine-resolution analyses were effective for assessing threats to species that have relatively small ranges (Schwartz, 1999;Wiens, 1989).The results of this work can be integrated with important conservation tools, such as the IUCN Red List, that serve as valuable inputs into setting conservation priorities.Considering the key role that the IUCN Red List plays at the global and regional levels, it is imperative that the extinction risk assessments and their underlying data are updated regularly.The expanded understanding of disease-related declines in DSFS presented in this study can be used to improve the quality of data in IUCN Red List assessments for these species, many of which have been undetected for long periods and for which data are scarce.

Fig. 1 .
Fig. 1. a) Map showing gradients of climatic suitability for Batrachochytrium dendrobatidis (Bd) in Mesoamerica; b) the extent of suitable area for Bd (beige region) comprises 41.6% of the total area of Mesoamerica.

Fig. 3 .
Fig. 3. Niche dynamics of the dry forest robber frog Craugastor ranoides and the Golfito robber frog (C.taurus) assuming extinction of multiple populations driven by Batrachochytrium dendrobatidis (Bd) during the Bd enzootic period (i.e., based on post-2005 records).White areas represent the full environmental space of the dry forest robber frog and the Golfito robber frog.a) The current climatic niche of C. ranoides (light gray area) has contracted 91% compared to the historical niche (blue area).b) The climatic niche of Bd (beige area) expanded 91% within the environmental space of C. ranoides.c) The current climatic niche of C. taurus (light gray area) has contracted 76% compared to the historical niche (blue area).d) The climatic niche of Bd (beige area) expanded 87% within the environmental space of C. taurus.The blue region shows the climatic space occupied by C. taurus where Bd is predicted to be absent.For both focal species, the centroid has shifted in position and direction towards dry environmental conditions.Light blue dots show the spatial location of extinct populations and black dots the location of remnant populations.In Figs.3a and 3c, the light blue dots show the spatial location of extinct populations, and black dots the location of remnant populations.In Figs.3b and 3d, the purple dots show the spatial location of sites where Bd has been found and black dots show the location of remnant populations of DSFS.