The global decline of amphibians as a result of rapid growth of human populations and their expansive impact on natural environments is well documented (Petranka, 1998; Collins and Storfer, 2003; Becker et al., 2007), with many species facing rapid declines and extinction (Pechmann and Wilbur, 1994; Alford and Richards, 1999; Stuart et al., 2006; Wake, 2012). A major challenge in worldwide amphibian conservation efforts is that many amphibian species are characterized by a highly conserved morphology (Cherty et al., 1978; Wake, 1991; Stuart et al., 2006; Kozak and Wiens, 2010). This morphological conservatism in turn exacerbates the challenge of delimiting species, thus preventing meaningful assessment of their true conservation status. Inaccurate or incomplete taxonomy may exacerbate population declines and elevate extinction risks through the lack of necessary conservation interventions.

Species delimitation is one of the most difficult and controversial subfields of ecology and evolutionary biology (De Queiroz, 2007; Wiens, 2007), yet given trends in defaunation (Dirzo et al., 2014), the need for characterizing extant diversity is growing. Modern molecular tools, a more unified species concept (i.e., the generalized lineage concept; De Queiroz, 2007), and tree-based delimitation methods have led to the identification of a plethora of unique taxa from cryptic complexes, establishing them as invaluable tools in this endeavor (Bogart and Tandy, 1976; Hillis et al., 1983; Highton et al., 1989; Wynn and Heyer, 2001; Jockusch and Wake, 2002; Pauly et al., 2006; Elmer et al., 2013).

A common approach to facilitate species delimitation is to leverage phylogenetic species delimitation methods (e.g., Roy et al., 2014; Kotsakiozi et al., 2018; Titus et al., 2018). Based on, in part, the phylogenetic species concept, these methods use phylogenetic and coalescent models to guide the delimitation of species from empirical phylogenies (e.g., Pons et al., 2006; Yang and Rannala, 2010; Zhang et al., 2013; Yang, 2015). However, methods leveraging the multispecies coalescent have recently been demonstrated to often delimit population structure, not species (Sukumaran and Knowles, 2017). To avoid this concern, at least in part, methods reliant on models other than the multispecies coalescent may be used (e.g., Poisson Tree Processes [PTP; Zhang et al., 2013]). Species delimitation necessitates the use of multiple lines of evidence (e.g., molecular, morphological, ecological, and biogeographical). Although no single phylogenetic species delimitation method is invulnerable to the concerns voiced by Sukumaran and Knowles (2017), these approaches serve as a valuable complement to other sources of evidence when describing cryptic diversity.

Cryptic salamander complexes frequently occur in areas that contain a rich variety of climatic zones, geologic formations, and habitat diversity (e.g., California Floristic Province [Myers et al., 2000; Lapointe and Rissler, 2005; Rissler et al., 2006; Rissler and Apodaca, 2007; Reilly and Wake, 2015] and eastern and southern Mexican highlands [Rovito et al., 2013]). The Appalachian Mountains of the eastern United States display many of these attributes and are one of the oldest continuously exposed land masses on earth. They have a rich geologic history (Pickering et al., 2003), a diverse set of climatic regions, and a multitude of rare ecosystem types, many of which are globally threatened (Noss et al., 1995). The southern Appalachians are especially rich in biodiversity and endemism (Pickering et al., 2003) and are a center of diversity for many taxa.

Notably, salamanders of the family Plethodontidae are extremely diverse in this region and have high levels of both endemicity and cryptic diversity (Camp and Wooten, 2016). A classic example of non-adaptive radiation (Rundell and Price, 2009), salamanders of the genus Plethodon diversified rapidly in the eastern United States, yet are extremely morphologically conserved (Kozak et al., 2006). Hypothesized to be driven by phylogenetic niche conservatism (Kozak and Wiens, 2006), this radiation has produced multiple cryptic species that have been delimited, at least in part, using molecular data (e.g., Highton, 1995, 1999; Highton and Peabody, 2000; Kuchta et al., 2018).

The identification of cryptic species is especially important for species that are experiencing rapid population decline. One such susceptible species is the green salamander (Aneides aeneus), which is listed as Near-threatened with decreasing populations by the IUCN Red List (Hammerson, 2004), Vulnerable by NatureServe (Hammerson and Dirrigl, 2017), and has been petitioned for listing under the federal Endangered Species Act. Note, however, that both the IUCN and NatureServe threat categories are in need of revision, having last been assessed 14 and 11 years ago, respectively. Consequently, these threat categories likely do not reflect the present-day population trends.

Aneides aeneus is a habitat and microhabitat specialist, making the species particularly vulnerable to habitat loss and fragmentation. These salamanders are generally associated with crevices and hollows of rock outcrop, woody, and arboreal habitats within cove forests (Gordon, 1952; Cupp, 1991; Waldron and Humphries, 2005; Smith et al., 2017). These scansorial animals are frequently found in shady crevices in rock outcrops in mixed mesophytic forests from ∼300 to 1200 m elevation. Although they have been observed to exhibit arboreal tendencies, they typically are restricted to areas harboring rock outcrops (Gordon, 1952; Waldron and Humphries, 2005; Smith et al., 2017). A large number of nesting sites have been detected on rock faces (but see Pope, 1928), which may in turn lead breeding populations to be patchily distributed (Petranka, 1998; Corser, 2001). In the 1970s, green salamander populations collapsed in the disjunct Blue Ridge Escarpment (Snyder, 1983, 1991; Corser, 2001). However, those outside of the escarpment within the western portion of their range remained stable (Snyder, 1991). Corser (2001) studied trends in several populations in the Blue Ridge Escarpment throughout the 1990s and found a 98% decline in relative abundance since 1970; extensive surveys over the last 30 years by the North Carolina Wildlife Resources Commission (NCWRC) indicate that declines appear to have continued in these populations over the last decade and a half.

