Population isolation in the Plains spadefoot toad: causes and conservation implications

Study system

Spea bombifrons is a fossorial amphibian species with a range that extends throughout much of western North America, from southern Canada into northern Mexico (Bragg, 1965; Stebbins, 1985; Elliott, Gerhardt & Davidson, 2009). Although some range maps show a continuous distribution (e.g., Elliott, Gerhardt & Davidson, 2009; IUCN SSC Amphibian Specialist Group, 2015), other sources identify a completely isolated population in southern Texas that is separated from the rest of the range by over 400 km (Stebbins, 1985; Dodd, 2013) (Fig. 1).

In arid regions, S. bombifrons spend most of the year estivating underground, emerging only during monsoon rains in the summer to feed and reproduce (Bragg, 1965). The larval period is very short and can last as little as 2 weeks (Dodd, 2013). Little is known about movement and dispersal in this species, although dispersal is likely limited by their short active period and reliance on ephemeral ponds for moisture and reproduction. Previous population genetic work has shown that S. bombifrons has the highest genetic diversity in the central portion of its range, with lower diversity at the range margins, although southern Texas was excluded from this study (Rice & Pfennig, 2008). This pattern of decreasing genetic diversity at range margins is relatively common (reviewed by Eckert, Samis & Lougheed, 2008) and can result from recent range expansion from a central ancestral population towards the range edges as seen in the spadefoot toads (Rice & Pfennig, 2008) or from decreased population sizes at range margins (Gaston, 2003).

Geographic distribution

To map the known distribution of this species, we used historic museum collection records. Data were downloaded from VertNet (https://vertnet.org/, downloaded on 16 March 2018) or provided directly by museum curators. Records that included collection locality descriptions but were missing latitude and longitude data were assigned initial latitudes and longitudes using the web tool GeoLocate (http://www.geo-locate.org/web/default.html, Rios & Bart, 2010). All records were then manually verified and the precision of the coordinates was estimated to the nearest kilometer following the best practices guidelines established in Chapman & Wieczorek (2006). All records with an accuracy measure of greater than 10 km were discarded. This distance excluded records that were too vague to be useful in identifying the collection locality (e.g., records listing only a state or other large geographic region) while maintaining enough records with reasonably precise locality data (e.g., the level of a specific town or city) to ensure good geographic coverage across the species’ entire range. Additionally, we downloaded research grade observations from iNaturalist for Spea bombifrons (downloaded though GBIF.org (2024) on 07 June 2024) to further demarcate the observed distribution of the species.

Each locality was then mapped in ArcGIS 10.8.2 (ESRI, 2020). The Central Range of S. bombifrons was demarcated by a convex hull polygon added around all locality records outside of South Texas with a 20 km buffer added to account for the fact that the species likely occurred beyond the recorded localities. The South Texas range was similarly demarcated with a convex hull polygon and a 20 km buffer.

As a comparison group to test for the effects of sampling effort and to establish whether other species showed range disjunctions that could represent a geographic barrier to gene flow, we identified 27 additional frog and toad species that occurred in similar habitats to S. bombifrons in Texas. Locality records from those species were downloaded from VertNet (between March and April of 2018) and overlaid on the spadefoot toad distributional map. Of these species, those that did not have any records within at least 20 km of any S. bombifrons locality (indicating potential differences in habitat use) and those species with less than 10 unique localities (indicating rarity) were discarded. Twenty species remained. To simplify the visualization of these species’ ranges, records that fell more than 150 km outside the IUCN defined range of S. bombifrons were removed. Each of these species was then coded for their occurrence within each of the three regions of our study: 1) the Central Range of S. bombifrons; 2) the South Texas range; and 3) the area of hypothesized disjunction between these two regions.

Species distribution modeling

Population genetic analyses

