Anticipating the potential impacts of Batrachochytrium salamandrivorans on Neotropical salamander diversity

Emergent infectious disease caused by the fungal pathogens Batrachochytrium dendrobatidis (Bd) and Batrachochytrium salamandrivorans (Bsal) represents one of the major causes of biodiversity loss in amphibians. While Bd has affected amphibians worldwide, Bsal remains restricted to Asia and Europe, but also could be a major threat for salamanders in the Western hemisphere, including the 320 bolitoglossine species described. Here, we predict the suitable areas for Bsal in the Neotropics and assess its potential impact on bolitoglossine diversity. For this, we determined the geographic patterns of taxonomic, phylogenetic, and functional diversity for bolitoglossines and modeled the potential distribution of Bsal in the Neotropics. We identified which species and regions could be at risk from an eventual introduction of Bsal in the region, quantified the degree of overlap between regions of high diversity and the suitable conditions for the pathogen, and considered species IUCN Red List status, and geographic range size. We found that regions of high taxonomic, phylogenetic, and functional diversity are concentrated in the Trans‐Mexican Volcanic Belt, Sierra Madre Oriental, the southern portion of Sierra Madre del Sur and the mountains of Oaxaca in México, as well as the Chiapan‐Guatemalan highlands, and the Cordilleras of Costa Rica and Panama. Alarmingly, the regions of high diversity for bolitoglossines and over 75% of the ranges of the more threatened species could be affected by Bsal. Given the unknown vulnerability of these species, we strongly recommend measures to avoid the introduction of Bsal in the continent.