In this study, we examined spatial genetic patterns in A. aeneus through the integration of population genetic and phylogenetic methods. We analyzed 12 polymorphic microsatellite loci from specimens collected from across the range of A. aeneus in western North Carolina, including samples obtained from Virginia to represent the northern Appalachian portion of their range. In doing so, we sought to define geographic patterns of genetic variation and population structuring across the range of A. aeneus in western North Carolina. We additionally conducted phylogenetic analyses using 1) targeted mitochondrial sequences (Cyt b and 12S rDNA) of 74 samples across the range of A. aeneus and 2) a nuclear multilocus SNP dataset for ten individuals across the range obtained using 3RAD-seq (Bayona-Vásquez et al., 2019: bioRxiv:205799). Population assignments obtained using microsatellites are concordant with the identification of clades in the phylogenetic component of this study. Preliminary morphological comparisons among near-topotypic A. aeneus and the clade restricted to the Hickory Nut Gorge (HNG) of western North Carolina reveals cryptic, but unequivocal differentiation that is concordant with molecular evidence. As such, we describe this lineage as Aneides caryaensis, a new species, the first cryptic species to be described within the subgenus Castaneides (eastern Aneides: Dubois and Raffaëlli, 2012).

MATERIALS AND METHODS

Population sampling.—For the microsatellite analysis, we collected tail tips (<5 mm) from each specimen and immediately preserved them in 95% ethanol. We then stored these in a –70°C freezer at Tangled Bank Conservation until extraction. We sampled in the Blue Ridge Escarpment (BRE: 197 samples) and Hickory Nut Gorge (HNG: 28 samples) areas of western North Carolina. An additional 17 samples collected from western Virginia were donated to the study by the University of Virginia's College at Wise and were stored in a like manner. In total, 242 samples were collected for microsatellite analysis.

For analysis of mtDNA, we employed an iterative strategy of geographically dispersed sampling, analysis, targeted re-sampling, and re-analysis (e.g., Jockusch and Wake, 2002; Morando et al., 2003; Moritz et al., 2009). From the late 1990s to the middle 2000s, we collected tail clippings from specimens at 72 different populations (74 total samples) and immediately placed them in 95% ethanol and maintained them at –80°C until extraction.

Microsatellite discovery.—Total genomic DNA was extracted from eight individuals using Qiagen DNeasy kits (Qiagen, Valencia, CA) according to the manufacturer's protocols. An Illumina paired-end shotgun library was then prepared and sequenced on an Illumina HiSeq. Resulting reads were analyzed using the program PAL Finder v0.02.03 (Castoe et al., 2012). This program extracted reads containing di-, tri-, tetra-, penta-, and hexanucleotide microsatellite repeats. Resultant positive reads were then analyzed with Primer3 (version 2.0.0; Koressaar and Remm, 2007; Untergasser et al., 2012) to design primer sequences. Only those loci with sequences that were identified once or twice in our total set of reads were selected for optimization in order to avoid primer sequences that may amplify in multiple locations across the genome.

Primer testing and validation.—From the set of putative microsatellite loci described above, we selected 48 loci for amplification and polymorphism across 20 individuals from three populations according to the protocol described by Eschbach and Schöning (2013). This protocol does not necessitate individual genotypes. Instead, presence/absence of alleles in pooled DNA samples are used to estimate within and among population variability in allele frequency. Loci were then scored according to the degree of this variation using the scoring procedure described by Eschbach and Schöning (2013).

Multiplex PCR reactions were conducted using Qiagen Multiplex PCR kits (Qiagen, Valencia, CA) in accordance to the manufacturer's protocols, albeit conducted in 10 µL volume reactions. Each locus was amplified individually using labeled primers. PCR cycles consisted of an initial 15 min denaturing step at 95°C, followed by 35 cycles of 30 sec denaturing at 94°C, 90 sec annealing step at 57°C, and a 60 sec extension at 72°C. A 30 min extension step at 60°C followed the 35 cycles, prior to cooling the product to 4°C. Resultant amplification products were then sized on a ABI 3730 DNA Analyzer at Florida State University using a GS500HDROX size standard. Results were then scored using Geneious version 6.1.6 (Kearse et al., 2012). Of the 48 loci screened, 12 were sufficiently polymorphic within and among populations and amplified consistently with easily scored genotypes (Table 1).

Table 1

Results of the analysis of molecular variance (AMOVA).

Microsatellite characterization.—With the 12 polymorphic microsatellite loci identified using the method of Eschbach and Schöning (2013), we then validated these loci by genotyping the 242 samples collected for this study. Procedures for amplification, genotyping, and scoring of these samples was the same as previously described. Once scored, observed (H_O) and expected (H_E) heterozygosities and number of alleles per locus (k) were calculated using GenoDive v2.0b25 (Meirmans and van Tienderen, 2004). The presence and frequency of null alleles at each locus were tested for using the R package Genepop (Rousset, 2008).

Microsatellite analysis.—To identify how populations are structured on the basis of geospatial patterns of genetic variation, we used TESS 2.3.1 (Chen et al., 2007). TESS is a program that is applied in a Bayesian framework, building Voronoi tessellations to group individuals through the identification of geographical discontinuities in allele frequencies, assuming that geographically clustered individuals will be more likely to have similar allele frequencies than more distant individuals. The number of populations (K) was estimated by running five repeat models for K = 2 through K = 6, in which 1,000,000 iterations with a burn-in of 1,000 were conducted. As suggested by François and Durand (2010), Deviance Information Criterion (DIC: Spiegelhalter et al., 2002) values were averaged between repeat runs, plotted against their respective Ks, and the value of K at which DIC plateaued was chosen as our best supported model. Populations that were defined according to this approach were then used in later phylogenetic analyses.

Once populations were defined using TESS, measures of population genetic diversity were calculated in GenoDive v2.0b25 (Meirmans and van Tienderen, 2004). Specifically, we calculated Nei's inbreeding statistic, G_IS (Nei, 1987), to quantify within-population genetic diversity. To quantify among-population diversity, we took two approaches. We performed a locus-by-locus analysis of molecular variance (Excoffier and Laval, 2005) using 10,000 permutations and a ploidy-independent mutation model. We then calculated pairwise values of F_ST using 1,000 permutations for significance, followed by a Bonferroni correction. These measures served to complement our phylogenetic reconstructions, providing an indication of the extent of gene flow occurring between recovered clades.