We used a complementary S. bombifrons population genetics dataset to determine whether the population genetic signatures are consistent with the South Texas population being genetically distinct from the Central Range. To do this, we combined previously published cytochrome b (cyt b) sequences from S. bombifrons populations sampled within the Central Range (Rice & Pfennig, 2008) with newly generated cyt b sequences from 13 S. bombifrons sampled within the South Texas range. Because we no longer had access to tissue or DNA from the individuals sampled for Rice & Pfennig (2008), sequencing additional markers would only have been possible in our novel South Texas population, and therefore would not have been informative. While sequence analysis of a single mitochondrial marker is not ideal for the characterization of population genetic structure (Toews & Brelsford, 2012), the available genetic data provide an independent but complementary analysis of the S. bombifrons range, together with our species distribution models. We used our expanded S. bombifrons cyt b sequence dataset to construct a haplotype network, calculate pairwise Phi_ST values among all populations, and test for signatures of isolation by distance.

To add the South Texas S. bombifrons population to previously sampled Central Range populations, we obtained tissue samples from 13 new S. bombifrons individuals. Toe clip tissue samples were collected from 11 juvenile S. bombifrons in May–June 2013 (Brooks County, TX, USA; collected by Gerry Salmon under permit SPI 1097-912) and from one adult S. bombifrons in September 2013 (Hidalgo County, TX, USA; collected by Seth Patterson under a valid Texas State hunting license with Herp Stamp). Toe clip tissue samples were preserved in 70% ethanol. We also obtained one additional Hidalgo County, Texas S. bombifrons tissue sample from the University of Texas at Arlington’s Amphibian and Reptile Diversity Center.

We grouped the Central Range sequences from Rice & Pfennig (2008) plus the new South Texas sequences into 12 populations, each defined by sampling locations within 200-km diameter circles, defined using ArcGIS 10.8.2 (Table S1). Unlike all the other populations, where adult or juvenile toads were sampled, the sequences from the Southeast Arizona population (Rice & Pfennig, 2008) came from tadpoles. Because the vast majority of tadpoles will not survive to reproduce, including all 191 (Rice & Pfennig, 2008) Southeast Arizona cyt b sequences would likely provide an inaccurate representation of population genetic variation. We therefore subsampled this population as follows. The tadpoles had been collected from 16 temporary pond sites, with a maximum of three different haplotype sequences in any individual pond. We chose four samples from each of the pond sites, making sure that we sampled every haplotype. If a pond contained sequences from only a single haplotype, all four samples that we chose had that haplotype. If a pond had three different haplotypes, we sampled all three and chose the fourth from whichever haplotype had the highest frequency among that pond’s samples. We therefore used 64 of the 191 available cyt b haplotype sequences for the Southeast Arizona population. The number of samples from the other 11 populations ranged from 3 to 51, giving us a total of 202 S. bombifrons cyt b sequences (Table S1). Detailed information about all collection locations across both the Central Range and the South Texas range, along with any museum catalog and GenBank accession numbers, is available in Table S1.

We extracted DNA from the new South Texas tissue samples using Qiagen DNeasy Blood and Tissue Kits (Qiagen, Hilden, Germany), following the manufacturer’s protocol. Next, we amplified and sequenced 663 bp of the mitochondrial cyt b gene. We used a forward primer designed from a S. multiplicata cyt b sequence (SCB1-F: 5′-TCCCAACCCCATCTAACATC-3′) and a reverse primer designed from a Xenopus laevis cyt b sequence (XCB2-R: 5′-GAGGGCTAGATTAGGATGGATA-3′), both previously published (Rice & Pfennig, 2008). For PCR amplification, we followed the protocol from Rice & Pfennig (2008), using an Eppendorf Mastercycler Pro gradient thermal cycler (Eppendorf, Hamburg, Germany). After PCR amplification, we purified the PCR products with ExoSAP-IT (Affymetrix, Santa Clara, CA, USA), and submitted sequencing reactions to GENEWIZ (South Plainfield, NJ, USA) for Sanger sequencing. We aligned the new South Texas S. bombifrons cyt b sequences together with the previously published cyt b sequences (Rice & Pfennig, 2008; Table S1) from the 11 other populations using MEGA v6 (Tamura et al., 2013). The new South Texas S. bombifrons cyt b sequences were submitted to GenBank, and the resulting accession numbers are listed in Table S1.