| INTRODUC TI ON
Amphibians are one of the most threatened groups of vertebrates (Monastersky, 2014), with approximately 56% of described species experiencing population declines and even extinction in the last 40 years (Stuart et al., 2008;Wake & Vredenburg, 2008). One widely studied cause of amphibian declines is chytridiomycosis (Berger et al., 1999), an emergent infectious cutaneous disease caused by the fungal pathogens Batrachochytrium dendrobatidis (Bd) (Longcore et al., 1999) and Batrachochytrium salamandrivorans (Bsal) (Martel et al., 2013). While Bd disrupts osmoregulation and eventually leads to death (Van Rooij et al., 2015;Voyles et al., 2009), Bsal erodes skin, promoting secondary bacterial overgrowth that may result in lethal septicemia (Bletz et al., 2018). The presence of Bd has been linked to drastic population declines and extinctions of many amphibian species around the globe (Lambert et al., 2020;Scheele et al., 2019); meanwhile, the recent emergence of Bsal has been associated with rapid population collapses in European salamanders (Martel et al., 2014;Stegen et al., 2017).
Current evidence suggests that both fungal pathogens have their origin in East Asia (O'Hanlon et al., 2018). While Bd has been found in all regions of the globe where amphibians occur (Lips, 2016;Olson et al., 2013), the current known distribution of Bsal remains restricted to Asia and Europe (Beukema et al., 2018;Laking et al., 2017;Lötters et al., 2020;Martel et al., 2014;Yuan et al., 2018). The detection of Bsal in the pet trade suggests that this activity might be the primary cause for its spread to naïve regions and continents (Cunningham et al., 2015;Gray et al., 2015;Laking et al., 2017;Martel et al., 2014Martel et al., , 2020Nguyen et al., 2017;Yuan et al., 2018).
Additionally, experimental studies have demonstrated that Bsal infections are lethal to several North American salamander species from the Plethodontidae and Salamandridae families and cause more erosive skin lesions than Bd Friday et al., 2020;Martel et al., 2014;North American Bsal Task Force, 2020). On the contrary, anurans might be more tolerant or resistant to Bsal, potentially acting as Bsal-reservoirs species (Martel et al., 2014;Stegen et al., 2017). Although some amphibian species may not become infected with Bsal and the knowledge on species-specific vulnerability to Bsal is still limited, the current evidence suggests that the introduction of Bsal in naïve ecosystems could be devastating, even in populations resistant or tolerant to Bd (Longo et al., 2019;Martel et al., 2014;Stegen et al., 2017). Therefore, the potential introduction and spread of this pathogen must be considered a latent threat for salamander-rich regions in continents where this pathogen is not currently present or still undetected (e.g., the Western hemisphere).  et al., 2000;Hanken & Wake, 1994;Parra-Olea et al., 2020). Their highly restricted distributions, together with significant ongoing threats, make species within this group highly prone to extinction (Ripple et al., 2017), with nearly 60% of its species categorized as Endangered (EN) or Critically Endangered (CR) on the International Union for Conservation of Nature Red List of Threatened Species™ (hereafter IUCN Red List of Threatened Species™; IUCN, 2021).
Currently, much of our knowledge of biodiversity is derived from analyses focused on the biological richness at the species level-that is, taxonomic diversity (TD) (Jarzyna & Jetz, 2016). Although important, TD measures alone ignore species relatedness and ecological functions, limiting the understanding of the underlying processes of diversity patterns (Cardoso et al., 2014). In recent years, this speciesbased approach has been complemented with the quantification of the phylogenetic diversity (PD) and functional (FD) dimensions of biodiversity. While PD provides insights on the evolutionary history of the biotic assemblage (Faith, 2013), FD captures the diversity of relevant traits for ecosystem functioning (Petchey & Gaston, 2006).
Efforts linking such dimensions have recently improved our understanding of the eco-evolutionary processes underlying the conformation of the rich amphibian biota in the Neotropics (Ochoa-Ochoa et al., 2020). A further required step is anticipating the potential impacts of contemporary threats (e.g., climate change, biological invasion, and emerging diseases) on such dimensions to develop effective and integrative conservation strategies.
In the specific context of emergent disease, identifying high-risk areas for the establishment of a pathogen and determining the most vulnerable host groups are essential to implement effective conservation strategies (Gray et al., 2015;Mendelson et al., 2019;Woodhams et al., 2011). Addressing these tasks from a framework that considers the multiple facets of salamander biodiversity and the potential areas the pathogen may colonize could ideally prevent or at least mitigate another collapse of the region's amphibian fauna. Here, we assessed the potential impact of Bsal outbreaks on bolitoglossines using two different approaches: 1. Assemblage-based analyses to formally describe the geographic patterns of different facets of bolitoglossine biodiversity (i.e., taxonomic, phylogenetic, and functional) and determine their potential overlap with suitable areas for Bsal, and 2. Clade-based methods to determine which species could be at risk from an eventual introduction of Bsal in the region. Using species distribution modeling, phylogenetic approaches, and comparative methods, we identified both the geographic regions where most species occur and the portions of the bolitoglossine phylogeny that might be most drastically pruned. The incorporation of this eco-evolutionary perspective aims to better guide conservation priorities for these Neotropical species.

| Geographic patterns of bolitoglossine biodiversity
We compiled range maps for 304 bolitoglossine species (~95% of the bolitoglossines in the region) from the IUCN Red List of Threatened Species™ (IUCN, 2021). Because range polygons are based on collection points instead of suitability surfaces, they might potentially underestimate or overestimate the geographic range of some of our study species, especially those that lack many collection points (deficient species or species that have restricted distributions). This uncertainty in geographic area estimation can affect some of our metrics.
However, to reduce geographic uncertainty, we used updated IUCN range polygons, which were generated by (1) plotting point data, (2) drawing a minimum convex polygon (MCP) around the points, (3) expanding the range considering the knowledge of habitat preferences, (4) removing areas known to be unsuitable (e.g., unsuitable habitat, elevation limits, or climate/temperature restrictions), and (5) smoothing polygons. Major modifications in range maps (steps 3 and 4) are based on IUCN protocols using environmental data and IUCN expert criteria (Bland et al., 2017;IUCN Red List Technical Working Group, 2018).
Additionally, we generated maps for other 16 bolitoglossines not included in the IUCN dataset, but for which occurrences were available in the literature or in the Global Biodiversity Information Facility (GBIF.org; Table S1), by applying 5-km radius buffers around each occurrence. By combining all maps, we were able to account for 100% of the region's known salamander diversity. To describe relevant indicators of bolitoglossine diversity across space, we quantified five metrics within three dimensions of biodiversity: (a) taxonomic diversity (species richness and endemism), (b) phylogenetic diversity (Faith´s phylogenetic diversity and evolutionary distinctiveness), and (c) functional diversity (body size).

