Representatives of all 19 species of the hirtus-group were collected from 2012 to 2014. The Tennessee Wildlife Resources Agency permitted work in Tennessee (permit #1605). The Georgia Department of Natural Resources permitted work in Georgia (permit #8934). The National Park Service permitted collection of P. hirtus from Mammoth Cave National Park (permit #MACA-2013-SCI-0008). P. shapardi specimens from Oklahoma were collected by Matthew Niemiller. Two outgroup species representing the other main lineages in the subgenus Adelops (P. cavernicola and P. consobrinus) were collected in Florida. All beetles were collected by hand, typically with an aspirator or moist brush, and stored in 95% ethanol at -20°C. Seven species were collected from their type locality and several other species were collected from sites <1 km from their type locality. Sampling localities and species names are detailed in Table 1.

We amplified and sequenced three regions of the mitochondrial genome (cox1, cytb, rrnL-trnL-nad1) totaling on average 2300 bp. These regions were previously used in studies of palearctic Leiodidae (Ribera et al. 2010, Fresneda et al. 2011, Rizzo et al. 2013, Cieslak et al. 2014). DNA extractions and PCR amplifications followed standard protocols (Dixon and Zigler 2011). The primers used were based on Ribera et al. (2010) but modified based on transcriptome sequence information available from Ptomaphagus hirtus (Friedrich et al. 2011) (Table 2). Both strands of successful PCRs were sequenced on an ABI3730. Sequences were aligned and edited in Sequencher (v. 4.9, GeneCodes Corp).

We sequenced the three gene regions from 48 individuals from 29 populations, including all 19 species of the hirtus-group and two outgroups (P. cavernicola and P. consobrinus) representing the other Adelops species-groups (Peck 1973). In 19 cases, we sequenced two individuals from the same cave. For five species, we sequenced individuals from more than one cave. All sequences have been submitted to Genbank (KT167394-KT167537; Table 1).

The lack of indels or stop codons in the protein-coding cox1, cytb, and nad1 regions allowed for unambiguous multiple sequence alignment across all sites. The ribosomal and transfer RNA coding rrnL-trnL sequences were aligned using MAFFT (Katoh et al. 2005). Similar alignments were obtained using MUSCLE (Edgar 2004) and CLUSTAL (Larkin et al. 2007). Bayesian analysis was conducted with MrBayes v.3.2.2 (Ronquist and Huelsenbeck 2003), using four partitions (cox1, cytb, rrnL-trnL, nad1). Evolutionary models were estimated prior to the analysis with jModelTest v.2.1 (Posada 2008). GTR+I+Γ was the preferred model for all partitions except nad1 (TIM1+I+Γ) by AIC. MrBayes ran for 60 million generations using default values and the GTR+I+Γ model for each partition, saving every 6,000^th tree. 10% of values were discarded as burnin. Similar results were obtained using codon-position partitions for the three protein coding regions and from a two partition approach (protein coding and non-protein coding sequences); both yielded similar or identical topologies with the only differences occurring at poorly-supported nodes.

Maximum likelihood searches were conducted using RAxML (Stamatakis 2014) with a GTR+Γ model used for each of the four partitions. Branch support was evaluated by rapid bootstrapping over 1,000 replicates (Suppl. material 3). Maximum likelihood analyses yielded similar topologies and node support to the Bayesian analyses.

Divergence time estimates were generated with BEAST v1.8.1 (Drummond et al. 2012) using a total chain length of 60,000,000 generations with sampling every 1,000 generations. A Yule type speciation model with a wide normal distribution of tree height (treeModel.rootHeight = normal, initial: 15, Mean 15, Stdev 6, 95% of distribution between 5.495 mya and 24.89 mya) was used to model the tree. (We obtained similar results for tree height when using a uniform prior distribution of tree height between 5-25 mya.) The four data partitions shared a lognormal relaxed clock model with a substitution rate of 0.01 substitutions/site/my (stddev 0.002) (Rizzo et al. 2013). Each partition was assigned its own GTR+I+Γ substitution model. Tree topology was constrained to require monophyletic clades recovered with greater than 0.90 posterior probability in both MrBayes trees to appear in the final tree. We confirmed the convergence of parameter estimates by examination of three simultaneous runs using Tracer v1.5. Resulting tree files were summarized using TreeAnnotator v1.8.1, discarding 25% of samples as burnin and visualized using FigTree v1.4.2 (http://tree.bio.ed.ac.uk/software/figtree/).

Adopting the approach of Ribera et al. (2010), we produced a molecular phylogeny for the hirtus-group based on sequences from five genes across three mitochondrial regions, applying the molecular clock developed for the same genes in this subfamily of beetles (Leiodidae: Cholevinae) to establish divergence times for the phylogeny. Using that approach, the hirtus-group is estimated to have diverged from the consobrinus- and cavernicola-groups around 15 mya (Fig. 3). Within the hirtus-group we observed three stages of diversification. First, we found strong support for the soil-dwelling and geographically separated P. shapardi as sister species to the 18 cave-obligate species of the hirtus-group, with the two lineages diverging ~10 mya (Fig. 3). This divides the sole member of the hirtus-group from the Ozarks from the rest of the species in the southeastern United States (Fig. 2A). In the second stage of diversification, seven geographically distinct lineages originated between 8.5 and 6 mya (Fig. 3) across Kentucky, Tennessee, Alabama and Georgia, spreading hirtus-group members across several ecoregions (Fig. 2A). We were unable to resolve the branching order of these lineages with high confidence (Fig. 3). Six of these lineages are currently represented by a single species - P. barri, P. fiskei, P. hirtus, P. hubrichti, P. walteri and P. whiteselli (Figs 2A, 3). The seventh lineage contains all the species from the southern Cumberland Plateau (Fig. 2B), which form a well-supported clade (Fig. 3). In the third stage of diversification, the South Cumberlands lineage expanded from one to 12 species beginning around 6 mya. The branching order in the South Cumberlands clade is well-resolved (Fig. 3).

Figure 3. 

Ultrametric tree for the hirtus-group. Bayesian tree estimated from combined partial mitochondrial sequence data. Branches supported by posterior probability >0.90 are labeled with values or, for branches with posterior probability of 1.0, an asterisk. Blue bars indicate 95% confidence intervals of estimated ages for the nodes. Taxa are labeled with species name and specimen identifier (Table 1). Scale at bottom indicates divergence times in millions of years as estimated by BEAST (Drummond et al. 2012). Branch colors correspond to those in Figure 2A.

Most specimens sampled for a single presumed species formed monophyletic groups in our molecular phylogeny. This was not the case, however, for the P. hatchi specimens sampled from four different sites. Our molecular phylogenetic analysis reveals that P. hatchi constitutes a polyphyletic clade with respect to six species (P. chromolithus, P. episcopus, P. fecundus, P. laticornis, P. torodei and P. valentinei) with which it has overlapping or adjacent distribution ranges on the southern Cumberland Plateau (Figs 2, 3). One well supported clade joins several populations of P. hatchi with P. episcopus, P. fecundus and P. chromolithus (Fig. 3). A second well supported clade joins P. hatchi with P. laticornis (Fig. 3). We further noted that the intraspecific (between cave) divergence between P. hatchi populations was high, ranging from 1.2 to 5.1% P-distance for the cox1 gene (Suppl. material 2). Intraspecific (between cave) distances for populations of other Ptomaphagus species are all < 2.0% (Suppl. material 2). In several cases, the intraspecific P. hatchi divergence is greater than the interspecific divergence observed between other Ptomaphagus species.

To assess genetic diversification within and between lineages, we also surveyed intraspecific variation at population (within cave) and species (between cave) levels. In 19 cases (representing 14 species), we sequenced two individuals from the same cave. Variation between individuals from the same cave was low, with a mean cox1 p-distance of 0.15% (N = 19, range = 0.00 - 0.58%) (Suppl. material 1). In all 19 cases individuals from the same cave were each other’s closest relatives (Fig. 3).

For five species, we sampled animals from multiple caves. Specifically, we sampled P. hubrichti, P. loedingi and P. longicornis from two caves, P. barri from three caves, and P. hatchi from four caves. Intraspecific variation between caves had a mean cox1 p-distance of 2.1% for all pairwise comparisons (N = 12, range = 0.2 - 5.1%) (Suppl. material 2). Intraspecific variation between caves across the 2300 bp of all five sampled loci was slightly lower (mean p-distance = 1.7%), consistent with the previous observation that the cox1 region evolves faster than the other mitochondrial regions in this group (Ribera et al. 2010) (Suppl. material 2).

Our molecular phylogeny provides a time-calibrated picture of the diversification of the hirtus-group, shedding new light on the origins of an important component of cave biodiversity in North America. The integration of molecular and biogeographic data now suggests that the hirtus-group diversified in three stages. The first stage, occurring mid-Miocene ~10 mya, separated the sole extant eutroglophile in the hirtus-group, P. shapardi, located in the central United States from the 18 exclusively troglobiotic species located in the southeastern United States (Figs 2A, 3). The second stage of diversification spawned seven lineages across the southeastern United States during the late Miocene, 8.5-6 mya (Figs 2A, 3). In the third stage, further fine scale diversification occurred in the southern Cumberland Plateau region during the Pliocene and Pleistocene over the past six million years (Figs 2B, 3).

During the second stage of hirtus-group diversification, the lineages that diversified into the 18 cave-obligate members of the hirtus-group were established in five distinct geographic regions (Fig. 2A): the Mammoth Cave region of Kentucky, central Tennessee, northwest Georgia, the southern end of the Sequatchie Valley in Alabama, and the southern Cumberland Plateau in Tennessee and Alabama. Although further investigation of the deeper nodes of the hirtus-group phylogeny is still warranted, it appears that the two species in central Tennessee (P. barri and P. hubrichti) do not form a monophyletic group, nor do the two species in northwest Georgia (P. fiskei and P. whiteselli). This suggests that distinct hirtus-group lineages invaded underground habitats as many as seven times. Any alternative explanation requires one or more episodes of up to 75 km long-distance dispersal across significant non-karst terrain separating current species distributions. We do not favor long-distance dispersal as an explanation for current species distributions as the dispersal ability of hirtus-group members appears to be quite limited. All species are wingless and small, and none have been collected in surface habitats. Even short-distance migrations appear unlikely, as two species have never been collected in the same cave, which is particularly notable in the case of the cave- and species-rich southern Cumberland Plateau. In many cases, uninhabited cave habitats are present within just one or a few kilometers of known hirtus-group populations (Peck 1973), indicating extremely low vagility across the hirtus-group.

Further consistent with minimal dispersal ability is the fact that the southern Cumberland Plateau lineage diversified into twelve species over the past 6 million years and no members of this lineage appear to have migrated from the southern Cumberland Plateau. Notwithstanding the strong circumstantial evidence, long-distance migration cannot be completely discounted as an explanation for some aspects of the species distributions we see today. For instance, long-distance dispersal has been proposed for one group of troglobiotic leiodids in the palearctic realm (Troglocharinus of Spain; Rizzo et al. 2013). Rizzo et al. (2013) proposed these beetles expanded their range in the early Pliocene from the central Pyrenees across a significant non-karst region to coastal karst habitats 60–70 km away, via stepping-stone migration across the surface during a permissive climate period, followed by subsequent isolation of the two lineages (Rizzo et al. 2013).

The limited dispersal capacity of troglomorphic (microphthalmic and wingless) Ptomaphagus species suggests that troglomorphy developed multiple times through convergent evolution in this group. This scenario and its implications could be further scrutinized at the molecular level. Previous transcriptome analyses on P. hirtus, for example, identified several genes in the ommochrome eye pigmentation pathway that are no longer expressed in P. hirtus, consistent with the lack of pigment granules in the highly reduced eyelets of this species (Friedrich et al. 2011). Similarly reduced and non-pigmented eyelets are found in all other cave-adapted species of the hirtus-group. A single cave colonization event during early hirtus-group evolution would predict shared lack-of-function mutations underlying the regression of eyes and eye pigmentation across the hirtus-group. Multiple independent cave colonization events, in contrast, would be reflected by a lack of shared lack-of-function mutations in eye pigmentation genes. This approach, which can be extended to other cave adaptive traits, will ultimately require genomic analyses to determine whether identical mutations are present in geographically distinct populations (consistent with a single cave invasion and the evolution of troglomorphy followed by long-distance dispersal) or whether different mutations are present in different lineages (consistent with multiple cave invasions followed by the convergent evolution of troglomorphy). Consistent with independent cave invasions and evidence for the feasibility of such an approach, distinct mutations have been observed in the pigmentation gene cinnabar and in the opsin genes of subterranean diving beetles (Leys et al. 2005, Tierney et al. 2015) and in rhodopsin genes of cavefish (Niemiller et al. 2013) in geographically distinct populations of those groups.

The Climatic Relict Hypothesis suggests that a species’ initial colonization of caves occurred when it sought refuge from environmental stressors in the surface environment (reviewed in Culver and Pipan 2009). The Climatic Relict Hypothesis is often proposed for troglobiotic taxa lacking close relatives on the surface, with those relatives presumably extinct as a result of the changing surface conditions that forced the proto-troglobiotic taxa underground. The Climatic Relict Hypothesis can be supported by evidence that the cave species originated or diversified during a period of climate change. A convincing case for the Climatic Relict Hypothesis was developed by Leys et al. (2003) for Australian diving beetles where numerous independent invasions of isolated calcrete aquifers occurred. Leys et al. (2003) used a molecular clock to show that the timing of these invasions correlated with increasing aridity moving across the region from north to south during the Pliocene. In North America, the Climatic Relict Hypothesis has often been suggested as a driver of cave colonization related to climate change during late Pleistocene glaciation events (Culver and Pipan 2009). Peck (1973, 1984) offered the hirtus-group as a possible example of cave colonization via the Climatic Relict Hypothesis. To explain the species diversity of the hirtus-group, Peck (1973, 1984) suggested that the first Ptomaphagus to invade cave habitats subsequently dispersed to nearby caves during recent Pleistocene interglacial periods via cool and moist habitats such as leaf litter, talus or moss mats, after which the species became more completely cave adapted and isolated.

Our time-calibrated molecular phylogeny rejects the hypothesis of middle to late Pleistocene diversification in the hirtus-group. We found the timing of most speciation events in the hirtus-group to be an order of magnitude greater than Peck hypothesized (Peck 1973, 1984). Our divergence time estimates are similar to those derived from molecular data in other troglobiotic groups (e.g., Leys et al. 2003, Faille et al. 2010, Ribera et al. 2010, Derkarabetian et al. 2010, Niemiller et al. 2012). The divergence times estimated by these studies indicate that the species in these groups are much older than the recent Pleistocene. With most diversification of the hirtus-group occurring in the Miocene and Pliocene, we can rule out the Climatic Relict Hypothesis related to Pleistocene glaciations as a general rule for the hirtus-group.

Now recognizing that the third stage of hirtus-group diversification, the radiation in the southern Cumberland Plateau, began ~6 mya, we hypothesize that cave-adapted Ptomaphagus populations distributed throughout the region were isolated by vicariant events as the southern Cumberland Plateau eroded and fragmented over millions of years. Barr and Holsinger (1985) suggest gene flow between cave systems can be reduced via “simple erosion that divides a karst area into isolated segments by cutting down into underlying noncavernous strata.” This mechanism is particularly applicable to terrestrial troglobionts, whose populations are unlikely to maintain connectivity via deep groundwater connections, as might be the case in aquatic troglobionts (e.g., Fenolio et al. 2017). The southern Cumberland Plateau is greatly dissected, with numerous isolated and peripheral mountains and ridges, and hirtus-group species are frequently limited to a single isolated ridge. For example, P. longicornis is limited to caves on Keel Mountain (in Madison and Jackson Counties, Alabama), P. julius is limited to caves on July Mountain (Jackson County, Alabama), and P. solanum is limited to Tater Knob (Jackson County, Alabama). Each of these ridges is an isolated remnant of the Cumberland Plateau (Fig. 4). Further consistent with the vicariance by erosion model (Fig. 5), the earliest lineages to diverge in the South Cumberlands clade (P. julius, P. hazelae + P. solanum, P. loedingi and P. longicornis) are found in peripheral and isolated ridges at the edge of the southern Cumberland Plateau, whereas later-diverging lineages are concentrated in the more intact central region of the southern Cumberland Plateau (Figs 3, 4).

Figure 4. 

Distribution of Ptomaphagus species on the southern Cumberland Plateau, overlaid on a digital elevation model. Higher elevations (to 500 m) are indicated by darker shades, lower elevations (to 180 m) by lighter shades. Ptomaphagus species diverging early in the South Cumberlands lineage are limited to isolated ridges and mountains on the fringes of the plateau. These species are P. loedingi (yellow), P. longicornis (dark gray), P. julius (blue), P. solanum (dark green) and P. hazelae (light blue). The colors used here correspond to those in Figure 2B.

Figure 5. 

Biogeographic and phylogenetic expectations for a ‘vicariance by erosion’ scenario as hypothesized for the southern Cumberland Plateau. A–C Karst (gray) erodes and fragments over time, leading to the isolation and divergence of cave populations (colored circles) in the remaining patches of karst D A phylogeny consistent with the vicariance by erosion process, with taxa that diverge early distributed at the periphery of the eroding region.

Although we lack a clear picture of the timing and pattern of erosion and fragmentation of the southern Cumberland Plateau, it is clear that the ~6 my over which the South Cumberlands lineage diversified is sufficient for significant erosion and cave development to have occurred. In support of this, on the western edge of the Cumberland Plateau in middle Tennessee, extensive stream incision, erosion and cave development occurred over a similar period of time. In this region, the oldest caves (now located 60–90 m above current river level) were estimated to be 3.5–5.7 my old based on the dating of radioactive cave sediments (Sasowsky et al. 1995, Anthony and Granger 2004, 2007). Assuming a similar degree of erosion in the southern Cumberland Plateau, it is reasonable to assume that the hirtus-group has been heavily impacted by fragmentation, resulting in cave population isolation events over the past 6 my.

The hirtus-group is widespread in the Mammoth Cave region of Kentucky and in the southern Cumberland Plateau of Tennessee and Alabama (Fig. 2). Measured as extent of occupancy (Bachman et al. 2011), the range extent of hirtus-group species in the two regions is similar: 3,952 km² in the Mammoth Cave region and 4,151 km² in the southern Cumberland Plateau. It is striking that, across these similar ranges, only one species (P. hirtus) is present in the Mammoth Cave region, whereas 12 species are present in the southern Cumberland Plateau. As discussed above, the diversification in the southern Cumberland Plateau may be related to more frequent occurrence of vicariance by erosion, in contrast to the more continuous karst of the Interior Plateaus in the Mammoth Cave region. Indeed, such patterns have been proposed in the past (e.g., Barr and Holsinger 1985 and references therein).

One prediction from this scenario is that P. hirtus evolved less genetic diversity across its range than we observed across the South Cumberlands lineage. Alternatively, P. hirtus may represent a collection of cryptic species, as has been observed in numerous cave lineages (e.g., Trontelj et al. 2009, Derkarabetian et al. 2010, García-Machado et al. 2011, Hedin 2015, Zhang and Li 2013). Our analysis presented here sampled only a single P. hirtuspopulation (Table 1). Thus, further phylogeographic and population genetic studies of P. hirtus are warranted to clarify whether current species definitions are consistent with the genetic diversity present in the region. Such studies are further recommended given the evidence of cryptic species in what has previously been considered a single taxon, P. hatchi.

With the exception of P. hirtus, P. hatchi has the largest range extent of any hirtus-group species. P. hatchi is also the only member of the hirtus-group with a range that overlaps those of other species in the group. We found that the current species definition of P. hatchi is polyphyletic with respect to six species (P. chromolithus, P. episcopus, P. fecundus, P. laticornis, P. torodei and P. valentinei) from the southern Cumberland Plateau lineage (Figs 2, 3). Peck’s (1983, 1984) study of species boundaries within P. hatchi relied on distinctions in the female spermatheca (Peck 1984) and evidence of reproductive isolation in various interpopulational crosses (Peck 1983). In the most recent taxonomic revision, Peck (1984) restricted P. hatchi to populations with ‘form I’ spermathecae and described P. laticornis (with ‘form II’ spermathecae) and P. chromolithus (with ‘form III’ spermathecae). He also described P. torodei as a close relative of P. valentinei and raised P. fecundus to full species status (Peck 1984).

Significantly, our mitochondrial sequence divergence data do not align with these distinctions. This could be due to a more dynamic variability of spermatheca form than previously envisioned. Alternatively, ‘form I’ spermathecae may represent an ancestral state that has been retained in some but not all descendant lineages (now called P. hatchi), leading to the erroneous support for a polyphyletic taxon based on the shared similarity of a plesiomorphic character state, i.e. symplesiomorphy. It is also possible that hybridization and introgression have confused the molecular phylogenetic picture of these lineages as we only analysed mitochondrial DNA. This, however, seems unlikely as two species of Ptomaphagus (or two spermathecal forms) have never been reported from the same cave (Peck 1984). At present, the available molecular data are too limited to recommend taxonomic revision, but further molecular analysis of these taxa across their distributions has now become essential to gain a reliable understanding of the hirtus-group diversification.

Exhibiting high levels of taxonomic diversity and endemism, the southern Cumberland Plateau is a hotspot for cave biodiversity. This region has been compared to other centers of cave biodiversity such as the Dinaric karst of Slovenia, the French Pyrenees and the Mammoth Cave region (Culver et al. 2006). This biodiversity peaks in a six-county region in northeast Alabama and south-central Tennessee. Totaling ~10,000 sq km, this region has more than 4400 known caves that support more than 150 troglobionts (Culver et al. 2000, Niemiller and Zigler 2013). Members of the hirtus-group compose nearly 10% of the troglobionts present in the southern Cumberland Plateau.

With a few exceptions, we lack time-calibrated molecular phylogenies for troglobionts from the southern Cumberland Plateau. This has limited our understanding of the development of biodiversity in this cave biodiversity hotspot. Important exceptions exist for several aquatic taxa including crustaceans, fish and salamanders (Buhay and Crandall 2005, and comment in Trontelj 2007, Niemiller et al. 2008, Niemiller et al. 2012). These studies indicate diversification during the Pleistocene for all of these taxonomically diverse groups. Here, we provide the first time-calibrated study of a diverse terrestrial troglobiont group from the region. We found that the diversification of the South Cumberlands lineage of the hirtus-group extended across ~6 my, beginning near the end of the Miocene and continuing through the Pliocene and Pleistocene (Fig. 3). Thus, it appears troglobionts have accumulated in the southern Cumberland Plateau since at least the late Miocene.

We therefore suggest that vicariance by erosion played a generally significant role in the accumulation of troglobiotic species in this biodiversity hotspot. Better models of the fragmentation of the southern Cumberland Plateau and the development of time-calibrated phylogenies for other species-rich terrestrial troglobiotic taxa from the region (such as pseudoscorpions and millipedes) would allow further evaluation of “vicariance by erosion” model (Fig. 5) for the southern Cumberland Plateau cave biodiversity hotspot.