Temporal Changes in Geographic Range and Population Size of the Santa Cruz Black Salamander Free

The status of Lissamphibia as a highly threatened vertebrate taxon is well documented (Wake 1991; Green 2003; Green et al. 2020). The reasons for this are diverse and include factors such as habitat alteration (Apodaca et al. 2012), introduced species (Falaschi et al. 2020), climate change (Blaustein et al. 2010), and disease (Rovito et al. 2009; Caruso and Lips 2013; Cowgill et al. 2021). Some aspects of the natural history of certain species increase their susceptibility to these factors, with restricted range size and low population size being especially pertinent factors correlated with an increased risk of decline and extinction in amphibians (Sodhi et al. 2008). Some range-restricted species also feature low rates of vagility, which can exacerbate their level of risk by decreasing the chances for recolonization of areas where localized population extinctions occur (Welsh and Droege 2001). Small ranges and populations can render amphibians vulnerable to proportionately severe changes in population size such as bottleneck events (Ellstrand and Elam 1993), which lower genetic diversity and reduce the capacity of a species to adapt to changing selective pressures such as environmental changes and emerging diseases (Frankham et al. 2002; Willi et al. 2006). Determining the level and distribution of genetic variation in a range-restricted species is an important step toward drafting an effective conservation management plan.

Knowledge of a species’ ecological niche, or the combination of environmental variables that allows a species to persist, is critical for understanding its basic ecology as well as its current distribution (Grinnell 1924; Wiens and Graham 2005). Ecological niche modeling produces maps of estimated suitability for a species that researchers can use to identify and prioritize geographic areas to protect, or areas to survey where they have yet to find the species (Peterson 2001; Elith and Leathwick 2009). For fossorial species the environmental variable used in the model may not define the niche of the species, but one can assume that they are correlated with the species niche (see Farallo et al. 2020). Nonetheless, the projection of a species’ estimated niche into the past can be useful for interpreting their current range boundaries and for explaining patterns of genetic diversity (Elith and Leathwick 2009; Reilly et al. 2015). Conversely, projections of a species niche into the future can help with conservation planning and management of a species by predicting range expansions, contractions, or shifts, and help focus limited resources to protect habitat that will be suitable for the species as the climate changes (Sinclair et al 2010; Alvarado-Serrano and Knowles 2014). Additionally, the identification of areas that will become increasingly important to the species in the near future can highlight critical dispersal corridors between regions to protect (Ramirez-Villegas et al 2014; Vasudev et al. 2015).

The Santa Cruz black salamander Aneides niger is a lungless, direct-developing species in the family Plethodontidae that is restricted to the Santa Cruz Mountains along the central coast of California (Figures 1A, 1B, and 2A). It is terrestrial, saxicolous, and fossorial and occurs in a number of distinct habitats including redwood Sequoia sempervirens forest, chaparral, and mixed evergreen forest, though researchers most commonly find it in rocky, streamside microhabitats (Figures 1C and 1D). The range size for this species (∼1,700 km²; Google Earth 2022) is considered small for a vertebrate. It is restricted enough that it was assigned the highest-ranking criteria for range size (10/10) in the species’ assessment as a California species of special concern (SSC) for which it was ranked a SSC Priority 3, defined as “…clearly at risk but likely are not experiencing a substantial and immediate threat of expiration, although the potential for this threat to develop exists if no management actions are undertaken.” (Thompson et al. 2016). Additionally, the species is listed as endangered (EN; code B1ab(iii)) by the International Union for the Conservation of Nature based on 1) the extent of occurrence, 2) its occurrence in threat-defined locations, and 3) the continuing decline in the extent and quality of its habitat (IUCN 2022). Wildlife biologists find this species at low densities across its range and encounter it less frequently than other sympatric Plethodontid salamanders such as the California slender salamander Batrachoseps attenuatus, the Yellow-eyed ensatina Ensatina eschscholtzii xanthoptica, and the congeneric arboreal salamander Aneides lugubris (B.R.K. and S.B.R. personal observations). The species (a subspecies at the time) was so difficult to find in the late 20th century that R. C. Stebbins indicated in his field guide that it had suffered from severe population size and range decline (Stebbins 2003). Although demonstrating a decline in cryptic taxa is challenging, historic accounts and museum collections indicate that the species was not always so difficult to find, further suggesting that a very recent population decline may have occurred (D. Wake [Museum of Vertebrate Zoology], personal communication; Van Denburgh 1895; Myers and Maslin 1948; Thomson et al. 2016).