| Taxonomic diversity
We calculated two measures of taxonomic diversity: species richness and endemism. With the full set of range maps, we created a presence-absence matrix at 50-km resolution and estimated species richness by summing all species co-occurring within each grid cell using the R package "letsR" (Vilela & Villalobos, 2015). Given that spatially restricted species may be more prone to extinction risk (Gaston, 1996;Lawler et al., 2003), we also quantified a proxy of endemism by dividing species richness by the mean range sizes of the species occurring within each cell to obtain an estimation of corrected weighted endemism (Crisp et al., 2001). This metric highlights regions of high endemism by assigning greater proportional weight to pools of spatially restricted species. For both metrics, we used all bolitoglossine species with available geographic data (320 species) and performed the analyses using the R package "raster" (Hijmans & van Etten, 2010).

| Phylogenetic diversity
We calculated two measures of phylogenetic diversity: Faith's phylogenetic diversity (PD) and evolutionary distinctiveness (ED).
Conservation policy has found that the phylogenetic dimension of biodiversity is a valuable way to incorporate evolutionary history into decision-making (Faith, 2018;Isaac et al., 2007). Therefore, we reconstructed the phylogeny of bolitoglossines based on the matrix data used in , which includes seven mitochondrial markers (NDI, COI, COII, cyt b, tRNAs, 12S, and 16S) and three nuclear markers (POMC, RAG1, and SLC8A3). We compiled new data available from Genbank (www.ncbi.nlm.nih.gov) for these same markers to enhance the species representation in the phylogenetic reconstruction. With this process, we added 33 species (Table S2) to the Rovito, Vásquez-Almazán, et al. (2015) dataset for a total of 267 species included in our phylogeny, representing 83% of known bolitoglossines.
Using the CIPRES data portal (Miller et al., 2010), we conducted a Bayesian inference analysis in BEAST 1.8.2 (Drummond et al., 2012) with four chains of 50 million generations, sampled every 1000 generations. We configured the input file under the uncorrelated lognormal relaxed clock, with a Yule process tree prior. We assigned the substitution model to each partition according to the results estimated in PartitionFinder v1.0 (Lanfear et al., 2012). We used Tracer v1.7.1 (Rambaut & Drummond, 2007) to assess convergence of the runs, appropriate burn-in, and ensure that effective sample size values were sufficiently high (>200). Finally, a maximum clade credibility tree with median heights was calculated with TreeAnnotator 1.8.2 (Drummond et al., 2012).
With this tree, we estimated two evolutionary metrics that account for phylogenetic diversity: (1) Faith's phylogenetic diversity (PD), which was calculated for the assemblage of all bolitoglossine species occurring in a given grid cell by summing up all their branch lengths (Faith, 1992). Because the branch lengths count the relative number of new features arising in a specific portion of a phylogenetic tree, the PD summarizes the evolutionary history accumulated by the set of species occurring in a given region. In contrast, (2) evolutionary distinctiveness (ED) is a metric that measures species uniqueness as determined by means of phylogeny (Redding & Mooers, 2006). ED scores are assigned to each species by applying a value to each branch equal to its length divided by the number of species composing the branch (Redding et al., 2008). With these species' values, we calculated means considering all species occurring in a grid cell. All phylogenetic analyses were conducted with the R package "ape" (Paradis & Schliep, 2019) over our newly generated phylogenetic hypothesis.

| Functional diversity
Finally, we calculated one value of functional diversity based on body size. Due to its capacity to reflect morphological, physiological, and ecological traits within biological communities, functional diversity provides relevant information to safeguard ecosystem functioning (Petchey & Gaston, 2006). Body size is a commonly used trait in functional ecology due to the availability of body measurements for many species and because it co-varies with many other species features (Gaston & Blackburn, 2008). To determine the variability of this trait, we estimated the standard deviation of adult body size based on the species assemblage occurring in each cell.
For this, we compiled data on body size (measured as maximum snout-vent length-SVL-) for 311 species (>95% of the total diversity of Neotropical salamanders) representing all genera included in this study. Maximum SVL data for each species were obtained from the literature and from direct measurements on museum specimens (Table S3).

