We collected tissue samples from 113 individuals among 32 sites across the Cumberland Plateau and Appalachian Valley and Ridge of central and eastern Tennessee and northeastern Alabama, USA (Table 1; Figure 1). Tissue samples were obtained from captured salamanders and comprised of tail tips (<5 mm) that were immediately preserved in 95–100% ethanol, then stored at −20°C at The University of Alabama in Huntsville until extraction.

We generated novel mtDNA sequence data for 106 individuals (nd4: GenBank accession nos. OM111022–OM111127; cytb: OM110924–OM111021). We also supplemented our novel mtDNA sequence library for both nd4 and cytb with datasets available on GenBank and supplemental data available from Johnson (2002) (AF489614–AF489629) and Patton et al. (2019) (OM678465–OM678536). The nd4 dataset included 42 haplotypes from 122 individuals among 39 sites in Alabama, North Carolina, Tennessee, and West Virginia. The best models of nucleotide substitution for first, second, and third codon positions were Hasegawa-Kishino-Yano + Γ (HKY + G), Hasegawa-Kishino-Yano + invariant sites (HKY + I), and Tamura-Nei + Γ (TrN + G), respectively. Both ML and Bayesian inference identified two major reciprocally monophyletic lineages in Tennessee and northern Alabama (Supplemental Figure S1): a southern lineage comprised of seven haplotypes from 60 individuals and 17 sites in northern Alabama and south-central Tennessee and a northern lineage comprised of 35 haplotypes from 61 individuals and 15 sites in Tennessee, as well as six sites in West Virginia (Johnson, 2002). The southern lineage includes populations from the southern Highland Rim and Cumberland Plateau of northeastern Alabama and south-central Tennessee (Figure 1). The northern lineage includes populations from the Appalachian Plateau of West Virginia and the northern Cumberland Plateau, Eastern Highland Rim, and Valley and Ridge of north-central Tennessee (Figure 1). A contact zone between the northern and southern lineages based on the nd4 dataset occurs along the Cumberland Plateau and Walden Ridge along a line running from south-central Warren and Van Buren counties eastward to central Bledsoe and Rhea counties.

The cytb dataset included 63 haplotypes from 170 individuals and 97 sites in Alabama, Georgia, Kentucky, Ohio, North Carolina, Pennsylvania, South Carolina, Tennessee, Virginia, and West Virginia. The best models of nucleotide substitution for first, second, and third codon positions were Jukes-Cantor + Γ (JC + G), Jukes-Cantor (JC), and Hasegawa-Kishino-Yano (HKY), respectively. Both ML and Bayesian inference again identified two major reciprocally monophyletic lineages in Tennessee and northern Alabama (Figure 2): a southern lineage comprised of 16 haplotypes from 80 individuals and 43 sites, and a northern lineage comprised of 16 haplotypes from 73 individuals and 45 sites. Topologies of the cytb ML and Bayesian genealogies were very similar to those presented in Patton et al. (2019) and the nd4 dataset. The southern lineage includes populations from the southern Highland Rim, Cumberland Plateau, and Valley and Ridge of northeastern Mississippi, northern Alabama, northwestern Georgia, and south-central Tennessee (Figure 1). The northern lineage includes populations from the Appalachian Plateau and Valley and Ridge of Kentucky, Ohio, Pennsylvania, Virginia, and West Virginia southward through the northern Cumberland Plateau, Eastern Highland Rim, and adjacent Valley and Ridge of north-central Tennessee (Figure 1). A contact zone between the northern and southern lineages based on the cytb dataset occurs on the Cumberland Plateau along a line from northern Van Buren and Bledsoe counties into central Rhea County.

For the nd4 dataset, the ABGD approach recovered three MOTUs with convergence of initial and recursive partitions occurring at prior intraspecific divergence (P) = 0.009 (Supplemental Figure S2). The three MOTUs corresponded to the major lineages in the nd4 phylogeny (Supplemental Figure S1): southern lineage, northern lineage, and Blue Ridge (the single sequence from Henderson Co., North Carolina included). Initial partitions stabilized at P = 0.002, and both initial and recursive partitions remained stable through P = 0.062. The mPTP approach identified four MOTUs: the same three identified by the ABGD approach as well as the population from No Business Branch in Campbell Co., Tennessee (site 30).