To determine whether South Texas is genetically distinct from the Central Range populations, we used the pegas package (Paradis, 2010) in R to construct a haplotype network using the Parsimony network method (Templeton, Crandall & Sing, 1992). We then used the haplotypes package (Aktas, 2020) in R to calculate pairwise Phi_ST values for our twelve populations, including 11 within the Central Range and 1 from South Texas. With the ape R package (Paradis & Schliep, 2019), we performed a principle coordinate analysis (PCoA) using the pairwise Phi_ST values (Borokini, Klingler & Peacock, 2021) to determine whether the South Texas population showed patterns of genetic differentiation from the 11 Central Range populations that would be consistent with population disjunction.

Finally, to assess whether the South Texas population showed greater pairwise genetic differentiation per unit of pairwise geographic distance than the 11 Central Range populations, which would also be consistent with a range disjunction, we assessed isolation by distance relationships both including and excluding all pairwise comparisons involving South Texas. We calculated pairwise geographic distances between the centers of the 200-km diameter circular sampling areas defining all 12 populations by measuring the straight-line distance between each pair of populations in ArcGIS 10.8.2 (Table 1). We then used the vegan package in R (Oksanen et al., 2022) to run Mantel tests on the pairwise geographic distance and pairwise Phi_ST distance matrices, using the Pearson correlation method and 9,999 permutations. All R analyses were done using version 4.4.1 (R Core Team, 2024). The cyt b sequence data and R code are provided as Supplemental Files.

Geographic distribution

After removing imprecise locality records, we were left with 5,583 S. bombifrons records from 1,613 unique localities from museums plus an additional 882 observational records from iNaturalist. The distribution map captures the described range of S. bombifrons, from southern Canada into northern Mexico and clearly shows a marked range disjunction between the Central Range and South Texas (Fig. 2A). In the comparison species group, the 20 selected species resulted in 18,323 records covering the entire range of S. bombifrons and the area of disjunction (Fig. 2B). Ten of the 20 species were found across the entire region. Nine species either fell both in the Central Range and region of disjunction (6/20) or in the South Texas range and region of disjunction (3/20). Only one species was restricted to a single region (the Central Range) (Table S2). With locality records from 19 of the 20 comparison species found between the South Texas and Central Range, we conclude that the absence of S. bombifrons is not due to a lack of sampling effort or a geographic barrier, but rather represents a true range disjunction. Furthermore, the fact that no other species shows a similar pattern of disjunction strongly suggests the range disjunction is not due to a common natural barrier to gene flow across many species; it is unique to S. bombifrons.

Species distribution modeling

After cleaning and spatially rarefying the museum locality data to the more precise time frame (1950–2015) and degree of precision (within 5 km), we were left with 923 unique sampling localities: 18 from South Texas and 905 from the Central Range region.

Current models

To evaluate model performance, we compared both AUC values and the omission rate of the independent samples across all models. The regularization multiplier determines how tightly fit the model is to the known presence localities; lower regularization multiplier values produce more tightly fit models which can lead to underprediction of suitable habitat (e.g., overfitting) (Phillips, 2017). Here, we found lower regularization multipliers resulted in higher AUCs and all models correctly identified the presence of most (at least 54 of the 63 points) of the independent population genetic samples and at least 397 out of the 513 points from iNaturalist that remained after spatial rarefication (Table S3). All models confirm that the core of the range is in the Central Range region, with a much smaller region of habitat suitability in South Texas (Figs. 3A–3C). This disjunction is clear across all regularization multipliers; at the default regularization multiplier (1.0), there is a complete gap between these populations where no moderate or high-quality habitat is found (Fig. 3A). This disjunction is even more obvious at the lowest value of the regularization multiplier (0.25); this was also the only model using all locality points to show any high-quality habitat in the south Texas region (Fig. 3B). Even much higher regularization multipliers (3.0) show only lower quality habitat in the region of disjunction (Fig. 3C). As predicted by the literature (Elith et al., 2011; Kramer-Schadt et al., 2013), the reduced variable model showed no substantive differences from the full model (Fig. S1).

Modeling Central Range (Fig. 4A) and South Texas (Fig. 4B) locality records separately further confirms that there is little overlapping habitat space between these two distinct populations. Using only locality data from the Central Range shows no suitable habitat within the entire South Texas region, and using only locality data from South Texas shows no suitable habitat within the Central Range. Additionally, both models also show only low-quality habitat throughout the region of disjunction. As models trained on localities within one region do not identify suitable habitat within the other, this indicates that the environmental conditions between the two regions are distinct and non-overlapping. This finding is further confirmed by our principal components analysis, which shows the two populations do not overlap in environmental space (Fig. 5).