Targeted sequencing, alignment, and model selection.—Genomic DNA was extracted from frozen tissues using a PuregeneTM 68 kit (Gentra Systems, Inc.). Mitochondrial Cytochrome b (Cyt b) was amplified by PCR using the Cytochrome b primers MVZ 15 (5′–GAACTAATGGCCCA CACWWTACG–3′) and Cyt b2 (5′–CCCCTCAGAATGA TATTTGTCCTA–3′; modified from Moritz et al., 1992). 12S ribosomal DNA was amplified using primers 12SJ-L (5′–AAAGRTTTGGTCCTRRSCTT–3 ′) and 12SK-H (5 ′–TCCRGTAYRCTTACCDTGTTACGA–3′) as described by Goebel et al. (1999). Amplifications began with 4 min at 94°C followed by 45 cycles at 94°C for 1 min, 50–55°C for 1 min, 72°C for 1 min, followed by 72°C for 5 min. PCR products were purified using the GeneClean II kit (Bio 101). These products were then sequenced in both directions using a 377 ABI automated sequencer and the ABI BigDye^TM Terminator v3.0 Ready-Reaction Cycle Sequencing kit, but with one quarter the recommended reaction volume. Mitochondrial DNA sequences were aligned using Clustal X (Thompson et al., 1997).

Using the concatenated aligned mitochondrial sequences, we then evaluated models of nucleotide substitution using jModeltest 2.1.10 (Darriba et al., 2012). Model selection was carried out using the Bayesian information criterion (BIC: Schwarz, 1978) to evaluate the relative model fits of 88 models of evolution. A single model of nucleotide substitution was assumed for the two genes as the mitochondrion is a single non-recombining locus. Considered models included those with equal or unequal base frequencies, those that included rate variation among sites (Number Rate Categories = 6), and those with or without a proportion of invariant sites. The base tree for fitting models of nucleotide substitution was generated via Maximum Likelihood, using Subtree Pruning and Regrafting as the tree topology search operator. Cyt b and 12s rDNA sequences from GenBank of A. hardii (NC_006338) and A. flavipunctatus (NC_006327) were used as outgroups (Mueller et al., 2004).

3RAD sequencing and bioinformatics.—Individual extractions were normalized and prepared using a 3RAD library procedure (Adapterama III; Bayona-Vásquez et al., 2019: bioRxiv: 205799). The three enzymes used during the digestion step were BAMHI, MSPI, and ClaI. Each sample was then quadruple-indexed, limited-cycled in PCR, and cleaned using speed beads (Rohland and Reich, 2012) following the 3RAD procedure. Finally, samples were pooled together, size selected for 500 bp on a Pippin Prep (Sage Science), and sequenced on an Illumina HiSeq with a PE 150 kit (Illumina Inc.) with 5 million reads per sample. The 3RAD sequence data was demultiplexed, quality assessed, clustered, consensus called, and assembled de novo, using ipyrad v0.7.28 (Eaton and Overcast, 2016). The params file used for ipyrad is included in the supplement (see Data Accessibility).

Sequencing was done using the Illumina HiSeq 2500 platform at the University of Georgia Genomics and Bioinformatics Core. The resultant 30,524 SNPs were then filtered in VCFtools v0.1.14 (Danecek et al., 2011). Resultant SNPs following our filtering met the following requirements: minimum and maximum number of alleles per site of 2, minimum mean depth of coverage of 5, minor allele frequency of 0.2 (to remove singletons), present in at least 50% of samples. Indels were removed. Using these resultant SNPs, we produced a second SNP dataset by excluding heterozygous sites.

We produced multiple sequence alignments of these concatenated filtered SNPs in fasta format using PGDspider v2.1.1.2 (Lischer and Excoffier, 2011) and custom scripts. When converting from variant call to fasta format, PGDspider produces two concatenated sequences per sample: one generated from the concatenation of the first allele across all sites, and one from the second allele across all sites. In this way, we produced two alignments—one from the first allele at each site using the homozygous only SNPs, and one produced by concatenating the haploid, concatenated sequences of the first and second allele (using all sites) as described by Margres et al. (2018).

Mitochondrial phylogenetic reconstruction.—To infer our phylogeny using both our Cyt b and 12S rDNA mitochondrial sequences using BEAST 2 (Bouckaert et al., 2014), we left clock-model, site-model, and trees linked across genes as they are inherited as a single non-recombinant locus. Parameters regarding the site model (gamma shape, substitution rate, etc.) for each gene were specified according to the best-supported model from jModeltest 2; note, however, that we allowed for six gamma rate categories as the default of four may lead to excessive discretization of rate heterogeneity (Jia et al., 2014). We used a strict-clock model as our clock prior, provided the well-documented clock-like pattern of evolution of the mitochondrial genome. Provided the 0.8% pairwise sequence divergence per million years between salamandrid genera Taricha and Notopthalmus (Tan and Wake, 1995), we used a clock rate of 0.004 substitutions/ site/MY []. Division by two is necessary to convert the pairwise sequence divergence to an estimate of per-branch divergence. A Yule model with a birth rate prior of 1/X was chosen to serve as an uninformative prior (Drummond and Bouckaert, 2015). All other priors were used in their default settings. The MCMC was allowed to run for 100,000,000 generations, with a burn-in of 10,000 generations, sampling every 10,000 generations.

Resulting trace files were visualized in Tracer v1.6 (Rambaut et al., 2014) to estimate Effective Sample Sizes (ESS) and to check for convergence of the chain. Trees were then uploaded into TreeAnnotater v2.3.1. A burn-in of 10% was set for the annotation, with the target tree set as the maximum clade-credibility tree. Node heights of the resultant phylogenetic tree were set using the mean heights across the posterior distribution of trees. Annotated trees were then visualized using FigTree v1.4.2 (Rambaut et al., 2014).