For the cytb dataset, the ABGD approach resulted in four MOTUs with convergence of initial and recursive partitions occurring at P = 0.007 (Supplemental Figure S2). The four MOTUs corresponded to the major lineages in the cytb phylogeny (Figure 2) with a few key differences. The southern lineage and A. caryaensis both were delimited as distinct MOTUs, which is consistent with the cytb phylogeny; however, the northern and Blue Ridge lineages were delimited as single MOTU and sample Aa70TN from Marion Co., Tennessee (site 85) was delimited as a distinct MOTU. Sample Aa59TN from Bledsoe Co., Tennessee (site 75) grouped in the northern + Blue Ridge lineage. Initial partitions stabilized at P = 0.003, and both initial and recursive partitions remained stable through P = 0.011. The mPTP approach delimited six MOTUs: the southern lineage and A. caryaensis were delimited as distinct MOTUs and the northern + Blue Ridge lineages were delimited as a single MOTU as in the ABGD analysis. Samples Aa70TN from Marion Co., Tennessee (site 85), Aa76TN from Franklin Co., Tennessee (site 89), and Aa59TN from Bledsoe Co., Tennessee (site 75) each were delimited as distinct MOTUs. The cytb sequence for sample Aa76TN includes a significant amount of missing data, likely influencing its delimitation as a distinct MOTU.

Nearly all measures of genetic diversity were higher for the northern lineage compared to the southern lineage for both mitochondrial loci (Table 3). Tests for population size changes based on summary statistics showed a trend toward population expansion in both the northern and southern lineages (Table 3). Mean pairwise distance (p-distance) between the southern and northern lineages is 8.7 ± 0.9% for nd4 and 8.5 ± 1.2% for cytb. Mean distances (±SE) within the southern lineage for nd4 and cytb were 0.3 ± 0.1% and 1.5 ± 0.3%, respectively, while mean distances (±SE) within the northern lineage for nd4 and cytb were 1.0 ± 0.2%, and 1.7 ± 0.3%, respectively.

A total of 451,117 SNPs were obtained following SNP calling in ipyrad. After filtering, 33,423 unlinked SNPs were included in the final SNP dataset for 83 individuals. We estimated overall F_st at 0.56, indicative of strong population structuring for the entire dataset. This was supported by fastStructure results, which revealed the best fit of K = 2 clusters (Figure 3). Remarkably, there is complete isolation between the two clusters, with all individuals exhibiting probability of assignment nearly equal to 1. We also conducted independent fastStructure runs (K = 1–7) within each cluster. In both cases, the best fit was K = 1, indicating no additional substructure among populations within lineages and substantial gene flow within each of the two clusters, even among geographically distant sites.

The DAPC had 113 unlinked SNPs and estimated two clusters (K = 2) based on BIC. The optimal number of PCs retained for the DAPC estimated by the α-score optimization was one. The clusters estimated by the DAPC were consistent with a southern and northern group. Cluster 1 included all samples from the southern sites (Jackson Co. in Alabama and Franklin, Marion, and Grundy counties in Tennessee; sites 1–10, 12, 14–16) except one sample from Ozone Falls SNA in Cumberland County, Tennessee (site 26) and one sample from Catoosa WMA in Roane County (site 28). Both individuals were included in the northern group in other analyses (fastStructure and IQ-TREE). Cluster 2 included the remaining samples from the northern sites in Cannon, Van Buren, White, Cumberland, Roane, Fentress, Anderson, Campbell, and Sullivan counties in Tennessee (sites 19–32; Figure 3).

Two well-supported, reciprocally monophyletic lineages were identified based on the maximum-likelihood IQ-TREE analysis (Figure 4): a southern clade including populations from Jackson County, Alabama as well as populations in Franklin, Marion, and Grundy counties in south-central Tennessee; and a northern clade including populations from the Valley and Ridge, northern Cumberland Plateau, and northern Eastern Highland Rim as far south as Cannon and Van Buren counties in Tennessee (Figure 1). The nuclear SNP phylogeny is similar to the mitochondrial nd4 phylogeny (Supplemental Figure S1), and there was no discordance in population assignment to major lineage (southern or northern) for populations included in both the mitochondrial nd4 and nuclear SNP phylogenetic analyses.