When evaluating the effect of individual variables on the model results, soil was the most important variable when using localities from across the geographic range. Soil had an importance value of 25.4% when including all environmental variables and 32.4% in the reduced variable model (Table 2). The regularization multiplier did not have a noticeable impact on the variable contribution. When only the Central Range localities were included, soil was the second most important variable (at 16.1%) while maximum temperature of the warmest month was the most important (20.6%). In contrast, when only South Texas localities were included, elevation was overwhelmingly the most important variable at 71.7%.

Table 2:

Variable importance for each of the models run with a regularization multiplier of 1.0.

Variable	Importance value (%) total range	Importance value (%) reduced model	Importance value (%) Central Range	Importance value (%) South Texas
Soil	25.4	32.4	16.1	7.0
Bio 5–max temp of warmest month	21.8		20.6	0.6
Bio 8–mean temp wettest quarter	13.1		12.2	6.5
Bio 10–mean temp of warmest quarter	11.1	38.6	14.7	2.7
Sr 6–solar radiation June	9.2	12.3	15.4	2.0
Sr 7–solar radiation July	8.9	7.7	6.6	1.7
Bio 13–precipitation of wettest month	5.9	6.2	4.8	1.5
Elevation	2.7	2.8	7.4	71.7
Bio 18–precipitation of warmest quarter	1.2		1.6	2.8
Bio 16–precipitation of wettest quarter	0.7		0.7	3.5

DOI: 10.7717/peerj.17968/table-2

Future models

Models using future climate conditions show a significant range retraction across the entire range when all localities are used in training the model across all regularization multipliers (Figs. 6A–6C). The core region of best habitat shows a northward shift with a reduction in overall area, and moderately suitable habitat also shows an extreme reduction. At the default regularization multiplier (1.0), there is the complete loss of any suitable habitat in the South Texas region (Fig. 6A). For both lower (0.25, Fig. 6B) and higher (3.0, Fig. 6C) regularization multipliers, only very small, isolated patches of moderately suitable habitat remain. These models suggest future extirpation of the entire South Texas population.

Population genetic analyses

The haplotype network (Fig. 7) revealed that the South Texas haplotypes are unique to that population, a result that is consistent with the geographic disjunction shown in the species distribution modeling results. The pairwise Phi_ST values (Table 3) also suggest that South Texas is genetically distinct from the other subpopulations. The PCoA based on these pairwise Phi_ST values shows that South Texas is differentiated from 10 of the Central Range populations along Axis 1 (Fig. 8), which had an eigenvalue equal to 0.68. Further, South Texas is differentiated from the remaining Central Range population (Southeast Arizona) along Axis 2 (Fig. 8), which had an eigenvalue equal to 0.31.

Our isolation by distance analyses suggest that South Texas shows greater genetic differentiation per unit of geographical distance than the Central Range populations, which is also consistent with the geographic disjunction shown by the species distribution modeling. Signatures of isolation by distance were significant, both when including (Mantel r = 0.72, p = 0.005) and excluding (Mantel r = 0.58, p = 0.04) pairwise comparisons involving South Texas. The slope of the line illustrating the correlation between geographic distance and genetic distance is steeper when the comparisons involving South Texas are included (Fig. 9). Further, the South Texas comparisons all fall above the line that illustrates the isolation by distance correlation among only the Central Range populations (Fig. 9).

Evaluating the probability of a true Spea bombifrons range disjunction

Here, we addressed two key questions: 1) Is the South Texas Spea bombifrons population geographically and environmentally disjunct and genetically distinct from the rest of the range? and 2) How will climate change affect the future distribution of this species? We find that all lines of evidence, including the distribution of known locality data (Fig. 2A), the species distribution modeling results (Fig. 3), and the population genetic analyses (Figs. 7–9, Table 3), strongly suggest that the South Texas population of S. bombifrons is completely isolated from the Central Range of the species, consistent with a true range disjunction. The distribution of other amphibian species from throughout the region demonstrate that the gap between populations has been well sampled, so the lack of S. bombifrons records from this region is not the result of under-sampling (Fig. 2B). Further, current climate models (Fig. 3A) and Principal Components Analyses (Fig. 5) suggest this geographic disjunction is likely linked to the environment. Consistent with these findings, our population genetic analyses indicate that S. bombifrons in South Texas are genetically distinct from populations in the Central Range (Figs. 7–9, Table 3).

