The complex geological history of the Appalachian Mountains in the southeastern U.S. has contributed to a number of phylogeographic breaks across different taxa (Kozak et al., 2006; Soltis et al., 2006; Lemmon et al., 2007a; Rissler and Smith, 2010; Ennen et al., 2017; Barrow et al., 2018). Several of these divergences occurred during the Pliocene and Pleistocene, with evidence pointing to the dynamism of glacial and interglacial episodes as the main cause of diversification within taxa (Williams et al., 2004; Soltis et al., 2006). One important consequence of glacial/interglacial periods is that sea levels and rivers undergo dramatic changes, with the potential to impose barriers to gene flow among populations (Jones et al., 2006; Lemmon et al., 2007a). Changes in river courses of the Appalachian Mountains have been proposed as an important driver of population divergence in many amphibian species from this region (Jones et al., 2006; Rissler and Smith, 2010; Barrow et al., 2018), which could be attributed to limited dispersal capabilities in some taxa (Pauly et al., 2007; Degner et al., 2010). Even though amphibian fauna of the Appalachian Mountains is relatively well studied (Soltis et al., 2006; Rissler and Smith, 2010; Barrow et al., 2018), divergence within species with reduced activity periods and restricted geographic distributions requires additional study.

The Mountain Chorus frog (Pseudacris brachyphona), first described by Cope (1889), is one example of an elusive and understudied species that ranges throughout the Appalachian Mountains as far south as Georgia and Mississippi in the United States (Hoffman, 1980 and references therein; Powell et al., 2016). This species can be difficult to study due to their patchy distribution, apparently small population sizes in many areas (Green, 1952; McClure, 1996; Mitchell and Pauley, 2005), and crypsis in their natural habitat (Schwartz, 1955). In addition, P. brachyphona is challenging to detect outside the brief breeding season, which occurs in late winter and early spring, with the peak of their season depending upon rainfall and temperature in the local area (Schwartz, 1955; Green, 1964). Breeding activity typically lasts for a few weeks (Forester et al., 2003), and the shortness of the season during which these frogs can be reliably detected has probably contributed to the paucity of information available for this species. Pseudacris brachyphona breed in roadside ditches, pools of standing water in roads, low wet areas on mountainsides, and in formerly forested areas of clear-cut fields and slopes (Green, 1938; Hoffman, 1980). Outside of the breeding season, they disperse to nearby upland habitats (Mitchell and Pauley, 2005).

Whereas some life history information and morphological data have been published for P. brachyphona (Walker, 1932; Barbour and Walters, 1941; Wright and Wright, 1949; Schwartz, 1955; Barbour, 1957), molecular data have only recently been applied to understand the evolutionary history of this species. Previous genetic studies included P. brachyphona in genus-wide phylogenetic analyses of Pseudacris based on albumin and allozyme data (Maxson and Wilson, 1975; Hedges, 1986), mitochondrial sequence data (Moriarty and Cannatella, 2004; Lemmon et al., 2007a, 2007b), and nuclear sequence data (Barrow et al., 2014; Banker et al., 2020). These studies determined that P. brachyphona is the sister species to Brimley's Chorus Frog (P. brimleyi) from along the Atlantic Coast, both of which have fast pulse rate calls. The clade formed by these two species is sister to the remaining taxa in the trilling frog clade, which is a subgroup of taxa with pulsed calls within Pseudacris (Lemmon et al., 2007a).

Although the phylogenetic position of P. brachyphona has been resolved, information is lacking regarding genetic structure within the species. The only intraspecific genetic data available to date are from a 25-individual mitochondrial sequence data set (Lemmon et al., 2007b). This study found a strongly supported (bootstrap = 100) geographic split between samples referred to as the Northern and Southern clades. The Northern clade included samples from Kentucky, Ohio, Tennessee, Virginia, West Virginia, and Northern Alabama, whereas the Southern clade included those from Alabama, Georgia, and North Carolina. Furthermore, Lemmon et al. (2007b) observed a mismatch between morphological/acoustic identification of some individuals and their mitochondrial haplotype, suggesting that P. brachyphona may hybridize with another chorus frog, P. feriarum, with which it overlaps in lower elevation portions of its range.

Despite evidence of a genetic split within P. brachyphona, Lemmon et al. (2007b) refrained from making taxonomic changes until additional lines of evidence regarding the distinctiveness of the Northern and Southern clades were assembled. Specifically, their study proposed examining ecological and behavioral data, as well as determining whether nuclear sequence data would be concordant with their mitochondrial data. To continue the work begun by Lemmon et al. (2007b), this study aims to: (1) investigate whether nuclear sequence data from >1800 loci concur with mitochondrial data regarding a deep genetic split between the Northern and Southern clades described by Lemmon et al. (2007b); (2) test for genetic evidence of hybridization between P. brachyphona and closely related congener P. feriarum; (3) assess whether clades of P. brachyphona occupy different ecological niches; (4) test for significant behavioral differences between clades with respect to male acoustic signals; and (5) test for morphometric differences between clades. We evaluate this information to define whether the clades within P. brachyphona require taxonomic changes and to infer possible geographical features that contributed to this phylogeographic break.

MATERIALS AND METHODS