Previous studies of Aneides species in western North America have discovered considerable cryptic genetic structure and diversity (Reilly et al., 2013, 2015; Reilly and Wake, 2015). Likewise, studies have uncovered potential cryptic diversity within the Castaneides species complex. Sessions and Kezer (1987) examined karyotypic variation in plethodontid salamanders, including the genus Aneides. At least three karyotypes occur in the Castaneides species complex, with two karyotypes reported originally in sympatry at one site “near Chattanooga, Tennessee.” However, the report of sympatry was an error and explained in Patton et al. (2019). Three specimens were collected from Buck Creek Cove in Franklin Co., Tennessee; one of these specimens was examined by Sessions and the other two sent to and examined by Kezer. Kezer reported that the two specimens examined had predominately bi-armed chromosomes (aeneus II karyotype in Sessions and Kezer, 1987). When Sessions double-checked his lab notes, he found that the specimen he examined also had the aeneus II karyotype and was not aeneus I (Sessions, personal communication to David Wake in Patton et al., 2019). The aeneus II karyotype also was found in northern Alabama (Sessions and Kezer, 1987). The aeneus II karyotype may be diagnostic of the southern lineage but additional karyological is needed to confirm this hypothesis. The aeneus I karyotype was identified from specimens collected in southern Pennsylvania, West Virginia, southwestern Virginia, and southern Kentucky (Sessions and Kezer, 1987) and may be diagnostic for the northern lineage. However, this karyotype also was detected in southwestern North Carolina from Hickory Nut Gorge in what is now recognized as A. caryaensis (Patton et al., 2019), which is likely the sister group to the northern lineage. Morescalchi (1975) reported a potential third distinct karyotype (aeneus III in Sessions and Kezer, 1987) from Tennessee. Unfortunately, detailed locality information is lacking (reviewed in Patton et al., 2019) but is suspected to have been collected from “a probable site somewhere in Cumberland County, Tennessee” by an unknown supplier (Sessions and Kezer, 1987). Based on historical occurrences, if this collection indeed was from Cumberland County, it most likely was collected from Ozone Falls (site 26), which grouped within the northern lineage in our analyses.

More recently, Patton et al. (2019) integrated population genetic, genomic, and species delimitation approaches of mitochondrial, nuclear microsatellite, and nuclear SNP loci with morphological analyses and found strong support for recognition of four lineages. In addition to A. caryaensis, which was described concurrently in their study, the three delimited lineages include: (i) a southern lineage from the southern Cumberland Plateau of south-central Tennessee, northeastern Mississippi, northern Alabama, and northwestern Georgia, (ii) a northern lineage from the Cumberland Plateau, Cumberland Mountains, Appalachian Plateau, and Appalachian Valley and Ridge from north-central Tennessee northward into eastern Kentucky, southwestern Virginia, West Virginia, southern Ohio, and southwestern Pennsylvania; and (iii) a lineage from the Blue Ridge Escarpment of western North Carolina, western South Carolina, and northeastern Georgia (Figure 2A in Patton et al., 2019). These lineages are allopatrically distributed, except northern and southern lineages may overlap to some degree on the Cumberland Plateau of Tennessee based on mtDNA analyses. Geographic sampling across the potential distributions of northern and southern lineages was limited (n = 4 in Virginia and n = 2 in Alabama) for population genomic analyses, but results were consistent with mtDNA species delimitations (Patton et al., 2019). The authors refrained from describing these three lineages but recommended they be recognized as ESUs in light of these findings until additional morphological and population genomic studies were conducted.