There is surprisingly little information on the basic biology of the Santa Cruz black salamander, despite this being critical for developing a management plan and carrying out conservation actions to preserve the species. There are no population-level surveys from which to estimate population sizes or trends, which was the main reason that Thompson et al. (2016) did not give the species a higher SSC priority ranking. Even basic genetic work is mostly lacking except for phylogenetic studies showing its distinctiveness from the rest of the black salamander species complex (Rissler and Apodaca 2007; Reilly and Wake 2015; Reilly and Wake 2019). Rissler and Apodaca (2007) created a present-day species distribution model for the Santa Cruz black salamander, but did not project this model to the past or future. Another study projected future climatic suitability for the northern members of the black salamander complex, but not for the Santa Cruz black salamander (Early and Sax 2011).

In this study we utilize sequence data generated by previous studies to perform a preliminary genetic assessment to address three key factors: 1) genetic variation within the species, 2) phylogeographic structure, and 3) genetic signatures of population size change. We then estimate species distribution models for multiple points in the past to look for evidence of range restrictions relative to the current range size, which likely coincide with population size reductions. We then projected species distribution models into the future (2081–2100) based on four climate models and two emissions scenarios with the goal of informing conservation planning and management of the species. We hope that this study, while not comprehensive enough to fully inform a conservation management plan, will serve as a starting point to guide future studies and surveys of the Santa Cruz black salamander.

We gathered sequence data (Table S1, Supplemental Material) for Santa Cruz black salamander samples from Reilly and Wake (2015) and Rissler and Apodaca (2007) including up to three mitochondrial genes (ND4, cytb, 12S) for 14 samples from six localities (Data S1, Supplemental Material), and one nuclear intron (POMC) plus 10 anonymous nuclear loci for seven samples from four localities (Data S2, Supplemental Material). Table S2 (Supplemental Material) contains GenBank accession numbers. We inferred haplotype sequences from the nuclear loci using the software program PHASE v2.1 (Stephens et al. 2001; Stephens and Scheet 2005) using the default threshold of 0.9. We used DNASP v5 (Librado and Rozas 2009) to estimate nucleotide diversity and to analyze each locus with four similar methods to detect signals of directional or balancing selection as well as population expansion or contraction including Tajima’s D, Fu’s F_S, Fu and Li’s D*, and Fu and Li’s F*. We plotted the mismatch distribution for all mitochondrial DNA (mtDNA) sequences comparisons in R 4.2.2 (R Core Team 2022) as a further test of population change.

We imported the concatenated mitochondrial alignment of 1,732 bp (Data S1, Supplemental Material) into IQTREE v2.1.3 where maximum likelihood phylogeny estimation was coupled with automatic model selection (Kalyaanamoorthy et al. 2017; Minh et al. 2020). We partitioned the data by gene and the models of sequence evolution chosen by Bayesian information criterion for each partition were 1) the Tamura and Nei (1993) model with empirical base frequencies and a discrete gamma model of rate heterogeneity with four rate categories (TN + F + G4) for ND4, 2) the Hasegawa-Kishino-Yano (HKY, Hasegawa et al. 1985) model with empirical base frequencies and a discrete gamma model of rate heterogeneity with four rate categories (HKY + F + G4) for cytb, and 3) a transition model where AC = CG and CG = GT with unequal, empirical base frequencies and a discrete gamma model of rate heterogeneity with four rate categories (TIM3 + F + G4) for 12S. We assessed node support by 1,000 ultrafast bootstrap replicates (UFBoot; Hoang et al. 2018) and 1,000 SH-aLRT branch tests, where UFBoot ≥95% and SH-aLRT ≥80% indicate robust node support.

We undertook a time-calibrated Bayesian approach to obtain a rough timeframe of diversification events using BEAST v2.5 (Bouckaert et al. 2019) on the mitochondrial dataset (Data S1, Supplemental Material). The HKY model of sequence evolution was the highest ranking among the various models implemented by BEAST. As such, we conducted two independent runs under the HKY model, the Yule tree model prior, and a strict molecular clock with a rate of 0.01, which approximates the average mitochondrial mutation rate across all vertebrates (Allio et al. 2017). We ran the program for 50 million generations, sampling trees every 1,000 generations. The XML file used to run this analysis is in Data S3 (Supplemental Material). After checking both runs for stationarity with TRACER v1.7 (Rambaut et al. 2018) by visually inspecting individual parameter trace plots and ensuring all effective sample size values were >200, we removed the first 10% of saved trees from each run and combined them to produce a maximum clade credibility tree.