Nuclear phylogenetic reconstruction.—SNPs obtained via 3RAD-seq were then used to produce two unrooted phylogenies, one using SVDquartets (Chifman and Kubatko, 2014) in PAUP* v4.0a157 (Swofford, 2003) as in Margres et al. (2018), and one using RAxML v8 (Stamatakis, 2014) on the CIPRES portal (Miller et al., 2010). SVDquartets uses SNPs to infer relationships among quartets of taxa under a coalescent model in which each site is assumed to have its own genealogy. Thus, the choice of concatenating each allele is compatible with this model. As input for SVDquartets, we used the concatenation of all alleles as input. In contrast, the concatenated sequence of the first allele at every homozygous site was used for RAxML.

Using SVDquartets, all quartets were estimated under the multispecies coalescent model (expecting matrix-rank 10), and these quartets were assembled using the QFM algorithm. Confidence in tree topology was quantified through nonparametric bootstrapping. A consensus tree was produced by summarizing across bootstrapped trees using the SumTrees program as implemented in DendroPy v4.3.0 (Sukumaran and Holder, 2010) using –force-unrooted and –min-clade-freq = 0.25. Consensus trees were visualized using FigTree v1.4.3 (Rambaut et al., 2014).

Using RAxML, phylogenies were inferred under the GTR model disabling rate heterogeneity, assuming each SNP evolves independently. Bootstrapping (100 replicates) was conducted using a rapid bootstrap analysis following a search for the best scoring ML tree (–f a option). The best tree and bootstrap nodal support values were visualized in FigTree v1.4.3 (Rambaut et al., 2014). We implemented the Felsenstein acquisition bias for invariant sites by specifying the number of sequenced sites not included in our final set of SNPs (Leaché et al., 2015). Because the exclusion of heterozygous sites led to the removal of a large number of sites informative to distinguishing individuals from the Hickory Nut Gorge, we also inferred a phylogeny without an ascertainment bias, thus including heterozygous sites. These results may be found in the supplement (see Data Accessibility).

Morphological materials collection.—Specimens used for this study included seven adults (3 males, 4 females) from the Hickory Nut Gorge (see below) and a series of nine adults of Aneides aeneus (6 males, 3 females) in the Museum of Vertebrate Zoology (MVZ) collected from SE Fort Payne along Little River Canyon, De Kalb County, Alabama (34.39°N, 85.62°W). Herein, we explicitly tested the hypothesis that specimens from the Hickory Nut Gorge differ from topotypic A. aeneus.

The following measurements were used for preliminary morphological comparison: snout to posterior angle of vent (standard length, SL), head width (HW), snout to gular fold (SG), head depth at posterior angle of jaw (HD), eyelid length (EL), eyelid width (EW), anterior rim of orbit to tip of snout (ES), horizontal eye diameter (ED), anterior rim of orbit to external naris (EN), interorbital distance between angle of eyes (intercanthal distance, IC), interorbital distance between eyelids (IO), snout to forelimb (SF), distance separating external nares (internarial distance, IN), snout projection beyond mandible (SP), shoulder width (SW, distance between axilla across dorsum), snout to anterior angle of vent (SAV), axilla–groin distance (AX), number of costal interspaces overlapped by adpressed limbs (negative when overlapped, positive when not overlapped [limb interval, LI]), forelimb length (FLL), hind limb length (HLL), hand width (HAW), foot width (FW), length of third (longest) toe (T3), and length of fifth toe (T5). Measurements were made using digital calipers. Tooth counts are based on direct counts of clearly ankylosed teeth. Institutional abbreviations follow Sabaj (2016). Color information was derived from photographs of living specimens.

Preliminary morphological comparative analyses.—Measurements of all 18 traits were compared among topotypic A. aeneus and samples from the Hickory Nut Gorge using a two-tailed Welch two-sample t-test in R version 3.4.2 (R Core Team, 2017). We omitted tail measurements from animals with regenerated tails for statistical comparisons. Additionally, to test the hypothesis that the two species are distinguishable on the basis of morphology, we leveraged both Principal Components Analysis (PCA) and Linear Discriminant Analysis (LDA). For these analyses, we excluded tail measurements, as well as counts of pre-maxillary and maxillary teeth due to the presence of missing data. We also removed costal groove counts as this character was near invariable (excluding one sample from the Hickory Nut Gorge) and thus uninformative. The remaining 14 characters were visually inspected for normality and, where necessary and possible, transformed to improve normality. Longest toe and 5^th toe were log-transformed, and adpressed limbs were converted to absolute value (all values were negative) and subsequently log-transformed. Due to our low sample size, LDA was performed as follows. First, 11 of 16 samples (approximately two-thirds) were used for training, whereas the remaining five samples were to assess prediction accuracy. We then repeated this process 10,000 times, thus producing a distribution of prediction accuracies that we subsequently used to produce an estimate of mean prediction accuracy. Due to the low sample size in these analyses, we emphasize that results are preliminary in nature.

Phylogenetic species delimitation.—We chose to supplement our phylogenetic and population genetic analyses with species delimitation using Poisson tree processes (Zhang et al., 2013). The method Poisson tree processes (PTP) is not built around the multispecies coalescent, instead modeling the speciation process as a function of the substitution process. Briefly, PTP attempts to delimit the boundaries at which the substitution process transitions between the inter- and intraspecific substitution process as a function of the branching process, namely branch lengths. Notably, this method does not necessitate time-calibrated phylogenies, thus allowing us to use it with both of our molecular phylogenies. Further, the PTP has recently been demonstrated to perform with similar accuracy as both the Generalized Mixed Yule Coalescent (GMYC) and Bayesian Phylogenetics and Phylogeography (BPandP) and is more conservative than the GMYC (Luo et al., 2018).