Results of our phylogenetic, population genomic, and species delimitation approaches of mitochondrial and nuclear loci are consistent with the lineages delimited by Patton et al. (2019). Our mtDNA genetic and species delimitation analyses, which included the cytb dataset of Patton et al. (2019), identified the northern and southern lineages with strong support. These lineages were also strongly supported in our nuclear genomic analyses (i.e., fastStructure and IQ-TREE). The contact zone between the northern and southern lineages based on the nd4 dataset occurs along the Cumberland Plateau and Walden Ridge along a line running from south-central Warren and Van Buren counties eastward to central Bledsoe and Rhea counties in Tennessee. The putative contact zone derived from the cytb dataset, which included sequences from Patton et al. (2019), is similar to that of the nd4 dataset. However, a population from Morgan Co., Tennessee (site 68; Aa41TN in Patton et al., 2019) grouped with the southern lineage. Two additional sequences from Patton et al. (2019)—Aa79TN (site 92) and Aa77TN (site 90)—were discordant with lineage based on geographic proximity. Aa79TN from Cannon Co., Tennessee did not group with other samples we collected nearby at Short Mountain (site 20) within the northern lineage, while Aa77TN from Hamilton Co., Tennessee did not group with samples collected nearby (site 13) that grouped within the southern lineage. When the data from Patton et al. (2019) are excluded, we found no evidence of discordance between geography and lineage. Our nuclear genomic analyses also found little evidence of discordance between geography and lineage. Likewise, there was little discordance in population assignment to lineage for populations included in both the mitochondrial (excluding data from Patton et al., 2019) and nuclear SNP analyses. The exception was two individuals—one from Ozone Falls SNA in Cumberland County (site 26) and the other from Catoosa WMA in Roane County (site 28)—that based on geography, mtDNA, fastStructure, and IQ-TREE analyses grouped within the northern lineage but grouped within the southern lineage in the nuclear DAPC analysis. This discordance may be an artifact of missing data, which was 98% for the individual from Ozone Falls SNA and 62% for the individual from Catoosa WMA.

Patton et al. (2019) hypothesized that climatic and ecological changes associated with the Appalachian Miocene uplift 20–10 mya (Gallen et al., 2013; Liu, 2014) may have led to the formation of at least two refugia after the Miocene optimum 15 mya: one in the Hickory Nut Gorge area and the other on the southern Cumberland Plateau that later gave rise to the northern, southern, and Blue Ridge lineages. Based on the divergence time analysis of Patton et al. (2019), the southern lineage split from the ancestor of the Blue Ridge and northern lineages ~9 mya, while the Blue Ridge and northern lineages diverged about 5 mya at the onset of the Pliocene. Climatic and ecological shifts during the Pliocene and Pleistocene likely resulted in concurrent range shifts in the northern and southern lineages whereby salamanders retreated to cooler, mesic forested coves along escarpments during warmer, drier periods and expanded during cooler, wetter periods. Demographic analyses support the existence of two refugia since the late Miocene and subsequent episodes of range expansion and contraction leading to the present-day distributions of the northern and southern in the Cumberland Plateau and Appalachian Valley and Ridge.