| Climatic suitability for Bsal in the Neotropics
We estimated the potential distribution of Bsal in the Neotropics using the Maximum Entropy algorithm (MaxEnt; Phillips et al., 2006).
To build the models, we gathered data of proven presence for Bsal  (Table S4). Further details on the climatic predictors used as well as on the model tuning and evaluation approaches are thoroughly described in Appendix S1.

| Regions and species potentially threatened by Bsal
We used our maps of bolitoglossine biodiversity and the potential distribution of Bsal to define priority regions for conservation. To this end, we overlaid the prediction of Bsal suitability with each of the six regional maps of diversity using ArcGIS v.10.2 (ESRI, 2010).
We then identified the highest-risk areas as those where high Bsal suitability (>0.7) and high bolitoglossine diversity overlapped.
We also used a cross-species approach to define species at risk by Bsal. In this case, we first converted the prediction of Bsal suitability (logistic output) to a binary presence-absence map of potential distribution for the pathogen using the equal training sensitivity and specificity threshold (Liu et al., 2005). Additionally, we used the minimum training threshold for a sensitive analysis. To keep the more conservative approach, we downstream analyses using the equal training sensitivity and specificity threshold. This threshold is considered conservative and suitable for the purpose of delimiting priority conservation areas based on the precautionary principle (Liu et al., 2005). Then, using ArcGIS v.10.2 (ESRI, 2010) we estimated the overlap among Bsal and each species range polygon to quantify the proportion of the range that is predicted to be suitable for the pathogen. We also mapped species-specific information on conservation status according to the IUCN Red List (IUCN, 2021)

| Geographic patterns of bolitoglossine biodiversity
The highest levels of bolitoglossine species richness were associ-