We applied the PTP to both our mitochondrial and nuclear phylogenies on the PTP web server ( https://species.h-its.org/). For each phylogeny, the analysis was run for a total of 500,000 generations, sampling every 500 generations for a total of 1,000 from the posterior. The first 10% of samples were discarded as burn-in. For the mitochondrial dataset, outgroups were retained as a reference for interspecific branching processes and to minimize false positives.

RESULTS

Microsatellite characterization.—We identified, validated, and optimized 12 microsatellite loci for use in conservation studies of subgenus Castaneides and other members of the genus for which microsatellites have not been developed. We found these loci to possess varying levels of diversity (Table S1; see Data Accessibility). Although the prevalence of null alleles (>10%) is not trivial (Table S1; see Data Accessibility), their presence is distributed across populations (Table S2; see Data Accessibility). However, we do not believe this to be of great concern, as null alleles have only been found to inflate, not generate, population structure or impact hypothesis testing (Chapuis and Estoup, 2007; Carlsson, 2008). Thus, the presence of null alleles in our dataset is unlikely to generate the patterns of population differentiation documented herein.

Population structuring and differentiation.—Our analysis of population structuring in TESS (Chen et al., 2007) found that K = 4 was the best-supported model according to the method of François and Durand (2010; Fig. S1; see Data Accessibility). From here on, these populations will be referred to as the Blue Ridge Escarpment (BRE) populations, Hickory Nut Gorge (HNG) population, and Virginia (VA) population. Of the samples, 197 were assigned to the two populations comprising the BRE, 28 were assigned to the HNG, and 17 were assigned to the VA population (Fig. 1). Interestingly, whereas K of 3 reveals unambiguous assignment of all individuals collected from the BRE to one population, the best fit model (K = 4) reveals apparent population structure within the BRE (Fig. 1). However, among these two populations there is appreciable admixture. Although we discovered population structuring within the BRE, individuals from VA and the HNG were always unambiguously assigned to their own populations (Fig. 1).

Fig. 1

Population structuring in Aneides aeneus (Castaneides) as identified using 12 microsatellite loci and TESS. (A) Results of running TESS assuming the number of populations (K) ranges from two through six. Each column represents an individual, and the colors in each bar indicate the population each individual is assigned to. The height and color of each bar represents the admixture proportion for that individual. Results for a single run are shown. (B) Representation of best-fit (K = 4) TESS results in geographic space. Samples are shown as pie charts, with colors corresponding to the populations identified in A, and the fractions corresponding to each sample's admixture proportions.

Results of both AMOVA (Table 1) and pairwise F_ST (Table 2) further revealed the extent of differentiation among populations identified by TESS. Analysis of molecular variance showed that whereas 71% of variance could be explained within populations, the remaining 29% could be explained by among-population differences (P = 0.001). Pairwise F_ST revealed that whereas the two BRE populations are little differentiated (F_ST = 0.076), the HNG population is highly differentiated from both BRE subpopulations (BRE-1: 0.289, BRE-2: 0.305) despite their close geographic proximity. Likewise, the HNG and VA populations are highly differentiated (F_ST = 0.282). Lastly, the VA population exhibits an intermediate degree of differentiation from the two BRE subpopulations (BRE-1: 0.139, BRE-2: 0.192). Importantly, all comparisons are highly significant (P = 0.001).

Table 2

Pairwise F_ST values for populations of Aneides aeneus (Castaneides) identified in this study. All values are statistically significant at the <0.001 level.

Calculations of the inbreeding coefficient (Table 3), G_IS, ranged greatly between populations. Notably, the first BRE population appears to be harboring an excess of heterozygosity (G_IS = –0.05, P = 0.001). In contrast, the second BRE subpopulation (G_IS = 0.103, P = 0.001) and HNG population (G_IS = 0.250, P= 0.001) are quite inbred, with the HNG faring particularly poorly. The VA population does not appear inbred (G_IS = 0.035, P = 0.191).

Table 3

Genetic diversity and inbreeding coefficient, G_IS, at each population in this study. H_O is the observed heterozygosity and H_E is the expected heterozygosity within populations. Parenthetical values are the corresponding P-values for estimates of G_IS.

Mitochondrial sequence data and model selection.—Sequencing of mitochondrial loci Cytochrome b and 12S rDNA yielded aligned sequence lengths of 433 and 415 nucleotides, respectively. For our concatenated alignment, we found that TrN+G was the best supported model according to BIC as implemented in jModeltest 2. Model parameter estimates (Gamma = 0.2960, ncat = 6, Rate AC = 1, Rate AG = 4.9971, Rate AT = 1, Rate CG = 1, Rate CT = 8.9766, and Rate GT = 1.0) were thus implemented as our site model.

Phylogenetic inference.—Our Bayesian phylogenetic reconstruction using concatenated Cyt b mtDNA and 12S rDNA sequence data recovered four well-supported, reciprocally monophyletic lineages within A. aeneus (Fig. 2). Provided our strict molecular clock rate is accurate, the node heights (coalescent times) may be interpretable in terms of millions of years. Our tree thus indicates that the tMRCA (time to most recent common ancestor) between A. aeneus and our outgroups is roughly 30 MYA. The Hickory Nut Gorge lineage was inferred to be sister to all other lineages within the A. aeneus complex, with their tMRCA being estimated at around 11 MYA. The southern Appalachian lineage last shared a common ancestor with the lineage comprising the BRE and northern Appalachian lineages approximately 9 MYA. Lastly, the BRE and northern Appalachian lineages last shared a common ancestor around 5 MYA. Our estimates of divergence times are concordant with those of Shen et al. (2016: between 27.2 and 32.3 MYA), who used several methods to estimate the divergence times between A. aeneus and A. hardii.

Fig. 2

Mitochondrial phylogenetic lineages within Aneides aeneus (Castaneides). (A) Localities of samples used in the mitochondrial phylogenetic analysis. The type locality (Nickajack Cave) is in Tennessee, at the border with Alabama; specimens from this vicinity are nested in the southern Appalachian clade. Red stars are locations used in the nuclear phylogenetic reconstruction. (B) Bayesian consensus tree of Cytochrome b and 12S rDNA sequences as produced by BEAST 2. Node labels represent posterior probabilities for the four main lineage splits within A. aeneus (Castaneides). Outgroup contains one voucher specimen each of A. hardii and A. flavipunctatus, obtained from GenBank. Numbered stars indicate clades identified by PTP to comprise unique species across the majority of the posterior distribution. Support values for clades are as follows: 1) 59%, 2) 91%, 3) 71%, 4) 93%. See Data Accessibility for tree file.

Following SNP calling in iPyrad, a total of 30,524 SNPs were obtained. Filtering in VCFtools yielded final SNP dataset containing 3,550 SNPs across the ten samples. Reduction of these sites to only those homozygous genotypes for ascertainment bias correction reduced the number of SNPs to 2,634 in RAxML.