No obvious barriers to dispersal exist along the escarpments of the Cumberland Plateau and Walden Ridge in Tennessee between the northern and southern lineages. Interestingly, a contact zone in the southern Cumberland Plateau has been documented in some other taxa. In the Mountain Dusky Salamander (Desmognathus ochrophaeus) species complex, the Cumberland Dusky Salamander (D. abditus) occurs, albeit patchily distributed, throughout much of the southern Cumberland Plateau in Tennessee, contacting the Mountain Dusky Salamander (D. ochrophaeus) just south of the Cumberland Mountains section of the Cumberland Plateau in Morgan County (Anderson and Tilley, 2003; Drukker et al., 2018), which is ~50–60 km north of the contact zone between the northern and southern lineages in Aneides. Previous studies have noted the similarities in distribution and habitat between Aneides salamanders and the lampshade spiders of the genus Hypochilus (Catley, 1994; Hedin, 2001; Patton et al., 2019). Hypochilus thorelli ranges along the escarpments of the Cumberland Plateau from northeastern Alabama and northwestern Georgia northward through Tennessee into southeastern Kentucky and southwestern Virginia (Hedin and Wood, 2002; Keith and Hedin, 2012) where it can be found alongside A. aeneus on exposed sandstone and carbonate outcrops in mesic forests (Beatty pers. comm. in Fergusson, 1972). Hedin and Wood (2002) conducted a population genetic study of H. thorelli based on the mitochondrial cytochrome oxidase subunit I (COI) locus across 19 sites finding substantial genetic structure and evidence for very limited female-based dispersal, even among some sites within 5 km. Although spiders were morphologically similar across sites, three mitochondrial clades were identified: (i) a northern clade in the Cumberland Mountains of Lee Co., Virginia and Campbell and Morgan cos., Tennessee; (ii) a central + Kentucky clade from the central Cumberland Plateau in Tennessee from southern Cumberland, Bledsoe, and Van Buren counties but also two sites in Pickett County, Tennessee and Whitley Co., Kentucky; and (iii) a southern clade from the Cumberland Plateau of northeastern Alabama, northwestern Georgia, and southern Tennessee, including Franklin, Grundy, Marion, Hamilton, and Rhea counties. However, divergence times among H. thorelli clades appear to be younger than those of the Castaneides species complex estimated by Patton et al. (2019), with splits occurring within the last two million years (Ciaccio, 2020). Interestingly, the contact zone between the southern and central H. thorelli clades (Hedin and Wood, 2002) is strikingly similar to the contact zone between the northern and southern lineages in Aneides, suggesting similar biogeographical factors, such as historical cycles of climate change during the Pleistocene, have shaped present-day genetic structure and distributions in these co-occurring species.

In addition to describing A. caryaensis, Patton et al. (2019) recommended the recognition of the Blue Ridge lineage as an ESU to facilitate its conservation. However, the authors noted that species delimitation methods grouped Blue Ridge populations with the northern lineage, which also was consistent with results of our mtDNA species delimitation. We did not sample Blue Ridge populations in our nuclear SNP dataset and recommend that such analyses be conducted before recognizing Blue Ridge populations as a distinct ESU from the northern lineage. Patton et al. (2019) also noted that mtDNA phylogenetic and species delimitation approaches supported the distinction of the southern lineage, but they refrained from recommending its designation as an ESU until additional investigations were conducted. Our mitochondrial and nuclear SNP analyses unequivocally support the recognition of the southern lineage as an ESU sensu Moritz (1994) genetically distinct from the northern lineage. We found no evidence of admixture between northern and southern lineages including for nuclear loci, which recovered reciprocal monophyly for the southern lineage. The recognition of the southern lineage as an ESU is consistent with several other proposed definitions and criteria including Avise (1994); Vogler and Desalle (1994); USFWS (U.S. Fish Wildlife Service) NMFS (National Marine Fisheries Service) (1996); Fraser and Bernatchez (2001), and Funk et al. (2012). Moreover, we fully expect that the southern lineage will meet other definitions of ESU (e.g., Dizon et al., 1992), likely constituting an undescribed species pending additional morphological study.

Both the southern and northern lineages in the Cumberland Plateau face a plethora of threats that are currently or may in the future impact populations (Niemiller et al., 2020), particularly habitat loss and fragmentation from logging, construction of open corridors (e.g., roads and power line rights-of-way), and past and current mining operations. Other significant threats include water impoundments, over-collection associated with the pet trade, emergent amphibian diseases, and climate change. Life history traits, such as slow maturity and reduced fecundity (Waldron and Pauley, 2007), increase intrinsic vulnerability to threats. Unfortunately, the ecology and life history of populations of both the southern and northern lineages on the Cumberland Plateau are poorly understood (Miller and Reynolds, 2011; Niemiller et al., in press). In particular, very limited historical data is available for most populations to infer population trends. We do not know whether populations on the southern Cumberland Plateau have experienced declines similar to those reported for A. caryaensis and the Blue Ridge Escarpment populations (Snyder, 1983, 1991; Corser, 2001; Patton et al., 2019). Encouragingly, recent surveys in Tennessee and northeastern Alabama confirmed the presence of species at 16 historical sites and discovered 34 new sites (Niemiller et al., 2020, in press).