Sample selection and DNA extraction.—Samples of P. brachyphona and outgroups were collected across their range for genetic analysis in the early 2000s. Assistance in these collections was provided by a network of herpetologists across multiple states (see Acknowledgments). A total of 36 chorus frog specimens (33 P. brachyphona and outgroups: two P. feriarum and one P. brimleyi) were collected from 29 counties (33 localities) across the eastern United States (Fig. 1, Table 1). Appropriate state scientific collecting permits and IACUC protocol approvals were obtained prior to specimen collection. Frogs were either toe-clipped or dissected for liver, leg muscle, and heart tissue. Tissues were frozen in liquid nitrogen or preserved in tissue buffer or 95% ethanol and stored at –80°C. Voucher specimens were deposited into either the Texas Natural History Collection (University of Texas at Austin Biodiversity Center), the Fort Hays State University's Sternberg Museum of Natural History, or the University of Florida Museum of Natural History. Genomic DNA was extracted from tissue samples using the OMEGA bio-tek E.Z.N.A. Tissue DNA kit.

Fig. 1

Samples of the Mountain Chorus Frog collected for this study and used in genetic analyses. Frogs were collected at several locations (solid dots) throughout the Appalachian Mountains in Eastern U.S. As outgroups, specimens of P. brimleyi and P. feriarum (solid and hollow triangles, respectively) were collected.

Table 1

Genetic sampling of Pseudacris across eastern United States. ID number refers to the ID of the genetic sample processed at Florida State University's Center for Anchored Phylogenomics. An “n/a” refers to samples without a voucher (tissue sample only). Field collection abbreviations are as follows: ECM–Emily Moriarty Lemmon; JA–J.J. Apodaca; JTC–Joseph T. Collins; MHP–Museum of the High Plains (Fort Hays Sternberg Museum of Natural History).