Using SVDquartets, four well-supported groups were identified (>95% BS support: Fig. 3A) that roughly corresponded to the four clades identified in the mitochondrial phylogeny (BRE, northern Appalachians, southern Appalachians, HNG). These findings are further substantiated by our nuclear phylogeny estimated by RAxML (Fig. 3B). This phylogeny recovered three well-supported groups (BRE, southern Appalachians, and HNG), with the northern Appalachians forming a single monophyletic group despite lower support values. Notably, the HNG and southern Appalachian lineages are separated by long branches, indicating substantial genetic divergence of each clade from the rest of the A. aeneus species complex. Removal of heterozygous sites led to the inference of zero-value branch lengths among samples from the HNG, a fact that reiterated the low genetic diversity harbored within this population.

Fig. 3

Nuclear phylogenies as estimated by (A) SVDquartets and (B) RAxML. Filled circles denote nodes with support values >95%, empty circles denote nodes with support values >90%, and all other nodal support values are indicated as text. Numbered/colored circles enclosing monophyletic groups in the trees indicate species as identified by PTP when applied to the best tree identified by RAxML. Support for these species are as follows: 1) 44%, 2) 82%, 3) 50%. Low support values for the HNG are due to zero branch lengths (see Fig. S4; see Data Accessibility). Corresponding circles are similarly shown on the tree inferred by SVDquartets. However, PTP was not applied to this phylogeny as SVDquartets does not yet infer branch lengths. See Data Accessibility for tree files.

Species delimitation.—We applied the PTP to the phylogenies obtained using BEAST and RAxML. Using the mitochondrial phylogeny, the PTP identified four species with a posterior support greater than 50% (Fig. 2): southern Appalachians (59%), BRE (91%), northern Appalachians (71%), and the HNG (93%). In contrast, when applying the PTP to the nuclear tree produced by RAxML (Fig. 3), we identified three species; however, support values were low overall, likely owing to the low number of substitutions observed within species using the reduced dataset for RAxML (excluding heterozygous sites). Specifically, we recovered the following as species: southern Appalachians (82%), northern Appalachians/BRE (44%), and the HNG (50%). The low support for the HNG is almost certainly a consequence of there being zero-value branch lengths separating the two samples. This limits the ability of PTP to test hypotheses that the substitution process differs within and among species. Analysis of a phylogeny inferred by RAxML without an ascertainment bias (and thus including heterozygous sites) confirms this. Doing so, we identified three species with support values greater than 50% (Fig. S4; see Data Accessibility): southern Appalachians (97%), northern Appalachians (57%), and the HNG (95%).

Species description.—Herein we describe the Hickory Nut Gorge lineage as a distinct species on the basis of molecular and preliminary morphological differentiation. We adopt the use of subgenus Castaneides (Dubois and Raffaëlli, 2012) when referring to the species complex of green salamanders, as this etymology was initially established to refer to eastern Aneides which was thought to harbor at least three cryptic species.

Table 4

Means of untransformed morphological characters for each Aneides aeneus and A. caryaensis. Sample sizes are indicated for each species, as are the results of a two-tailed t-test comparing means among species assuming heteroscedasticity. As with our analysis using PCA and LDA, data were normalized where possible and necessary (longest toe, 5^th toe log transformed; adpressed limbs converted to absolute value and then log transformed). Units are specified below the character being measured. Standard deviations are reported alongside means. Comparisons were made among untransformed characters.

Fig. 4

Images representing Aneides aeneus (A) and A. caryaensis (B). Aneides caryaensis is characterized by smaller and less connected lichen-like patches of bright green to yellowish-green pigment.

Fig. 5

Results of Principal Components Analysis (PCA) and Linear Discriminant Analysis (LDA) using 14 morphological characters. For both analyses, data were normalized where possible and necessary (longest toe, 5^th toe log-transformed; adpressed limbs converted to absolute value and then log-transformed). (A) PCA with points and normal data ellipses colored by species. (B) Density plot of prediction accuracy by LDA across 1,000 permutations of samples for the training and prediction sets.

DISCUSSION

In this study, we integrated population genetic, phylogenetic, preliminary morphological, and species delimitation approaches to delimit cryptic species in the green salamander complex (A. aeneus). We find strong, consistent support for the recognition of the Hickory Nut Gorge lineage as a unique species. There is also strong support for the recognition of up to three other lineages (Figs. 2, 3). However, provided the undetermined boundaries between the northern and southern populations, which apparently overlap geographically, we recommend that these be recognized and managed as evolutionarily significant units (ESUs) for the purposes of status assessments and conservation actions, and that further research be conducted to analyze the validity of these lineages as full species.

Our phylogenies, both mitochondrial and nuclear, are concordant with previous research (Sessions and Kezer, 1987) that identifies distinct lineages within the A. aeneus complex. Our mitochondrial phylogenetic analysis supports the recognition of A. caryaensis as sister to all other members of this complex, sharing common ancestry with other lineages an estimated 12 MYA. That said, these estimates must be regarded as overestimates of the divergence times among these three species, as coalescent times (inferred in the present study) are necessarily greater than divergence times (Nielsen and Slatkin, 2013).

In addition to the strong topological support, a visualization of the geographical distribution of our samples revealed the allopatric distribution of each lineage (Fig. 2). Further, each lineage is geographically assorted with the exception of the contact zone of the northern and southern clades. Comparable levels of cryptic differentiation have also been recently delimited within formerly widespread species complexes in western Aneides (Reilly and Wake, 2015; Reilly et al., 2015), and, like A. aeneus, the ambiguous nature of these secondary contact zones in both eastern and western U.S mountain ranges seems to be driven by the complex geologic history of tectonic activity, which we discuss below. In addition to the evidence presented earlier, molecular species delimitation methods consistently and strongly support the recognition of A. caryaensis as a distinct species.