We tested the mitochondrial dataset for signatures of isolation-by-distance by conducting Mantel tests on matrices of Euclidean geographic distances and Edwards’ genetic distances with the R package ADEGENET (Jombart 2008). We assessed statistical significance in the form of a P value by employing 100,000 permutations. We simplified the mitochondrial dataset to only variable sites (Data S4, Supplemental Material), and used GPS coordinates (Data S5, Supplemental Material) for each sample to compute geographic distances. For each sample pair we plotted the genetic distance against geographic distance and local point density, and visualized it as a colored cloud over the plot, which we calculated using two-dimensional kernel density estimation in the R package MASS (Ripley and Venables 2002). Data S6, Supplemental Material, contains the R code used for this analysis.

We inferred temporal population size changes from the full dataset (3 mtDNA loci + 11 nuclear DNA loci) using the extended Bayesian skyline plot (EBSP) method implemented in BEAST v2.5 (Heled and Drummond 2008; Bouckaert et al. 2019). We analyzed the mitochondrial data (Data S1, Supplemental Material) under the HKY model of sequence evolution and analyzed the nuclear genes (Data S, Supplemental Material 2) under the Jukes and Cantor (1969), or JC69, model. Under a strict clock model, the mitochondrial gene served as the reference due to its higher rate of mutation, while we inferred the mutation rates of the nuclear genes relative to mtDNA. The mitochondrial mutation rate was set to 0.01, which is the approximate average rate across vertebrates (Allio et al. 2017). We set a scale factor of 0.5 for the mtDNA and a factor of 2.0 for nuclear genes. We set the hyperprior on the mean of the population size distribution to a tight Gaussian prior centered on 1.0 with a standard deviation of 0.1. After multiple test runs checking for consistency across runs, we ran the program for 100 million generations, sampling every 5,000 generations (the XML file is Data S7, Supplemental Material). After checking for convergence by confirming that effective sample size values were >200 in TRACER v1.7 (Rambaut et al. 2018), we viewed the individual population function lines as well as the median and 95% central posterior density bounds using the plotEBSP command in R (Data S8, Supplemental Material).

Given the topography of the Santa Cruz Mountains and the proximity of locations for nuclear data (localities 1, 2, 4, and 6 from Figure 1A) we treated locality 1 + 2 as our “west” population and locality 4 + 6 as our “east” population for demographic parameter estimation with the software G-PHOCS (Gronau et al. 2011). We do not have nuclear data from localities 3 and 5 so did not consider those localities for this analysis. We analyzed only the nuclear data (Data S2, Supplemental Material) for this coalescent-based program that assumes the data derives from unlinked, neutrally evolving loci. We estimated the contemporary and ancestral effective population sizes (θ), population divergence time (τ), and migration rates (m). We set the priors for θ and τ to 1 for alpha and 10,000 for beta and migration priors were set to 0.002 and 0.00001 for alpha and beta, respectively. The control file used to set parameter values is in Data S9, Supplemental Material. We ran the program for 1,000,000 generations and confirmed parameter convergence by checking that effective sample size values were >200 in TRACER v1.7 (Rambaut et al. 2018). We removed the first 10% of saved estimates as burn-in and used the rest to estimate the mean and 95% posterior density for each parameter. For population size estimates we are looking for major changes between the ancestral population and the sum of the two extant populations. To convert parameters, we used a nuclear mutation rate (μ) of 1 × 10⁻⁹ mutations · site⁻¹ · y⁻¹ (Allio et al. 2017) and a 2-y generation time to get some rough time scale of the time over which population size changes may have occurred. Wildlife biologists do not know the generation time for this species and it likely has more to do with size rather than age, and as such we stress that our converted parameter estimates should be interpreted cautiously. For the purposes of this study, we interpreted migration rates as being absent, equal, asymmetric, or unidirectional rather than actual numbers of migrants per generation.