Data generation and assembly.—Samples were processed and analyzed at the Center for Anchored Phylogenomics at Florida State University ( https://anchoredphylogeny.com) following the methods described in Lemmon et al. (2012) and Prum et al. (2015). Genomic DNA was quantified via Qubit^™ fluorometer using a dsDNA HS Assay kit (Invitrogen^™). A 2% agarose gel was used to visualize the size and quality of the DNA. This information was used to determine the amount of time needed to shear the sample to size range of 175–575 bp on a Covaris E220 focused ultrasonicator (Covaris Inc.). Sheared DNA samples underwent genomic library preparation following the protocol in (Meyer and Kircher, 2010) in 96-well plates using a Beckman-Coulter Biomek FXp liquid handling robot.

Indexed samples were pooled (18 samples per pool) in equal concentrations such that each pool contained a total of 750 ng of library DNA in 3.4 µL of water. Hybrid enrichment was conducted on each library pool using the custom Pseudacris v.1 probes (Banker et al., 2020), produced in a SureSelectXT kit (Agilent Technologies). Enriched pools were quantified and pooled together in equimolar concentrations for sequencing. The sequence pools were size-selected further prior to sequencing using the electrophoresis-based PippinHT system (Sage Science, Inc.) for the size range 400–600 bp (including library adapters) in order to decrease sequencing bias towards short fragments containing mostly adapters. All sequencing pools were sequenced on an Illumina HiSeq 2500 lane in the Translational Science Lab, Florida State University's College of Medicine. The sequencing protocol included 200 cycles for each of the paired reads, in addition to the dual 8 bp indexing reads.

After low quality reads were removed using the CASAVA high chastity filter, reads were demultiplexed using the 8 bp indexing read with no mismatches tolerated. Adapter sequences were subsequently removed, and overlapping reads were merged using the Merge program described by Rokyta et al. (2012). Probe regions and flanks were recovered for each sample using the quasi-de novo assembly approach described in Prum et al. (2015) and Hamilton et al. (2016). Probe region sequences from both P. feriarum and P. nigrita were used as references in the assembly. To reduce the effects of rare indexing errors and cross contamination, we discarded assembled contigs derived from less than 10% of the typical number for that locus (compared to other homologs for all individuals). Orthology was assessed for each locus using pairwise divergence across the homologs for each locus (see Hamilton et al., 2016 for details). Alignments were constructed using MAFFT v7.023b (Katoh and Standley, 2013) for each set of orthologs and subsequently masked and trimmed following Hamilton et al. (2016), but with the following parameters: MINGOODSITES = 14, MINPROPSAME = 0.5, and MISSINGALLOWED = 18. These alignments were manually inspected in Geneious R9 (Biomatters Ltd.; Kearse et al., 2012) to ensure that misaligned regions were not present.

Phylogenetic inference.—Gene trees for each sequence alignment were generated by using RAxML v8.2 (Stamatakis, 2014) assuming GTR + Γ as substitution model and including P. feriarum and P. brimleyi as outgroups. We performed 100 replicates for each gene tree in order to get bootstrap support and to estimate a species tree using ASTRAL-III v5.6 (Zhang et al., 2018). Given that contracting low-support branches may improve overall inference when using coalescent methods (Zhang et al., 2017), we contracted branches with bootstrap support <10% in each of the gene trees. The species tree was then estimated by performing multi-locus bootstrapping over the 100 replicates previously generated for each gene tree. In addition to the ASTRAL-III species tree, we estimated a phylogeny in RAxML v8.2, this time using a concatenated version of the 1,669 loci alignments, assuming GTR + Γ as substitution model, and with bootstrap support based on 100 replicates. Results were visualized in TreeGraph v2.14 (Stöver and Müller, 2010).

Bayes factor species delimitation.—To test for the existence of independent lineages within P. brachyphona, we calculated Bayes delimitation factors (BF) for three different models within the phylogeny of P. brachyphona (Leaché et al., 2014). The first model lumped all clades of P. brachyphona in one single clade, which is concordant with current taxonomy (single species). The second model split the phylogeny of P. brachyphona in two clades as supported by geographic subdivision (see Results). Finally, we tested a null model in which individuals were randomly assigned to the two clades from model 2. Marginal likelihood estimates (MLE) were obtained for each model using 48 sampler path steps of 100,000 MCMC repetitions. After a 10% burn-in, the highest MLE from each model was used to calculate Bayes factors as BF = 2*(Model Highest MLE–Alternative Model MLE). The BF analysis was performed in SNAPP v1.5 (Bryant et al., 2012) as implemented in BEAST v2.6 (Bouckaert et al., 2019). Parameter files for SNAPP inference are available as supplemental material (see Data Accessibility).

Genetic clustering and population structure.—We assigned each P. brachyphona to genetic clusters with Structure v2.3 (Pritchard et al., 2000) based on a 21,702 SNP data set, as well as a subset of this data set obtained by extracting one SNP per locus (see Data Accessibility to obtain these data sets and Structure parameter file). For these analyses, we included the two samples of P. feriarum as outgroups: one originating from outside the geographic distribution of P. brachyphona (Colleton County, South Carolina), and one from the southernmost portion of the distribution of P. brachyphona (Macon County, Alabama). In both data sets, we tested for the most likely number of populations from K = 1 to K = 5 using 100,000 steps for parameter estimation, and an additional 100,000 steps as burn-in. We assumed correlated allele frequencies and no linkage among loci. To improve accuracy in estimation of K, we adjusted the ALPHA parameter to 0.5 (Wang, 2017). Runs were automated with structure_threader (Pina-Martins et al., 2017), performing 20 runs per each K and selection of best Ks based on the Evanno Test (Evanno et al., 2005).

To provide additional support for our Structure analysis, we performed a discriminant analysis of principal components (DAPC) to our SNP data set as implemented in the R package adegenet (Jombart, 2008; R Core Team, 2018). First, the most likely number of clusters was assessed using the find.clusters function to estimate the Bayesian information criterion (BIC) for each K from K = 1 to K = 10. With the most likely K, we performed cross-validation (xvalDapc function) to identify the highest number of informative principal components (PCs) while avoiding overfitting of the cluster membership predictions. Finally, discriminant functions were estimated with those retained PCs. Possible hybrid individuals (P. brachyphona x feriarum) were identified by estimating their hybrid index (h) in GenoDive (Meirmans and Van Tienderen, 2004).

Tests of isolation-by-distance.—We tested for the correlation between geographic distance and genetic distance using our SNP data set. A genetic distance matrix was calculated for all P. brachyphona using the function dist.genpop in adegenet (Jombart, 2008) after removal of hybrid individuals. The geographic distance matrix was calculated using the “Distance matrix” tool in QGIS v3.6 (QGIS Development Team, 2017). Estimation of correlation between the matrices was made in R using the ade4 package (Dray and Dufour, 2007) with its function mantel.randtest. In addition, we tested for isolation-by-distance within the Northern and Southern clades of P. brachyphona.

Morphometric analyses.—We measured 191 P. brachyphona kept preserved at five museum collections (Supplemental Table 1; see Data Accessibility). The split between Northern and Southern clades occurs at approximately the Tennessee River and southern portion of the Blue Ridge Mountains (see Results). We used these geographic features to assign individuals from 25 counties and 6 states to the Northern and Southern clades. Measurements were taken by a single person and using an Insize digital caliper with resolution of 0.01 mm attached to its automatic data entry device. The measurements taken are defined in Duellman (2001): snout–vent length (SVL), head width (HW), head length (HL), tympanum diameter (TD), eye width (EW), snout length (SL), femur length (FL), tibia length (TL), foot length (FoL), and snout angle (SA; [arcsine((head width/2)/head length)] x 2). To identify overall differences in size of these features, we performed a principal component analysis (PCA) based on the regression residuals of each variable against SVL to control for body size effects. With the same residuals, we performed a series of permutation tests (1,000 replicates) for the mean residual difference between Northern and Southern clades individuals. All analyses and graphs were made in R and its packages ggplot2, ggfortify, and dplyr (Wickham et al., 2015; Tang et al., 2016; Wickham, 2016).

Species distribution models and tests for niche overlap.—To test for ecological divergence within P. brachyphona, we constructed species distribution models (SDM) that considered all known distribution of P. brachyphona, only the Northern clade, or only the Southern clade. In addition to the associated geographic coordinates from our specimens, we collected occurrence records from the Global Biodiversity Information Facility database (GBIF, 2019) that had a precision of at least five decimals and relied on specimens stored in museum collections. The complete list of specimens can be found in the supplemental material (Supplemental Table 2; see Data Accessibility).

We collected the 19 bioclimatic variables as provided with a resolution of 30 s (∼1 Km²). These variables are derived from 30 years of temperature and precipitation measurements (Fick and Hijmans, 2017). We obtained elevation data from DIVA-GIS ( https://www.diva-gis.org/gdata), which is provided in the form of country-level 30 s resolution rasters from the SRTM30 data set. Spatial data for vegetation was collected from LANDFIRE in the form of 30 m resolution Existing Vegetation Type (EVT) rasters (U.S. Department of Agriculture and Department of the Interior, 2014). All spatial layers (Bioclim, elevation, and EVT) were clipped to the same extent, aligned, and converted to ASCII format in QGIS v3.6.

To avoid overfitting of species distribution models, we tested for correlations among the environmental variables using the layerStats function from the raster package in R (Hijmans and van Etten, 2017). Variables that were highly correlated (r < –0.7 or r > 0.7) with two or more variables were discarded from subsequent analyses. We excluded EVT from this correlation analysis due to the categorical nature of those data. This resulted in a set of ten variables that was used for the generation of all species distribution models.

We used the maximum entropy algorithm implemented in MaxEnt v3.4 (Phillips and Dudík, 2008) to generate species distribution models. To that end, we produced three different partitions of the occurrence data and produced SDMs with them: (1) all the occurrences of P. brachyphona, (2) only the Northern clade occurrences, and (3) only the Southern population occurrences. We ran ten replicate models and assessed the model's power by means of 10-fold cross-validations and calculation of the Akaike information criterion (AIC). Following Morales et al. (2017), we allowed environmental variable responses to be modeled after quadratic or threshold features, and with output probability of occurrence following a logistic function.

Asymmetric background tests were performed to ascertain whether or not niche overlap occurs between the two clades' SDMs. In these background tests, Schoener's D and Warren's I niche overlap indexes for the Northern clade SDM were compared to a distribution of 100 D and I values from SDMs generated from randomly selected points from the Southern clade distribution range. The 100 replicate MaxEnt SDMs were run in ENMtools and rmaxent (Baumgartner et al., 2017; Warren et al., 2019) as implemented in R. We also tested for niche overlap in the opposite direction, namely, D and I values of the Southern clade SDM vs. the random distribution from the Northern clade.

Divergence in acoustic signals within P. brachyphona.—Acoustic signal data were recorded from two locations in Tennessee and Alabama corresponding to the Northern and Southern clades, respectively (Supplemental Table 3; see Data Accessibility). No recordings were analyzed from locations where brachyphona x feriarum hybrids were likely to occur. A Sennheiser ME67 shotgun microphone was used to record calls digitally onto a Fostex FR2LE portable recorder with a SuperMod upgrade (Oade Brothers Audio, Thomasville, GA), using a sampling rate of 44,100 Hz with a sample size of 16 bits. Audio files of frog calls in .wav format were analyzed using SoundRuler (reviewed by Bee, 2004;  http://soundruler.sourceforge.net/) using the following settings: Spectrogram FFT length 1024, Hanning window size 1024, amount of overlap between FFT samples 900, and power spectrum FFT length 2048. Call characters were either extracted from the raw data output or calculated from the raw data following Lemmon (2009). Call characters examined included pulse rate (PR; 1/time from 10% maximum amplitude [pulse onset] to 10% maximum amplitude [onset] of next pulse within call), pulse number (PN; number of pulses), and call dominant frequency peak (DFP; call dominant frequency at the call maximum amplitude). Pulse rate was significantly correlated with temperature and was thus corrected to a common temperature of 14°C using a species-specific slope from linear regression analyses following the equation: ((14°C – recording temperature) * slope) + pulse rate, where the slope = 6.47. The slope of the linear regression was estimated from a broader set of recordings beyond what is presented in this paper (EML, unpubl. data).

To compare acoustic differences between the Northern and Southern clades of P. brachyphona, we conducted randomization tests and Welch t-tests between the two groups. We randomized the individual mean values of PR, PN, and DFP for 1,000 times, and each time calculated the mean difference between the Northern and Southern clades of each parameter. To perform a two-tailed test, we counted the number of times the empirical mean difference was greater than or equal to the randomized mean difference and estimated its probability of occurrence (P-value). The t-tests were also performed for each individual acoustic trait. All analyses were carried out in R.

Institutional abbreviations follow Sabaj (2020).

RESULTS

Phylogenetic relationships within P. brachyphona.—The species tree and concatenated ML trees showed evidence of two main clades within P. brachyphona (Fig. 2). These clades are geographically separated and correspond to the Northern and Southern clades of P. brachyphona described by Lemmon et al. (2007b). Local posterior probabilities in ASTRAL and bootstrap support from RAxML were 100% for each clade. Within each clade, support for internal nodes varied, but when ASTRAL showed local posterior probabilities of >95%, RAxML also showed similarly high bootstrap support. ASTRAL and RAxML newick files are available as supplemental material (see Data Accessibility). Bayes factor analysis showed strong support towards the two-independent lineages model (MLE = –121446.29) when compared to the single-lineage model (MLE = –135936.16). The comparisons against the model with the highest MLE (two-independent lineages) yielded estimated BF = 28979.74 for the single-lineage model, and BF = 29022.24 for the randomized null model (Table 2).

Fig. 2

Species tree for P. brachyphona as estimated in ASTRAL. Numbers on the branches are local posterior probabilities/boostrap support values from RAxML (A). The identifiers are field numbers, followed by the “I” number (DNA sample identifier, see Table 1). State abbreviation and county where the sample was collected are also shown. The tree shows two well-supported clades, corresponding to P. brachyphona from the Northern and Southern distribution ranges, respectively. Pseudacris brimleyi and P. feriarum were included as outgroups. Branch lengths (B) are shown in the inset tree generated with RAxML.

Table 2

Results of Bayes factor delimitation (BFD) based on 48 sampler path steps with 100,000 MCMC iterations each. Marginal likelihood estimates (MLE) were used to rank models and to compute Bayes factors (BF) by comparing each of the models against the model with the highest marginal likelihood (two species). The large BF estimates give support to the two-species model over the alternative models.

Genetic clusters and population structure.—Based on the 21K SNP data set, the highest ΔK among the five tested K values corresponded to K = 2 (ΔK = 3489.9, lnL = –371679.7, SD = 36.0, Table 3), with P. feriarum (outgroups) being clustered together with P. brachyphona originating from the species' northern distribution range. Following Prunier and Holsinger (2010), we selected K = 3 as the most biologically relevant choice given that its ΔK was the second highest (ΔK = 991.1) and the log likelihood reaches a plateau at that point (lnL = –298822.0, SD = 73.6, Supplemental Fig. 1; see Data Accessibility). With K = 3, we observed Southern clade P. brachyphona in the first cluster, Northern clade P. brachyphona and a putative hybrid P. feriarum x brachyphona in the second cluster (I27685), and pure P. feriarum (outgroups) in the third cluster (Fig. 3A, Supplemental Fig. 2; see Data Accessibility). The hybrid index for the introgressed individual (I27685) was h = 0.66 (95% CI: 0.654–0.672), which provided further evidence for its hybrid identity. The selection of K = 3 was also in agreement with our phylogenetic analyses, with clusters representing the Northern and Southern clades. The results obtained by Structure from the data set containing only one SNP per locus were consistent with results from the 21K SNP data set above (Supplemental Fig. 2; see Data Accessibility).

Table 3

Likelihood estimates for each possible K during the Structure analysis and their Delta K (ΔK) for Evanno Method assessment (Evanno et al., 2005).

Fig. 3

Genetic clustering based on 21,702 SNPs from AHE loci for P. brachyphona. (A) The Structure plot shows two genetic clusters: Northern (n = 17, green) and Southern (n = 16, blue) as observed in the map (B). The pie charts correspond to the proportion of admixture as estimated in Structure. A blue arrow marks the samples from Hale County, Alabama. Two individuals assigned to the Northern cluster showed introgression from P. feriarum (“out,” gray). (C) The discriminant analysis of principal components (DAPC) required one discriminant function to explain 67.4% of the genetic variance. Lines at the bottom of the distributions represent one individual. The DAPC also shows two clusters in agreement with the Northern and Southern clusters from Structure. The location of the Tennessee River is shown.

DAPC also showed support for two genetically differentiated clusters within P. brachyphona (Northern and Southern; Fig. 3B, C). One discriminant function was sufficient to account for 67.4% of the genetic variance. The distribution of individuals along this discriminant function showed two narrow peaks with low differentiation within each cluster (Fig. 3C). Differentiation between the two clusters was high, with FST = 0.57 (95% CI: 0.567–0.580). Using BIC to test for further differentiation within the Northern cluster also suggested the existence of a single group (K = 1, BIC = 129.6 vs. K = 2, BIC = 129.9). Similarly, within the Southern cluster, the lowest BIC was observed for K = 1 (K = 1, BIC = 117.5 vs. K = 2, BIC = 118.0). Hence, no additional groups were defined within the Northern and Southern clusters (Supplemental Fig. 3; see Data Accessibility).

Tests of isolation-by-distance.—The Mantel test for correlation between genetic and geographic distances showed a pattern of isolation-by-distance among our samples (r = 0.76, P = 0.001). Nevertheless, inspection of pairwise comparisons of genetic and geographic distances yielded two separate point clouds (Fig. 4A), suggesting more genetic differentiation than expected solely under isolation-by-distance. A Mantel test for the Northern DAPC cluster yielded evidence of isolation-by-distance (r = 0.58, P = 0.002; Fig. 4B). Similarly, the Mantel test for the Southern DAPC cluster supported isolation-by-distance to explain divergence among individuals (r = 0.69, P = 0.001; Fig. 4B).

Fig. 4

Matrix correlation between pairwise genetic (Nei's) and geographic distances among individuals of P. brachyphona, after exclusion of hybrid individuals. (A) Mantel correlation test for all P. brachyphona showed significant correlation (r). The two clusters of points represent comparisons within and between genetic clusters (Northern and Southern). (B) Significant correlations were also observed for each cluster separately (Northern = black, Southern = gray).

Morphological variation between the clades of P. brachyphona.—The PCA of morphometric variables showed almost complete overlap among individuals from the Northern and Southern clades (Fig. 5). Subtle morphological differences were observable in PC1 (explained variance ∼26.8%) due to mostly to FoL, TL, HL, and FL eigenvectors, with significant differences only observed in FoL (average FoL North = 13.64 mm vs. South = 12.96 mm, P = 0.001) and HL (average HL North = 8.40 mm vs. South = 8.42 mm, P = 0.009; Supplemental Table 4; see Data Accessibility). Permutation tests also showed significant differences for HW (average HW North = 9.60 mm vs. South = 9.55 mm, P = 0.017) and EW (North = 3.05 mm vs. South = 3.12, P = 0.001; Supplemental Table 4; see Data Accessibility). It should be noted that well-known changes in size and shape during fluid preservation of voucher specimens (Simmons, 2014) could introduce morphological variation within and/or between species, thus either enhancing or obscuring real biological differences.

Fig. 5

Principal component analysis plot based on residuals from regressions of each morphometric variable against snout–vent length. Abbreviations for the labels on the PC loadings (arrows) are the same as in the text (see Materials and Methods). Convex hulls are shown for each clade.

Species and population-specific distribution models.—After adding records from GBIF, our data set included 228 occurrences of P. brachyphona. Species distribution models fitted with MaxEnt showed a high elevation distribution along the Appalachian Mountains in the eastern U.S., ranging from southwestern Pennsylvania to central Alabama (Fig. 6A). Elevation contributed 19.8% to the combined species model's predictive power, followed by temperature during the driest trimester (∼September–November) contributing 17.2%. Precipitation in the warmest trimester (∼June–August) and temperature of wettest trimester (∼April–July) contributed 15.6% and 10.9% respectively. Vegetation type contributed 16.8% to the prediction (Table 4).

Fig. 6

Maximum entropy species distribution models (SDMs) for all P. brachyphona (A), Northern (B), and Southern clades (C). Models represent average probability of occurrence based on ten replicates. We applied a 10^th percentile threshold rule to the predicted suitability, thus gray pixels in these maps represent suitable habitat. Hollow circles show the location of occurrence records used in the models, and the dashed line indicates the position of the Tennessee River.

Table 4

Variables with the highest contribution (%C > 10) to the predictive power of the MaxEnt distribution models for P. brachyphona and its Northern and Southern clades. The range of values for which the predicted occurrence was higher than 50% is shown. Results of model validation were assessed by using the Akaike Information Criteria (AIC) and its standard deviation.

The clade-specific models showed that Northern P. brachyphona inhabits southwestern Pennsylvania to Tennessee and northern Alabama (Fig. 6B). The Southern clade is distributed from central Alabama up to northwestern Georgia and southwestern North Carolina, possibly into southeastern Tennessee (Fig. 6C). For the model of Northern P. brachyphona, the highest contribution came from vegetation type (47.8%), followed by elevation (19.5%) and temperature during the driest trimester (10.4%). The variables with the highest contribution to the Southern clade model were temperature of wettest trimester (25.7%), precipitation during the driest trimester (19.8%), and temperature during the driest trimester (27.0%, Table 4).

Clade-level ecological niche description and niche overlap.—In the predicted niche of Northern P. brachyphona, the predominant vegetation was classified within the Eastern cool temperate ruderal/deciduous, South-Central mesophytic, and Allegheny/Cumberland dry oak forests (Supplemental Table 5; see Data Accessibility). In addition, Northern clade P. brachyphona was predicted to inhabit areas occurring at elevations from ∼230–830 masl, and with average temperatures around 0°C during the driest months of the year (Fig. 7, Table 4).

Fig. 7

Response curves of variables with highest contribution (>10%) to the Northern (black line, dark gray shadow) and Southern (gray line, light gray shadow) clade SDMs. The solid lines represent average probability of occurrence based on ten model replicates, while shadowed areas equal to ± standard deviation. The dotted line indicates a 50% occurrence probability, with environmental values above this threshold indicating suitable habitat. Each clade shows a peak of habitat suitability (>50% occurrence probability) at different values of each variable, suggesting ecological divergence.

In the case of the Southern clade, frogs were predicted to occur in areas where the mean temperature ranged from 0–10°C during the wettest trimester, precipitation between 270–500 mm, and temperature around 22°C during the driest trimester (Fig. 7, Table 4). In contrast to northern P. brachyphona, our model predicted the Southern clade to occur in mixed evergreen and deciduous forests, with pines being a common element in these vegetation categories. Common forest categories for the Southern clade SDM included the Southern piedmont dry pine forest, East Gulf coastal plain dry upland hardwood forest, and the Central interior and Appalachian riparian forest (Supplemental Table 5; see Data Accessibility). The Southern and Northern clade habitats shared the Allegheny/Cumberland dry oak forest category. Our models also predicted suitable habitat for Southern clade frogs to exist on elevations ranging from 100–320 masl (Supplemental Table 5; see Data Accessibility).

The background tests for niche overlap tests using Schoener's D and Warren's I showed that environmental models were significantly distinct. Comparing the D and I estimates from the Northern clade model with a distribution of the same estimates from models based on random points from the Southern clade yielded no niche overlap for D (DNORTH = 0.22 vs. DDIST = 0.62, P < 0.05), nor for I (INORTH = 0.47 vs. IDIST = 0.86, P < 0.05; Fig. 8A, B). The complementary test (South vs. North random distribution) also resulted in no overlap based on D (DSOUTH = 0.22 vs. DDIST = 0.37, P < 0.05) or I (ISOUTH = 0.48 vs. IDIST = 0.66, P < 0.05; Fig. 8C, D). Response curves from our species-level model of P. brachyphona showed bimodal distributions for two of the five variables with the highest contributions (Supplemental Fig. 4, Supplemental Table 5; see Data Accessibility), which also supports the ecological divergence between Northern and Southern clades.

Fig. 8

The distribution (solid gray line) and mean (gray dashed line) of D and I (A, B) niche overlap estimates obtained from MaxEnt models, originating from randomly selecting occurrences (and associated environmental values) in the Southern population. Asymmetric background tests were used to compare these distributions with the empirical estimate from the Northern population (black dashed line). The same procedure was carried out for the reciprocal comparison, Southern vs. Northern (C, D). In both cases, the empirical estimates significantly differed from the randomized distributions, suggesting the absence of niche overlap.

Divergence of acoustic signals within P. brachyphona.—Our randomization tests showed significant differences for pulse rate (PR) and dominant frequency peak (DFP) between P. brachyphona recorded in Tennessee (Northern clade) and Alabama (Southern clade). Mean PR for Tennessee was 84.2 pulses/s and 88.5 pulses/s for Alabama (mean diff = –4.3 pulses/s, P < 0.01). Mean DFP was 2,456.3 Hz for Tennessee and 2,716.3 Hz for Alabama (mean diff = –260 Hz, P < 0.01). We did not find significant differences in PN between Tennessee (mean PN = 25.0 pulses) and Alabama (mean PN = 24.5 pulses) according to our randomization tests (Fig. 9). Our t-tests performed for the same variables showed agreement with the randomization tests. Differences were found in PR (t = 3.0, P = 0.007) and DFP (t = 3.6, P = 0.002), but not for PN (t = –0.8, P = 0.452) between Tennessee and Alabama.

Fig. 9

Differences in pulse rate (PR), number (PN), and dominant frequency peak (DFP) in acoustic signals recorded in P. brachyphona from Tennessee (TN) and Alabama (AL). Locality-specific distributions of PR, PN, and DFP (A–C). The thick bar in these boxplots represents mean value. Lower and upper edges of each box are 25^th and 75^th quantile, respectively, and lower and upper whiskers represent the minimum and maximum limits of the interquartile range. Outliers are represented as hollow points. Randomization tests (1,000 reps) were performed for each of these variables (D–E), yielding significant differences for PR and DFP between Tennessee and Alabama, but not for PN.

Is P. brachyphona more than one species?—In light of our results, we contend that the taxonomic status of the Northern and Southern clades of P. brachyphona should be updated. The arguments that support splitting this taxon and recognizing both clades at the species level are: (1) Northern and Southern clades show strong genetic divergence, which is indicated by the monophyly of each clade (supported also by Bayes factor delimitation), high FST, clear genetic clustering by clade, and correlation of genetic structure with geography. The high genetic divergence leads us to think these clades are allopatric as they occur on either side of the Tennessee River. (2) There is evidence of ecological divergence as observed in species distribution models with differences related to vegetation and climatic variables; however, the evidence points to this being an event of allopatric speciation. This result is also concordant with the geographic separation provided by the Blue Ridge Mountains and surrounding high elevation areas. Finally, (3) there are significant differences between clades in regard to acoustic signals, which are critical for species recognition in frogs. We do not claim that differences in call preclude hybridization (which occurs with P. feriarum); however, differences in acoustic signals can further develop into strong prezygotic barriers. Given these multiple lines of evidence, we therefore propose the split of P. brachyphona and the elevation of the Southern clade to the species level. The type locality of P. brachyphona is “in west Pennsylvania, near the Kiskiminitas River” of Armstrong or Westmoreland County, Pennsylvania, USA (Cope, 1889; Frost, 2020), which occurs within the range of the Northern clade. Thus, populations from the geographic area of the Northern clade (IN, KY, OH, PA, WV, TN) retain the name P. brachyphona and the Southern clade receives a new name. The type locality of the new species is assigned to Hale County, Alabama (see below; Fig. 3B). This site was chosen because we obtained a reasonable sized collection of vouchers, tissues, morphological measurements, and acoustic recordings from this site (n = 10) for characterizing the new species.

Fig. 10

Dorsal and lateral views of type specimens of P. collinsorum (A–D) and P. brachyphona (E–H). (A–B) Male P. collinsorum from Hale County, Alabama (UF 190162). (C–D) Male P. collinsorum from Hale County, Alabama (UF 190167). (E–F) Female P. brachyphona from Preston County, West Virginia (NCSM 100109). (G–H) Male P. brachyphona from Preston County, West Virginia (NCSM 100110). Black bars equal to 1 cm. For a color version of this figure, please refer to the online version of this article.

Fig. 12

Oscillograms (first and second columns), spectrograms (third column), and power spectra (fourth column) of advertisement calls from (A) P. collinsorum (recorded at 10.2°C) and (B) P. brachyphona (recorded at 12.6°C). The second columns represent a single call extracted from those in the first column. In the case of P. collinsorum (A), two males were recorded while duetting and are identified with numbers.

Fig. 11

Live specimens of P. collinsorum (left, photo by EML) from Lawrence County, Alabama and P. brachyphona (right, photo by Suzanne Collins) from Harrison County, West Virginia. Specimens of P. collinsorum often lack a distinct dorsal coloration pattern. For a color version of this figure, please refer to the online version of this article.

DISCUSSION

Previously, Lemmon et al. (2007a) highlighted the need for analysis of additional lines of evidence to determine whether the Mountain Chorus Frog, P. brachyphona, included cryptic lineages. In this study, we have not only provided evidence from the behavior and ecology of this species, but also exploited the power of targeted hybrid enrichment of genome-wide sequence data to define genetic breaks (Lemmon et al., 2012; Lemmon and Lemmon, 2013). In agreement with the mitochondrial data set from Lemmon et al. (2007b), our nuclear phylogenetic analysis supported a clear split between Northern and Southern clades of P. brachyphona. Both clades were statistically supported by coalescent and maximum likelihood methods. Genetic clustering using >21K SNPs also supported the existence of two population clusters with significant genetic divergence, for which divergence is not explained simply by an isolation-by-distance process. There was, however, evidence of hybridization with the congeneric P. feriarum. Finally, we found evidence for divergence in acoustic signals between the two clades (based on one recording location per clade), as well as divergence in their ecological niches as they showed almost no overlap. We acknowledge that the genetic and acoustic data sets are somewhat limited in sample size; however, the agreement shown by the different lines of evidence led us to propose the reconsideration of the taxonomic status of these two clades and to describe a new species of chorus frog from the southern Appalachian Mountains, P. collinsorum.

Genetic evidence supporting divergence within P. brachyphona.—Previous studies on the phylogenetic relationships within Pseudacris pointed to the necessity of additional genetic data to support a split within P. brachyphona (Lemmon et al., 2007b). Based on mitochondrial sequences, Lemmon et al. (2007b) found two well-supported clades geographically corresponding to the Northern and Southern portions of the distribution of P. brachyphona. Our phylogenetic analysis based on nuclear sequences were consistent with this finding. The Northern clade, referred to henceforth as the Appalachian Mountain Chorus Frog (P. brachyphona), was represented in our phylogenetic analysis by samples originating from the states of Ohio, West Virginia, Kentucky, Tennessee, and northern Alabama. The Southern clade, referred to hereafter as Collinses' Mountain Chorus Frog (P. collinsorum), included individuals south of the Tennessee River mostly from Alabama and Georgia. Although more comprehensive sampling is required to accurately detect hybridization between these species, no evidence for admixture was detected in our data set (Fig. 3A, B), and a high coefficient of differentiation between species was observed (FST > 0.5). This pattern suggests the existence of a significant barrier to gene flow, which we propose may be the Tennessee River and adjacent higher elevations of the Blue Ridge Mountains.

Although P. collinsorum reaches into the lower elevation Piedmont of the southeastern U.S., both P. brachyphona and P. collinsorum are primarily distributed throughout the Appalachian Mountains, hence their common names (Mountain Chorus Frogs). Nevertheless, within this range, these species are not known to occur in the Great Smoky Mountains (Hoffman, 1980). Since divergence between the Mountain Chorus Frogs and their sister taxon, P. brimleyi, occurred in the Pliocene well after the uplift of the Appalachians (Lemmon et al., 2007a), we rather propose that the dynamism of climatic regimes during the Pleistocene and the concomitant generation of geographic barriers (Mayden, 1988; Rissler and Smith, 2010) resulted in divergence between the two species of Mountain Chorus Frogs. Because the phylogeographic break separating P. brachyphona and P. collinsorum coincides with the Tennessee River, it is possible that divergence between Mountain Chorus Frogs occurred after formation of this river. At the beginning of the Pleistocene, the Tennessee River formed from the confluence of the hypothetical Appalachian River and the Lower Tennessee River (Mayden, 1988; Kozak et al., 2006). This split adds an additional example to the many documented genetic breaks occurring in the vicinity of the Cumberland Plateau (Corser, 2008).

Further support for the idea that a barrier contributed to the genetic division observed between P. brachyphona and P. collinsorum comes from our observation that isolation-by-distance did not explain the break. Our analysis shows that divergence is higher between specimens of each species than between those within each species. We found that isolation by distance occurs within each species, as is expected for organisms with low vagility such as amphibians (Hillman et al., 2013).

Evidence of hybridization.—Although hybridization between P. brachyphona and P. collinsorum was not observed in our data set, the fact that we found evidence of hybridization between P. brachyphona and the more distantly related P. feriarum suggests that it could possibly occur between the two mountain chorus frogs in areas where they may overlap (Fig. 3). In fact, Lemmon et al. (2007a) previously observed evidence of introgression between P. brachyphona and P. feriarum, and between P. collinsorum and P. feriarum based on mitochondrial data. This finding is not surprising since most species in the trilling chorus frog clade (sensu Moriarty and Cannatella, 2004) occasionally hybridize with other close relatives at the periphery of their ranges (Lemmon et al., 2007a, 2007b). Furthermore, hybrids between P. brachyphona and P. feriarum have been created and raised to sexual maturity under laboratory conditions, suggesting that hybrids could also reach the same stage in nature (Mecham, 1965). Nevertheless, it should be noted that a larger sampling is required, especially throughout the Tennessee River basin, to accurately assess the extent of hybridization between P. collinsorum and P. brachyphona, if it occurs.

Niche divergence between clades of P. brachyphona.—As expected merely from the latitudinal differences between P. brachyphona and P. collinsorum, niche overlap tests yielded no similarity among SDMs. Pseudacris brachyphona has been a species traditionally associated with areas in proximity to woodlands (Barbour, 1957; Forester et al., 2003), which holds true for both clades according to our SDMs. Nevertheless, with respect to vegetation some differences were observed. Predicted habitat for P. brachyphona included forest categories with higher overall diversity of deciduous trees in which oaks (Quercus spp.), maples (Acer spp.), and hickories (Carya spp.) were a dominant element. Our SDMs for P. collinsorum showed that forest categories where this species occurs contain a higher richness of pine species than in habitat of P. brachyphona. The forest habitats of Pseudacris collinsorum were predicted to be generally drier than those of P. brachyphona. Pinus spp., Juniperus virginiana, and oaks were predicted to be common in habitat of P. collinsorum.

Response curves of SDMs generated by placing P. brachyphona and P. collinsorum as a single species showed bimodal distributions for the mean temperature during the driest and wettest seasons (Supplemental Figure 4; see Data Accessibility). The two habitat suitability peaks reflect the differences in habitat temperature and vegetation between P. collinsorum and P. brachyphona, and although do not necessarily imply local adaptation and ecological divergence, further research is warranted to test if fitness of each species is conditioned to their own habitat characteristics. It is also noteworthy that geographical overlap between the SDMs was not observed and that the break in distribution approximately aligned with the southern portion of the Blue Ridge Mountains. Although the Tennessee River represents the most likely vicariant barrier leading to divergence of the species, the paleogeographic importance of the Blue Ridge Mountains in the divergence between P. brachyphona and P. collinsorum merits further study.

Diversification of acoustic signals.—Significant differences in pulse rate and dominant frequency of calls recorded from Tennessee and Alabama suggest the presence of premating reproductive isolation between P. brachyphona and P. collinsorum (Coyne and Orr, 2004). Although these differences are based on single localities representing each clade, we contend the pattern provides further evidence of the distinctiveness of these species. Even small differences in pulse rate, for example, can be sufficient for frogs to discriminate between species (Platz and Forester, 1988; Gerhardt and Huber, 2002). The differences in dominant frequency, although compelling, are difficult to interpret in the light of unknown factors like size of the frogs, which predicts dominant frequency of the call and which is influenced by ecological selection (Lemmon, 2009). Pulse rate, however, is known to be critical for species recognition in other Pseudacris, and divergence of this call character, even among populations of the same species, can lead to reproductive isolation and incipient speciation (Lemmon, 2009).

Summary.—Our study supports the elevation of the Southern clade of P. brachyphona (sensu Lemmon et al., 2007b) to species level, and we name this new species Collinses' Mountain Chorus Frog (P. collinsorum). The diagnosis of this new species is supported by phylogenetic data, population genetic data, acoustic mating signals, and ecological information. Pseudacris collinsorum occurs at elevations of 100–320 masl across northern Alabama, northeastern Mississippi, and northwestern Georgia in the U.S.

DATA ACCESSIBILITY

Supplemental material, data set, and analysis parameter files are available at  https://www.copeiajournal.org/ch2020009.