We find A. caryaensis to be sister to the A. aeneus complex, the product of an ancient divergence that was followed by further population subdivision in the remainder of the complex. More detailed sampling will be required to elucidate the phylogenetic history of the complex as a whole. Accordingly, in combination with our population genetic, phylogenetic, and preliminary morphological analyses, we have formally described A. caryaensis as a distinct species. Although populations from the Hickory Nut Gorge occur within ∼25 km of populations in the Blue Ridge Escarpment, there are high levels of microsatellite differentiation (Table 2), indicating minimal gene flow. This is consistent with other members of Plethodontidae found within the Hickory Nut Gorge. For example, P. amplus is endemic to the Hickory Nut Gorge and a distinct form of P. yonahlossee that displays genetic and morphological differentiation is found within the Hickory Nut Gorge (Guttman et al., 1978; Highton and Peabody, 2000; Apodaca, unpubl. data). Formal recognition of A. caryaensis is urgent considering its small geographic distribution, which is estimated to be approximately 3,625 ha (∼35 sq. km) and comprises fewer than 25 known localities. Further, the recent population declines (Corser, 2001) of the nearby BRE lineage and the intensifying human footprint within the Hickory Nut Gorge imply that A. caryaensis may be facing similar risks of decline.

Zoogeography of Castaneides.—Arachnologists (Catley, 1994) have long recognized the striking similarities between salamanders in Aneides and syntopic rock outcrop-inhabiting relict lampshade spiders (Hypochilus). These ecologically comparable species (=cryophilic syndrome) share allopatric, disjunct portions of their range in California and the southern Rockies, as well as five eastern species having closely overlapping distributions with eastern Aneides. This includes H. coylei confined to the Hickory Nut Gorge, a microendemic spider described on subtle genital characters (Huff and Coyle, 1992) and molecularly upheld by Keith and Hedin (2012) as a sister clade to the Cumberland Plateau H. thorelli. Hypochilus coylei may have shared a common Miocene ancestor with a formerly more widely and easterly distributed H. thorelli (Catley, 1994; Hedin, 2001), and similar dynamics may have accounted for the circumscription of A. caryaensis within a single valley on the very edge of the southeastern Highlands. Additionally, the previously described species P. longicrus (now P. yonahlossee) and P. amplus are restricted to the Hickory Nut Gorge and surrounding Swannanoa Mountains as well. Although it is unclear at present whether P. longicrus is a distinct species or simply a color morph of P. yonahlossee, it seems apparent that whatever led to the diversification of A. caryaensis in the Hickory Nut Gorge also affected a number of other species.

We have noted that no obvious surface barrier separates the lower elevation populations of A. caryaensis from those in BRE, occurring just 25 km to the west (Figs. 1, 2). However, three Mesozoic-aged faults (Snipes et al., 1986) lie precisely in the region between A. caryaensis and the BRE clade and H. coylei/H. pococki. These lower hills are known as the Hendersonville bulge and lie just east of the Brevard fault zone—the generalized physiographic edge of the Blue Ridge province—belonging geologically to the lower relief of Piedmont Province (Hack, 1982). Rejuvenated global and Appalachian Miocene uplift 10–20 MYA (Gallen et al., 2013; Liu, 2014) altered the climate and ecosystems towards a more grassland- dominated terrestrial ecosystem at the expense of the mixed mesophytic forests (Kürschner et al., 2008). The faults are thought to have been the loci of small differential movements, most recently during the late Miocene ∼12 MYA that affected the gross morphology of the terrain sufficiently to act as a barrier to gene flow for these rather sedentary species (Hack, 1982; Gallen et al., 2013).

It is possible that Castaneides previously had a more easterly Piedmont distribution which subsequently was lost as forests retreated to protected upland gorges after the Miocene optimum ∼15 MYA (Kürschner et al., 2008). This hypothesis could explain why these salamander and spider populations became relictualized in the outlying Hickory Nut Gorge, surrounded as they are to the west by a much more mountainous terrain housing widespread congeners but which became extensively fragmented in the southern uplands. This scenario predicts the existence of another forested refuge on the southern Cumberland Plateau, and it was from here that the remainder of the Castaneides were subsequently founded as forest cover waxed and waned since the late Miocene. Martin et al.'s (2016) ancestral area analysis also hinted that the ancestor of Castaneides might have been a more widespread lowland and potentially arboreal form.

Karyotypic variation.—Aneides differs from other plethodontid genera in the extent of chromosomal variation (Sessions and Kezer, 1987). With few exceptions, plethodontids have 13 or 14 pairs of biarmed chromosomes and closely resemble each other. While all species of Aneides have 14 pairs of chromosomes, at least eight distinct karyotypes are found, three of them in the A. aeneus complex. One of these is “Morescalchi's aeneus III” (Sessions and Kezer, 1987), potentially important taxonomically but unavailable because of the lack of detailed geographic information (obtained by long-deceased Morescalchi from an unnamed animal dealer and suspected of being from the vicinity of Cumberland County, Tennessee). Morescalchi's materials, stored in Italy, were destroyed in the Naples earthquake in 1980. This region lies between sites where “aeneus I” (to the north) and “aeneus II” (to the southwest) have been recorded. In their map (fig. 13), Sessions and Kezer (1987) show aeneus I and II in sympatry (p. 24: “at one locality, near Chattanooga, Tennessee, a single specimen with the aeneus I karyotype has been detected among mainly aeneus II individuals”). A specimen from Hickory Nut Gorge has the aeneus I karyotype.