We used the historical (1970–2000) bioclimactic data from WorldClim version 2.1 (Fick and Hijmans 2017) to generate a species distribution model, which we used for making projections to future climates. We used data from four different climate models hosted on the WorldClim website (https://www.worldclim.org/data/worldclim21.html) to gauge how future climate change could impact the Santa Cruz black salamander. Specifically, we used future climate projections from the 2081–2100 time frame from the CMCC ESM2 (Lovato et al. 2022), MRI ESM2 (Yukimoto et al. 2019), EC Earth3 Veg (Döscher et al. 2022), and MIROC 6 (Tatebe et al. 2019) climate models. For each of the models, we used data from the lowest (shared socioeconomic pathway 126) and highest (shared socioeconomic pathway 585) emissions projections to assess climate suitability for the best-case and worst-case scenarios. We also used the historical (1960–1990) bioclimatic data from WorldClim version 1.4 (Hijmans et al. 2005) to generate a species distribution model that was projected to paleo climates, because the WorldClim paleo climate data were calibrated relative to the WorldClim 1.4 data (https://www.worldclim.org/data/v1.4/paleo1.4.html). We used data for the last glacial maximum (∼22,000 years ago) from the CCSM4 model and data for the last interglacial period (120,000–140,000 years before present) from Otto-Bliesner et al. (2006). We used climate data at a resolution of 2.5 min, except for the data for the last-interglacial period, which was only available at 30 arc-s resolution. We cropped the raster data to a region encompassing the San Francisco Bay Area based on the western extent of Sonoma County, the northern extent of Napa County, and the southern and eastern extent of Monterey County (minimum latitude = 35.78898, maximum latitude = 38.86424, minimum longitude = −123.5339, maximum longitude = −120.214). We repeated this methodology with an expanded region to test the robustness of our results by expanding the northern extent to the top of Mendocino County, and the southern and eastern extent to the southern and eastern boundaries of San Luis Obispo County. We have displayed these over the California shape file to put the data in geographic context (Data S10, Supplemental Material).

We used a subset of the 19 bioclimatic variables to avoid having our species distribution models be overly restrictive due to overfitting. We used the following four variables in our models, which were not highly correlated with one another (all correlations less than or equal to 0.56): maximum temperature of the warmest month (BIO5), mean temperature of the coldest quarter (BIO11), precipitation seasonality (coefficient of variation; BIO15), and precipitation of the wettest quarter (BIO16). These variables are relevant to the biology of salamanders in the Bay Area, which are dependent upon the moisture levels of their environment and are active during the fall and winter when there is rainfall.

We modified R code from Smith (2022) to generate our species-distribution models (Data S11, Supplemental Material), which used the R-packages raster (Hijmans 2022), rgdal (Bivand and Keitt 2022), and rgeos (Bivand and Rundell 2021) for manipulating GIS data, the legendary package (Smith 2021) for examining correlations among the climate variables, and the dismo (Hijmans et al. 2021) and maxent (Phillips 2021) packages to generate the species distribution models. We used 69 unique locality points to establish sites where Santa Cruz black salamanders were present based on museum specimen localities downloaded from Arctos museum database (https://arctos.database.museum), with manual removal of locality points deemed to be erroneous (Data S12, Supplemental Material). We randomly sampled the climate at 1,000 background points to represent the climate of the Bay Area. We then used maxent trained with these background points and the presence points for the Santa Cruz black salamander to generate a maxent species distribution model (Phillips et al. 2006). We then projected the species distribution model to paleo climates (our WorldClim 1.4 model) and to future climates (our WorldClim 2.1 model).

The mitochondrial alignment of 1,732 bp contained 32 variable sites with nucleotide diversity of 0.0054 (Table 1). When considering west (localities 1–3) and east (localities 4–6) populations the nucleotide diversity was approximately twice as high in the east (0.0049) as in the west (0.0024). Four of the nuclear loci (POMC, SR4, SR8, SR17) contained no variable sites and the other nuclear loci had nucleotide diversity values ranging from 0.0163 to 0.0004. There do not appear to be any significant signals of directional or balancing selection, or major population growth/contraction as only one neutrality test, Fu and Li’s D*, returned a significant value (P < 0.02) for locus SR2. The mismatch distribution shows a skewed unimodal distribution with most values close to zero, which is indicative of a recent population expansion (Figure S1, Supplemental Material).

The maximum likelihood tree displays nine unique mitochondrial haplotypes within A. niger (Figure 2B). The tree is rooted by a haplotype from locality 6 (SH-aLRT = 97.8, UFBoot = 100), with two sister groups (SH-aLRT = 84.1, UFBoot = 83), one containing haplotypes from localities 1 and 4 and the other containing haplotypes from localities 2, 3, and 5 (Figure 2B). The time-calibrated Bayesian phylogeny estimates the mean age of the most recent common ancestor at ∼0.66 million years ago (Ma) (Figure 2C). There is a significant signature of isolation-by-distance within the mitochondrial data (simulated P value = 0.00059; Figure S2, Supplemental Material; Figure 2D).

The EBSP estimates a generally stable population size from ∼2.5 Ma until 1 Ma, after which a population decline is inferred with a low point centered around 0.2 Ma (Figure 3A; Figure S3, Supplemental Material). Population growth is detected starting around 0.2 Ma until present time, though the individual trace lines and lower 95% central posterior density show a slight signal of population decline near the present day. The demographic parameters estimated from the coalescent-based G-PHOCS program suggest that the panmictic ancestral population existed until thousands or tens of thousands of years ago (Table S3). These estimates are considerably more recent than those generated by our EBSP. We estimate the individual western and eastern populations to have effective population sizes that are roughly equal to each other and are each less than 10% of the ancestral effective population size (Figure 3B, Supplemental Material). Migration estimates are generally low (rates <0.5 in either direction are considered low; see Shaffer and Thompson 2007) and roughly equal in each direction with 2Nm values ∼0.26 (Table S3, Supplemental Material).

The current distribution models using both the WorldClim 1.4 and 2.1 datasets produced comparable climactic ranges that both highlight the central and southern portions of the Santa Cruz Mountains as highly suitable (Figures 4A and 4B). Other suitable areas not inhabited by Santa Cruz black salamanders include the western portions of Alameda and Contra Costa counties, and a patch in south-central Sonoma County. The projection to the last glacial maximum time period (∼22,000 y ago) estimates that the most suitable habitat would have been along the lower-elevation foothills of the southeastern portion of the Santa Cruz Mountains (Figure 4C). However, the overall climatic suitability was substantially reduced during the last glacial maximum compared to the present (note the big difference in color scale of suitability between panels of Figures 4A and 4C). A projection to the last interglacial time period (∼120,000–140,000 y ago) estimates that the most suitable habitat would have been more restricted to the southwestern part of the Santa Cruz Mountains, primarily within Santa Cruz County (Figure 4D), but that climatic suitability was much more similar to the present than the last glacial maximum (note color scale).

Future projections of climatic suitability based on the lowest emissions scenario (shared socioeconomic pathway 126) vary between the four climatic models we assessed. Two models (CMCC-ESM2, EC-Earth3-Veg) forecast an overall restriction of habitat, with reductions in habitat suitability occurring in the eastern and middle parts of the species range within the Santa Cruz Mountains (Figure 4E; Figure S4A, Supplemental Material). In contrast, two models (MRI-ESM2-0, MIROC6) forecast that the Santa Cruz Mountains would continue to have high suitability and that habitat suitability would increase in coastal Monterey County (Figure 4G; Figure S4C, Supplemental Material). Future projections based on the highest emissions scenario (shared socioeconomic pathway 585) are all similar across the four climate models and show a drastic restriction of habitat, with northern coastal areas of the species range maintaining the highest suitability (Figures 4F and 4G; Figures S4B and S4D, Supplemental Material). The duplicate set of analyses with an expanded geographic region showed very similar results for each model (Figure S5, Supplemental Material).

This study found signatures of population contraction and growth since the late Pleistocene when examining all of the available genetic data (Figure 3). We would not have detected this finding within the individual loci on their own as they did not show significant signatures of population change in the locus-specific tests that we conducted (Table 1). Reconstructions of the species distribution at the last glacial maximum and the last interglacial both showed a much smaller area of available habitat than for the present day, with habitat suitability being particularly restricted during the last glacial maximum (Figure 4). This suggests that the range size of the Santa Cruz black salamander, and presumably the population size, would have been smaller than at present, particularly during glacial periods, which regularly occurred for over 2 million years (Hewitt 2000). Previous studies also detected genetic evidence of range expansion during this time period within the closely related speckled black salamander Aneides flavipunctatus in central Mendocino (Reilly et al. 2012) and Klamath black salamander Aneides klamathensis in the Klamath Mountains (Reilly et al. 2013). This restricted range and population scenario would be consistent with the dramatic dip in population size detected with our EBSP around the timeframe of 100,000–300,000 y ago. Habitat patches for Santa Cruz black salamanders in the past were located along the coast during the last interglacial and restricted to the inland eastern-margin foothills during the last glacial maximum. The mitochondrial data are consistent with the projected restriction of the species range during the last glacial maximum, because the eastern populations contain more deeply divergent haplotypes and higher nucleotide diversity than western populations, which is consistent with the eastern foothills being refugia for the Santa Cruz black salamander. The projections to paleo climate thus suggest that populations either slowly migrated with the shifting habitat suitability, or that population size fluctuated greatly with population growth in one region of the Santa Cruz Mountains mirrored by population decline in the other regions of the range.

The Santa Cruz black salamander features low genetic diversity (all reported π are for mtDNA; π = 0.005) and a general lack of genetic differentiation among populations in comparison to other regional salamanders with similar range sizes such as the yellow-eyed ensatina Ensatina eschscholtzii xanthoptica Bay Area clade (π = 0.029; Kuchta et al. 2009), the garden slender salamander Batrachoseps major northern clade (π = 0.052; Martínez-Solano et al. 2012), and the Mount Lyell salamander Hydromantes platycephalus northern clade (π = 0.033; Rovito 2010). However, the nucleotide diversity in A. niger is slightly greater than some lineages such as the California Giant Salamander Dicamptodon ensatus southern clade (π = 0.003; Callahan 2010) and the red-bellied newt Taricha rivularis (π = 0.0004; Reilly et al. 2014). This likely places the Santa Cruz black salamander in a vulnerable position in terms of its capacity to adapt to a rapidly changing environment (Willi et al. 2006). This lack of diversity is especially concerning when considering the rapid climatic change forecast for the near future.

We estimated the future distributions of the Santa Cruz black salamander under both the lowest and highest emissions scenarios across four different climate models to gauge the range of possible outcomes under global climate change (Figure 4; Figure S4, Supplemental Material). Two climate models under the lowest emissions scenario and all models under the highest emissions scenario forecast reductions in habitat suitability. Conservation management plans for the Santa Cruz black salamander should therefore account for likely reductions in range size in the future, especially because the low emissions scenario is based upon optimistic projections of socioeconomic change. The higher emissions scenarios suggest that habitat along the western margin of the Santa Cruz Mountains will maintain suitability, so this region is a high priority for conservation to protect against a worst-case climate scenario. All of the high-emissions forecasts suggest that the more northern portion of the range and San Francisco Peninsula will become increasingly suitable. These areas are currently unoccupied by the species but represent the only reasonable area for them to expand their range into, and thus represent an important asset for persisting through future climate changes. However, macroclimatic data may not fully encompass the evolutionary dynamics of a fossorial species’ ecological niche. Wildlife managers should consider the propensity of salamanders to rely heavily on microclimatic variables when interpreting our results (Farallo et al. 2020).

Our genetic data can be used to infer long-term population size changes in the past (Figure 3), but they cannot reliably assess the anecdotal evidence of a historically recent (<100-y) population decline in the Santa Cruz black salamander unless it was extreme. If there has been a recent population decline it may be difficult to detect without a large genomic dataset, especially if there is a small sample size, because there is a time lag between the point at which a bottleneck occurs and when its signature shows up in the genes of the population (Nei et al. 1975; Peery et al. 2012). Thus, assessing the possibility of a further population decline in the Santa Cruz black salamander requires repeated population surveys going forward. One of the difficult aspects of proving that a secretive species has declined is that there are usually not sufficient survey efforts in the past to provide a baseline. This is the case for A. niger. Widespread declines of eastern North American salamanders in the genus Plethodon were only brought to light by the extensive long-term survey efforts conducted at hundreds of sites over a nearly 50-y period (Highton 2005).

While our genetic and range suitability data do not address a possible recent population decline in A. niger, certain recent activities may be implicated in a population decline. In the late 1800s clear-cut logging was extensive and the Santa Cruz Mountains were the epicenter of redwood harvesting in California (Evarts and Popper 2001). Logging changed the composition and proportions of tree species and heavily altered the fire ecology of the region. Slash fires killed entire sections of forest, sterilized the soil, and converted many forested regions into chaparral (Greenlee and Langenheim 1990). These actions would have dramatically altered the temperature and humidity of black salamander microhabitats across their range. Another potential culprit is the arrival of chytrid fungus to the Santa Cruz Mountains in the 20th century. Swabbing of museum specimens from the Santa Cruz Mountains detected chytrid fungus in California newts Taricha torosa starting in the 1940s (Chaukulkar et al. 2018), in California slender salamanders in the 1960s (Sette et al. 2015), and in arboreal salamanders in the 1940s (Cowgill et al. 2021). Chytrid swabs of Santa Cruz black salamander specimens from these time periods should be taken and examined.

This study has revealed that the Santa Cruz black salamander has low genetic variation in both mitochondrial and nuclear markers and limited phylogeographic structure. The genetic data showed a reduction in population size in the late Pleistocene, with a recovery towards the present. Consistent with this, projection of species distribution models into past climates suggests that the habitat area was smaller than today, particularly during the last glacial maximum. Projections of the range into the future highlight the need for protection of land in the western and northern parts of the Santa Cruz Mountains, along with the establishment of dispersal corridors connecting those areas to the currently occupied portions of the range. The Santa Cruz black salamander is in critical need of baseline population surveys as well as a genomic-scale conservation genetic study to inform a management plan that will ensure long-term persistence of this charismatic species.

Please note: The Journal of Fish and Wildlife Management is not responsible for the content or functionality of any supplemental material. Queries should be directed to the corresponding author for the article.

Data S1. Concatenated mitochondrial DNA sequence alignment of 1,749 base pairs for 26 Aneides samples including 14 Santa Cruz black salamanders Aneides niger in NEXUS file format. Samples were collected from Northern California and Southern Oregon between 1979 and 2012. Bases 1–734 are the ND4 gene, bases 735–1,172 are the cytb gene, and bases 1,173–1,749 are the 12S gene.

Available: https://doi.org/10.3996/JFWM-23-004.S1 (46 KB NEX)

Data S2. DNA sequence alignments for the 11 nuclear loci used in the analysis of the Santa Cruz black salamander Aneides niger collected from the Santa Cruz Mountains, California between 1979 and 2012. The format of this .aln file includes the number of loci on the first line followed by each individual gene alignment. Each gene alignment’s first line contains the gene name, the number of sequences included, and the gene length in base pairs.

Available: https://doi.org/10.3996/JFWM-23-004.S2 (32 KB ALN)

Data S3. Input file for BEAST phylogenetic analysis of the mitochondrial data for samples collected from the Santa Cruz Mountains, California between 1979 and 2012. This XML format input file was created using the software program BEAUti and contains all of the sequence data and parameter values needed to run the BEAST software program.

Available: https://doi.org/10.3996/JFWM-23-004.S3 (32 KB XML)

Data S4. Variable sites from the mitochondrial DNA alignment used for isolation-by-distance (IBD) analysis. Samples were collected from the Santa Cruz Mountains, California between 1979 and 2012. These data are in STRUCTURE (.stru) format where base calls for biallelic SNPs are converted to the numbers 1 or 2 based on which allele is found. The IndID refers to the sample number corresponding to the Data S5 (Supplemental Material) file containing the sample GPS coordinates.

Available: https://doi.org/10.3996/JFWM-23-004.S4 (3 KB STRU)

Data S5. Geospatial coordinates for Santa Cruz black salamander Aneides niger samples collected from the Santa Cruz Mountains, California between 1931 and 2012. Mitochondrial data used in the isolation-by-distance (IBD) analysis is in .csv file format. The header label “Field No.” corresponds to the sample number in the sequence data file Data S4 (Supplemental Material), the X label refers to the longitude and the Y label refers to the latitude in WGS84 coordinates.

Available: https://doi.org/10.3996/JFWM-23-004.S5 (1 KB CSV)

Data S6. The R commands used to generate the isolation-by-distance plot and Mantel test results for the Santa Cruz black salamander Aneides niger. These R commands were written in 2023. Input file names and the file locations in the R code should be changed to match the location and titles of files Data S4 and Data S5 (Supplemental Material).

Available: https://doi.org/10.3996/JFWM-23-004.S6 (1 KB R)

Data S7. Input file for the Bayesian skyline plot analysis (EBSP) of the Santa Cruz black salamander Aneides niger. Samples were collected from the Santa Cruz Mountains, California between 1979 and 2012. The file, created by the software BEAUti, contains the DNA sequence data and all parameter settings for the analysis, which is carried out within the software program BEAST.

Available: https://doi.org/10.3996/JFWM-23-004.S7 (155 KB XML)

Data S8. The R commands to plot the Bayesian skyline plot analysis (EBSP) results as seen in Figure 3a. These R commands were written in 2023.

Available: https://doi.org/10.3996/JFWM-23-004.S8 (6 KB R)

Data S9. The control file used to run the G-PhoCS analysis along with the nuclear data input file (Data S2, Supplemental Material). This file was prepared in 2023. The file defines the parameter settings and tells the program which samples belong to which populations. The text for the first parameter in the file, “seq-file”, should match the file name of the sequence file (Data S2, Supplemental Material).

Available: https://doi.org/10.3996/JFWM-23-004.S9 (2 KB CTL)

Data S10. Shape files for California geospatial layers downloaded in 2023 including state and county boundaries that were used in the species distribution modeling analysis.

Available: https://doi.org/10.3996/JFWM-23-004.S10 (505 KB DBF, PRJ, SHP, XML, SHX, QPJ)

Data S11. R code used to generate the species distribution models for the Santa Cruz black salamander Aneides niger. The R code was written in 2023 and the distribution data includes museum samples collected between 1931 until 2012. This code is to be used in conjunction with the shape files in Data S10 (Supplemental Material), the localities in Data S12 (Supplemental Material), and the climate data which are downloaded from WorldClim (the URLs for these data are in the main manuscript text and in the R code).

Available: https://doi.org/10.3996/JFWM-23-004.S11 (25 KB R)

Data S12. Santa Cruz black salamander Aneides niger localities from the Santa Cruz Mountains, California, collected between 1931-2012 used to generate the species distribution models. The first line includes data column labels including museum number, longitude, and latitude. The longitude and latitude values are in the WGS84 coordinate format.

Available: https://doi.org/10.3996/JFWM-23-004.S12 (3 KB CSV)

Table S1. Sampling information for all Santa Cruz black salamanders Aneides niger used in genetic analysis including museum number and collecting locality. All samples are from the Santa Cruz Mountains, California, collected between 1979-2012. Locality # refers to Figure 1A in the main manuscript.

Available: https://doi.org/10.3996/JFWM-23-004.S13 (10 KB XLSX)

Table S2. GenBank accession numbers for Santa Cruz black salamander Aneides niger sequence loci analyzed in this study. All samples were collected from the Santa Cruz Mountains, California, collected between 1979-2012. These accession numbers can be input into the GenBank website (https://www.ncbi.nlm.nih.gov/genbank/) to access the sequences along with their associated data.

Available: https://doi.org/10.3996/JFWM-23-004.S14 (11 KB XLSX)

Table S3. Converted demographic parameters for Santa Cruz black salamanders Aneides niger from G-PhoCS analysis assuming a 1 × 10⁻⁹ mutations · site⁻¹ · y⁻¹ mutation rate and a 2-y generation time. All samples were collected from the Santa Cruz Mountains, California, between 1979-2012. Formulas to convert parameter estimates to real-world values are in the supplemental materials of Gronau et al. (2011). The 95% HDP columns refer to the 95% posterior density values at the low and high end of the distribution, and give a sense of how robust the mean values are.

Available: https://doi.org/10.3996/JFWM-23-004.S15 (10 KB XLSX)

Figure S1. Mismatch distribution among the concatenated mitochondrial sequences of the Santa Cruz black salamander Aneides niger samples collected from the Santa Cruz Mountains, California, between 1979-2012. The number of nucleotide site differences was calculated for each pairwise comparison of sequences. The blue line is the mean density curve fitted to the histogram. This shape of density curve, called a skewed unimodal distribution, is consistent with a recent population expansion demographic scenario.

Available: https://doi.org/10.3996/JFWM-23-004.S16 (171 KB PDF)

Figure S2. Mantel test for isolation-by-distance (IBD) from a matrix of Euclidean geographic distances and a matrix of Edwards’ genetic distances (mtDNA) of Santa Cruz black salamanders Aneides niger. All samples were collected from the Santa Cruz Mountains, California, between 1979-2012. Histograms along the x-axis are permuted values of the correlation between matrices when spatial genetic structure is absent, and the black diamond indicates the original value of the correlation between matrices. Significant spatial structure (IBD) results in the original value occurring outside of the range of the reference value, which is seen here (P = 0.00059).

Available: https://doi.org/10.3996/JFWM-23-004.S17 (5 KB PDF)

Figure S3. Histogram of 95% of gene tree-branching events through time (millions of years) from the extended Bayesian skyline plot (EBSP) analysis for Santa Cruz black salamanders Aneides niger. All samples were collected from the Santa Cruz Mountains, California, collected between 1979-2012. Note that most branching events have occurred relatively recently (<200,000 y), and that almost no branching events occurred beyond 0.8 Ma.

Available: https://doi.org/10.3996/JFWM-23-004.S18 (6 KB PDF)

Figure S4. Future range (years 2081–2100) suitability for Santa Cruz black salamanders Aneides niger under alternative climate models and greenhouse gas emission scenarios. The localities of all samples are from the Santa Cruz Mountains, California, collected between 1931-2012. Panels A and B were estimated with the EC-Earth3-Veg model and panels C and D were estimated with the MIROC6 model. Panels A and C are assuming a low greenhouse gas emissions scenario (shared socioeconomic pathway 126 [ssp 126]) while panels B and D are assuming a high greenhouse gas emissions scenario (shared socioeconomic pathway 585 [ssp585]).

Available: https://doi.org/10.3996/JFWM-23-004.S19 (1,520 KB PDF)

Figure S5. Santa Cruz black salamanders Aneides niger range suitability models for the current time period (A, B), the last glacial maximum (C), the last interglacial period (D), and in the future (E–L). This analysis is based on niche models estimated using samples collected from the Santa Cruz Mountains, California, between 1931-2012. These are the same models as in Figure 4 and Figure S3 (Supplemental Material) except the geographic bounds have been extended to the northern Mendocino County, and to the southern and eastern San Luis Obispo County boundary.

Available: https://doi.org/10.3996/JFWM-23-004.S20 (4,634 KB PDF)

We thank Mitch Mulks, Jon Hirt, and Alison Davis-Rabosky for help with field collection of tissues. We thank the late Barry Sinervo, the late David Wake, Sharyn Marks, and W. Bryan Jennings for their support of this research. S.B.R. thanks Save-the-Redwoods League for funding previous studies that resulted in the genetic dataset analyzed here. We also thank the editors and two anonymous reviewers for their comments, which helped to improve and refine the manuscript.

Any use of trade, product, website, or firm names in this publication is for descriptive purposes only and does not imply endorsement by the U.S. Government.