Current climate models show high habitat suitability throughout the Great Plains, into southern Canada, and in south Texas. Maxent models suggest a very limited potential for connectivity between the two populations if all localities from both the Central Range and South Texas are used (Fig. 3). This lack of connectivity is even more obvious when populations are modeled independently (Fig. 4). Models using both more and less restrictive regularization multipliers show similar patterns (Figs. 3B and 3C).

Several biologically important variables contributed to these results. Bulk soil density was particularly important as it had the highest percent contribution to the models (Table 2). Although soil variables can improve model performance in woody plant species, soil variables still remain rarely used in species distribution studies (Hageer et al., 2017). Given the important contribution of the soil variable in our results (Table 2), we would reiterate Hageer et al.’s (2017) suggestion that soil variables be more widely adopted in species distribution modeling, even for animal studies. The inclusion of soil variables is likely to be particularly important for fossorial animal species such as Spea bombifrons.

Several independent population genetics analyses are consistent with the disjunction of the South Texas population. The South Texas cytochrome b haplotypes are unique to this population (Fig. 7), and this population is differentiated from the remaining populations (Fig. 8, Table 3). Also, the isolation by distance analyses suggest that for the cyt b locus, South Texas shows greater genetic differentiation per unit of geographical distance than the Central Range populations (Fig. 9), again consistent with population disjunction.

As none of the other frog and toad species we examined show a similar pattern of disjunction (Fig. 2B), the isolation of the South Texas population likely occurred due to a rare long distance colonization event rather than a major habitat change causing an overall range retraction and a residual relict population. However, Gherghel & Martin’s (2020) study of glacial refugia in Spea spp. suggested that there may have been multiple refugia for S. bombifrons. Our work does not allow us to rule out the possibility that ancient habitat change caused a range retraction specific to S. bombifrons.

Predicting climate change impacts on future species distribution

Projecting the future range of S. bombifrons based on climate models shows troubling results. Under even the conservative representative concentration pathway chosen for future greenhouse gas emissions, range retractions are predicted across the entire range of the species. Furthermore, the isolated South Texas population will likely be severely reduced or even extirpated entirely (Fig. 6). While many species are responding to climate change by shifting ranges northwards (Parmesan, 2006; Angert et al., 2011; Chunco, 2014), the limited dispersal ability of this species combined with the unsuitable habitat underlying the range disjunction makes it unlikely that the South Texas population will be able to respond to climate change via dispersal.

As anthropogenic forces contribute to increasing range disjunction in many species as a result of habitat fragmentation and loss, predicting how individual populations will respond to climate change is increasingly important. The South Texas population of spadefoot toads experiences a significantly different climate than other populations and has two unique haplotypes not seen in other populations (Fig. 7). Future studies that attempt to detect signatures of selection across the entire genome for S. bombifrons from South Texas could provide some insight into whether and how local adaptation of the South Texas population may have contributed to its survival in this region.

Understanding the geographic distribution of a species and the cause of any disjunction is critically important in predicting the response of that species to climate change. For example, anthropogenically driven disjunctions due to land use conversion and the resulting habitat fragmentation can greatly limit a species’ ability to respond to climate change by restricting available habitat to shift into and because isolated populations have lower standing genetic variation and thus a decreased ability to adapt to climate change. Significant literature has focused on the impacts of anthropogenic habitat fragmentation on the abilities of species to respond to climate change (e.g., Honnay et al., 2002; Opdam & Wascher, 2004).

Natural (non-anthropogenic) range disjunctions can also be important considerations in evaluating and predicting a species’ response to climate change. For example, climate change will proceed at a different rate for populations occurring at different latitudes or elevations, and population isolation means that any given population will have a lower effective population size than an overall census would suggest. Thus, the combination of limited geographic range, low gene flow between populations, and demographic stochasticity within populations may mean that the extinction risk of disjunct populations is higher than currently recognized.