The report of Sessions and Kezer (1987) suggests that three karyotypes are found in eastern Tennessee, two of them in sympatry. In a recent communication to DBW (29 May 2018), Sessions provided additional details and corrected errors. For the published study, he examined three specimens of A. aeneus from “near Chattanooga”, which he now specifies as from 4 miles S Sherwood at Buck Creek Cove, Franklin Co., Tennessee (collected by Wayne Van Devender). This site is just a few km west of Nickajack Cave (now mainly inundated by a reservoir), the type locality of A. aeneus. These specimens were sent to Kezer in Eugene, Oregon. Sessions took one to Berkeley, California, for his studies, and Kezer examined the other two. In a letter to Sessions (dated 29 November 1981), Kezer reports that the specimens he examined had completely bi-armed chromosomes (i.e., aeneus II). Sessions double-checked his lab notes from the time and finds that his single specimen also was aeneus II. He cannot explain how an error was made in the published version, but he concludes: “there was no evidence of co-occurrence of aeneus I and II at that site”. Sessions sent DBW a complete list of specimens examined. There were no other specimens from Tennessee, and only the single specimen from Hickory Nut Gorge from North Carolina. They did have specimens from two sites in Alabama, both aeneus II. The map in their publication (fig. 13) is accordingly flawed and must be used with extreme care. Importantly for our current work, the specimen from Hickory Nut Gorge is aeneus I, which is found only at this site and from Kentucky and Virginia northwards to Maryland and Ohio, and nowhere else in any species of Aneides or any other plethodontid. The finding of the aeneus I karyotype in Hickory Nut Gorge is important in showing a sharp difference between this region and the region of the type locality of A. aeneus, very near the Tennessee–Alabama border.

A cladogram of chromosomal characters in Aneides shows the clade differing from the outgroup used (Plethodon) in two shared derived traits (Session and Kezer, 1987: fig. 15). Aneides aeneus differs further from other members of the genus in two more traits, and aeneus I and II are sister taxa. Karyotypic differences reported by Sessions and Kezer in the A. ferreus complex later proved to be species specific: ferreus I = Aneides ferreus, and ferreus II = Aneides vagrans (Jackman, 1998).

Conservation implications.—Aneides caryaensis faces pressing conservation concerns and need for management. This species has an extremely narrow range and is reported from fewer than 25 localities, despite extensive searches. The status of these populations is unknown, yet based on the personal observations of several of the authors, densities appear to be very low at all of the known localities, including those that had higher surveyed (by the North Carolina Wildlife Resource Commission) densities, often with only one to three individuals observed during a site visit. Additionally, these field observations are supported by our estimates of genetic diversity (Table 3), which indicate that populations of A. caryaensis are experiencing a high amount of inbreeding, a clear signature of small population sizes (Frankham et al., 2010). High levels of inbreeding can have detrimental effects on populations, including a decrease in fecundity that can in turn lead to lower population sizes and increased inbreeding, a feedback cycle often referred to as “the extinction vortex” (Frankham et al., 2010). Habitat loss and fragmentation resulting from the rapid rate of tourism, real estate development, and transportation and energy infrastructure within the Hickory Nut Gorge is likely the largest direct threat to the long-term survival of A. caryaensis. Habitat fragmentation and loss can easily occur for this species, as habitat features preferred by A. caryaensis tend to be distributed in patches across the landscape and are susceptible to disturbance from anthropogenic activities.

Similar to A. caryaensis, the BRE populations of A. aeneus also face a myriad of threats and conservation concerns. Although some populations in the BRE appeared to be recovering from the dramatic (98%) declines of the 1980s and 90s (Corser, 2001), they have been experiencing declines since at least 2005 (Apodaca and Williams, unpubl. data). The low present-day diversity observed for the BRE populations included in our study (Table 3) seems to reflect these declines. Importantly, the subpopulation structure within the BRE and apparent admixture between subpopulations attests to this complex demographic history. Although we have not explicitly tested demographic hypotheses for the BRE populations, our results imply either 1) secondary contact among previously isolated populations or 2) ongoing divergence. Future work studying the BRE populations would benefit from a thorough investigation of their demographic history.

Regardless, the low present day diversity of these populations in conjunction with their historical and continued declines emphasizes a great need to conserve and manage current populations and to develop novel approaches to mitigate threats. Necessary conservation actions will involve 1) protecting areas of known and potential habitat from certain types of logging, development, and collecting, and 2) maintaining a permeable matrix of habitat between such areas to facilitate dispersal.

In accordance with the mitochondrial and nuclear distinctiveness of this lineage, we propose the recognition of the BRE lineage of the A. aeneus complex as an ESU. This proposal is compatible with the criterion for recognizing ESUs outlined by Moritz (1994). In his article, Moritz (1994) defined an ESU as being “reciprocally monophyletic for mtDNA alleles and show significant divergence of alleles at nuclear loci.” We have demonstrated this for the BRE lineage of A. aeneus despite apparent population structuring. Although our OTU sampling for the SNP dataset is far less than that of our mitochondrial dataset, we have recovered reciprocal monophyly for the BRE populations. Further, the BRE lineage, as identified herein, is compatible with a number of other definitions of ESUs as outlined by Funk et al. (2012; e.g., Avise, 1994; Vogler and DeSalle, 1994; Fraser and Bernatchez, 2001). As such, we anticipate our recommendation for the recognition of the BRE populations of A. aeneus as an ESU to be uncontentious. Recognition of the BRE lineage as an ESU will facilitate its conservation, as ESUs are granted legal protection under the Endangered Species Act. We do note, however, that species delimitation methods included BRE populations with northern Appalachian lineages (Fig. 3). Interestingly, the PTP designated the southern Appalachian lineage as a distinct species. This is coincident with long branch-lengths as estimated by RAxML. Too little is known about these populations at present to designate this lineage as an ESU, but we anticipate that additional investigation of this lineage will uncover previously unrecognized differentiation.

Conclusions.—Herein, we have used multiple lines of evidence to identify four lineages within Castaneides. The newly described North Carolina endemic, A. caryaensis, is restricted to an extremely narrow geographic range and is in need of immediate protection. However, little is known about the ecology of this species, as most work done on Castaneides has been conducted on more westerly populations and the Blue Ridge Escarpment lineage (e.g., Gordon, 1952; Snyder, 1983, 1991; Cupp, 1991; Corser, 2001; Waldron and Humphries, 2005; Smith et al., 2017). Perhaps the most important research needs at present are population surveys, habitat management guidelines, landscape genetic analyses, and habitat use studies. We emphasize, however, that such research needs are not restricted to A. caryaensis. The BRE ESU of A. aeneus is also of great conservation concern and faces similar threats as A. caryaensis. In summary, we hope that in describing A. caryaensis, more work will be conducted to more rigorously describe the cryptic diversity harbored within Castaneides.

DATA ACCESSIBILITY

Supplemental information is available at  https://www.copeiajournal.org/ch-18-052.