| Climatic suitability for Bsal in the Neotropics
The best-fit species distribution model for Bsal had a parametrization of regularization multiplier of four with a linear-quadraticproduct (LQP) feature class combination (Table S6) occurrences, and a 0% omission rate for the equal training sensitivity and specificity threshold (logistic threshold = 0.271; Table   S6) applied. The MESS analysis also confirmed the similarity of climatic conditions between the current niche of Bsal in Europe and Asia with that of the Neotropics ( Figure S2). The climatic predictor with the highest contribution to the species distribution model was the maximum temperature of warmest month (BIO5), which explained 53.6% of the predicted distribution of Bsal.
The projected model into the Neotropics (Figure 2; Figure

| Regions and species potentially threatened by Bsal
We found a high concordance between high suitability areas for Bsal and high bolitoglossine diversity ( Figure 3). The mountains of Costa Rica (e.g., Cordillera de Talamanca) and Mexico (e.g., eastern  (Figure 3a). Similarly, the pathogen shows high expected suitability in regions of high phylogenetic diversity, including some of the previously cited speciose mountain ranges and other regions, such as the highlands of the northern portion of nuclear Central America (Figure 3b). Regions hosting high functional diversity for neotropical salamanders also are within areas with high suitability for Bsal (Figure 3c).
The potential distribution of Bsal coincides to some degree with the distribution of 88% of bolitoglossine species (Table S8). We found that at least 75% of the geographic ranges of 71 Bolitoglossa species are predicted to be suitable for Bsal (Figure 4a). We found a similar pattern in another 110 species of other genera of bolitoglossine salamanders (Figure 4b). In general, nearly a third of the studied species (90 out 309 excluding 11 which are considered possibly extinct) are categorized as EN or CR on the IUCN Red List, and over 75% of their ranges could be affected by Bsal (see Table S8 for the complete species ranking).

| DISCUSS ION
Our results suggest that the potential introduction of Bsal in the Neotropics could severely affect the taxonomic, phylogenetic, and

F I G U R E 3 Detail of overlap between
Bsal suitability and geographic patterns of the five metrics of bolitoglossine biodiversity in the most diverse regions. . These montane species usually show high levels of isolation and unique environmental adaptations that make them more prone to extinction in the face of current and future threats (Böhm et al., 2016;Cooper et al., 2008;La Sorte & Jetz, 2010).
In accordance with previous analyses at the genus level for salamander diversity across the Neotropics (Wake & Lynch, 1976), phylogenetic diversity is higher in Central Mexico and Costa Rica.
This reflects the heterogeneous evolutionary histories contained in these regions, which are among the most important diversification centers for bolitoglossines . Evolutionary distinctiveness is most evident in Costa Rica and portions of Mexico, but also in northern Guatemala and Honduras. Evolutionary distinct species represent highly isolated portions of the bolitoglossine phylogeny or representatives lacking close relatives (Redding & Mooers, 2006). data and the lack of equilibrium between the organisms and their environment in ongoing invasions (Elith et al., 2010;Katz & Zellmer, 2018;Václavík & Meentemeyer, 2012). To optimize our predictions and improve model transferability (Liu et al., 2020;Yates et al., 2018), we built our SDMs using exclusively proven presences of the pathogen, and given the few and clustered existing occurrence for Bsal, we used the combination of both presences from the native and invaded range to better depict the pathogen's niche. Moreover, we implemented a partition method that improve spatial independence of training and test data for a more robust spatial cross-validation and inspected the similarity between conditions in the calibration and projection regions to assess whether the transferability of our prediction in biologically plausible.
After carefully considering such hurdles implicit in modeling alien species, we found that in the worst scenario, all dimensions of boli-  (Faith, 2018), and functional traits (Petchey & Gaston, 2006) may have dramatic consequences for ecosystem functioning.
It is possible that microclimatic conditions associated with specific microhabitats (e.g., arboreal vs fossorial species) could constrain pathogen transmission, potentially allowing some species to coexist with Bsal. Given that no studies have shown differential vulnerability, we cannot predict which species could be more or less affected by Bsal; therefore, we recommend applying the precautionary principle. Future studies aiming to improve the predictions of Bsal in the Neotropics should also incorporate information of species susceptibility, spatial transmission (Malagon et al., 2020), and mechanistic data (Kearney & Porter, 2009 (IUCN, 2021). Several of these species also have extremely restricted distributions; thus, we suggest prioritizing the development of conservation strategies for these species. Additionally, if the introduction of Bsal occurs, the moretolerant species (e.g., most anurans) could act as reservoirs, thereby favoring the long-term persistence of the pathogen in the environment (Stegen et al., 2017). Therefore, studies to detect and quantify the presence of Bsal in wild amphibians, especially from the diversity hotspots where the effects of Bsal could be the most catastrophic, are urgently needed.
Our Bsal distribution model was developed based on the climate niche of Bsal, assuming equal vulnerability across hosts. Although some species could be more tolerant than others (Martel et al., 2014), by assuming equal vulnerability across hosts we followed the precautionary principle for conservation and decision-making purposes (Persson, 2016). A recent work testing Bsal susceptibility specifically in plethodontid salamanders showed a wide spectrum of responses from species experiencing rapid and high mortality to species with moderate infections (Direnzo et al., 2021). However, for bolitoglossines the evidence remains restricted to two species The potential distribution for Bsal also corresponds with some areas of the predicted distribution of Bd at both the regional level (Ron, 2005)

Recent studies have shown that simultaneous infections with Bd and
Bsal (1) may affect salamander populations more severely than Bsal infection alone (Longo et al., 2019;McDonald et al., 2020) and (2) that pre-existing Bd infections in some species may protect against Bsal (Greener et al., 2020). Given that both scenarios could be possible in the region, our suitability map also establishes a powerful predictive tool, especially in areas with current Bd infections in the Neotropics, as our knowledge on co-infection dynamics improves.
Our study suggests that most bolitoglossine diversity could be at risk if Bsal is introduced into the Neotropics. In this scenario, preventing the introduction of Bsal in the Neotropics is imperative for the