Research Article
Research Article
Geographically structured genetic diversity in the cave beetle Darlingtonea kentuckensis Valentine, 1952 (Coleoptera, Carabidae, Trechini, Trechina)
expand article infoOlivia F. Boyd§, T. Keith Philips, Jarrett R. Johnson, Jedidiah J. Nixon
‡ Western Kentucky University, Bowling Green, United States of America
§ Oregon State University, Corvallis, United States of America
Open Access


Cave beetles of the eastern USA are one of many poorly studied groups of insects and nearly all previous work delimiting species is based solely on morphology. This study assesses genetic diversity in the monotypic cave carabid beetle genus Darlingtonea Valentine 1952, to test the relationship between putative geographical barriers to subterranean dispersal and the boundaries of genetically distinct groups. Approximately 400bp of the mitochondrial cytochrome oxidase I (COI) gene was sequenced from up to four individuals from each of 27 populations, sampled from caves along the escarpments of the Mississippian and Cumberland plateaus in eastern Kentucky, USA. The 81 individuals sequenced yielded 28 unique haplotypes. Hierarchical analyses of molecular variance (AMOVA) within and among geographically defined groups tested two a priori hypotheses of structure based on major and minor river drainages, as well as genetic distance clusters defined a posteriori from an unrooted analysis. High genetic differentiation (FST) between populations was found across analyses. The influence of isolation by distance could potentially account for much but not all of the variation found among geographically defined groups at both levels. High variability among the three northernmost genetic clusters (FCT), low variability among populations within clusters (FSC), and low within-cluster Mantel correlations indicate the importance of unidentified likely intra-karst barriers to gene flow separating closely grouped cave populations. Overall phylogeographic patterns are consistent with previous evidence of population isolation among cave systems in the region, revealing geographically structured cryptic diversity in Darlingtonea over its distribution. The landscape features considered a priori in this study were not predictive of the genetic breaks among the three northern clusters, which are genetically distinct despite their close geographic proximity.


mitochondrial DNA, Mississippian Plateau, Pennyroyal, phylogeography, Southern Appalachians, troglobites, troglobionts


Variation within a species is usually not random, but structured in some way and typically forms a metapopulation with various levels of deviation from panmixia (Hanski 1999). Landscape features that correlate with intraspecific variation may represent boundaries reducing gene flow among discrete groups of populations. Alternatively, differences between populations of a species may increase linearly with physical distance, especially for less vagile organisms (e.g., Lee and Mitchell-Olds 2011, Goudarzi et al. 2019). The limestone karst regions of the Eastern United States support a remarkable diversity of cave-specialized animals (Barr 1985, Peck 1998, Hobbs 2012, and see White et al. 2019). Troglobionts, i.e., obligate and permanent cave inhabitants, can be predicted to demonstrate high levels of population genetic structure owing to a lack of gene flow between caves. Even long-term population isolation, however, may not yield diagnosable morphological differentiation due to phenotypic convergence in similar cave environments (Wiens et al. 2003, Derkarabetian et al. 2010, Hedin and Thomas 2010). Therefore, many troglobiotic taxa may harbor cryptic variation (Niemiller et al. 2012), and the biodiversity of cave-dwelling organisms may currently be underestimated.

Patterns of gene flow among caves in karst areas vary mostly in accordance with the geographical distribution of subterranean limestone (e.g., Caccone 1985, Katz et al. 2018). In limestone-rich parts of the Eastern United States (Fig. 1) where karst exposure is patchy, structurally fragmented, and discontinuous, caves are generally smaller and more isolated from one another (e.g., Currens 2002, Christman et al. 2005). One such region is the Appalachian Valley (AV), located primarily in eastern Tennessee and Virginia, which supports a high diversity of endemic cave beetles and other troglobites per unit area, many of which are limited in range to one or a few caves (Barr 1967, 1981, 1985, Christman et al. 2005, Niemiller and Zigler 2013). Conversely, troglobiotic invertebrates that inhabit large and highly interconnected cave systems which have permeated the large and uninterrupted exposures of limestone in the Mississippian Plateau (MP) region have comparatively broader ranges and less predictable distributional boundaries (e.g., Barr 1979). Species numbers and abundances differ among cave communities in the interior low (“Mississippian”) plateau (referred to below as “MP”) and Appalachians (Appalachian valley and ridge, referred to below as “AV”) regions; MP cave systems support larger and richer communities of troglobionts compared with those in the AV to the east (Barr and Holsinger 1985). With fewer endemics per unit area, cave species in the MP have been suggested as more likely to occur in sympatry than those inhabiting AV caves (Barr 1967, 1985, 2004). More recently though, Christman et al. (2016) presented contrasting evidence that despite the greater dissection of karst in the AV, cave species actually have lower rates of endemicity in the AV than in the MP.

The cave-rich limestone of the MP is bisected by the Cumberland Saddle, a low point in the Cincinnati Arch formation, which separates the MP into two regions: the MP-I to the west and the MP-II to the east (Fig. 1). Within both bands of the MP, cave interconnectivity has helped establish and maintain diversity by facilitating subterranean dispersal, leading to extensive range overlap and sympatry of species that were previously isolated, and linking populations together through gene flow which likely has reduced stochastic extinction events (Barr 1985, Barr and Holsinger 1985).

Isolating barriers between cave systems restrict gene flow and promote divergence among populations of cave organisms, effectively dividing parts of cave systems into subterranean islands (Culver 1970). Major waterways like the Cumberland and Ohio rivers serve as important fluvial barriers to dispersal of terrestrial troglobionts (Barr and Holsinger 1985, Barr 1985) and even some stygobionts (Niemiller et al. 2013), but smaller streams and rivers may actually promote their dispersal; Barr (1985) compared the “meander frequencies” of rivers dividing the distributions of cave beetle species, finding support for his hypothesis that the more turns a river takes over a given distance, the more often beetles washed out of caves will survive to encounter limestone outcrop karst refugia leading to an increase in distribution range via colonization of new cave systems.

Figure 1. 

(Adapted from Barr 1985, Figure 3) Map showing the major geologic features important for cave development in the southeastern United States: MP-I and MP-II (green) are western and eastern bands of the Mississippian Plateau. Dots indicate collecting records (see Figure 3).

Study species

Darlingtonea Valentine, 1952 is a monotypic genus of cave carabid beetle found in a narrow distributional band from north-central Tennessee (known from a single cave near the Kentucky border) extending northeastward into east-central Kentucky (mainly the northern part of “MP-II” in Fig. 1 and see Fig. 3). Like many of the other cave-specialized carabids of the subtribe Trechina, Darlingtonea are true troglobionts, with adaptations for subterranean life: they lack eyes and wings, possess enlarged mouthparts, lengthened appendages, and specialized sensory setae, and are depigmented compared with their epigean relatives (Fig. 2). Darlingtonea kentuckensis Valentine is usually abundant in caves within its range compared to many species of closely related Pseudanophthalmus (Valentine 1952). Molecular phylogenetic evidence from a 2012 study including representatives of all five eastern North American cave genera shows the genus shares common ancestry with a lineage of Pseudanophthalmus and is essentially derived from within the latter (Philips and Valkanas, unpublished). The close relationship of those genera together with Ameroduvalius Valentine, Nelsonites Valentine, and Neaphaenops Jeannel within the Trechoblemus series and within the Trechina is also strongly supported by Maddison et al. (2019).

Figure 2. 

Gravid female Darlingtonea kentuckensis photographed in Fletcher Spring Cave, Rockcastle County, Kentucky. Photo courtesy of Dr. Matthew Niemiller, University of Alabama, Huntsville.

Figure 3. 

Cave localities of currently known sites for Darlingtonea kentuckensis. White dots were the caves sampled for this study while black dots represent caves unsampled.

Regarding the origin and diversity of North American cave trechines, most authors have favored some version of a “Pleistocene-effect” model (Holsinger 1988). In contrast, Faille et al. (2015) puts the divergence times between two European trechine Aphaenops cave species around 9 my (with a credibility range of 4–17 my). Regardless of age, the proposed evolutionary scenario can be summarized as follows: As climate cycles associated with glacial advance and recession led to fluctuation of surface conditions, ancestral trechines followed cool, moist microhabitats from the deep soil which was abundant during glacial maxima to subterranean or montane refugia during warmer, drier glacial minima (Barr 1969, 1971, 1973, 1985). Periods of isolation in caves during warm intervals were punctuated by periods of introgression during cool intervals until a warm, stable post-Pleistocene climate restricted surface dispersal and promoted subterranean allopatric speciation (and see Jeannel 1948, 1949 for further details on the effects of glaciation).

Other authors have found isolation and divergence in allopatry to be an unsatisfactory model for cave colonization in other taxa, which may be better viewed as a parapatric ecological transition or ‘adaptive shift’ occurring in the presence of gene flow via diversifying selection (Niemiller et al. 2008). Further, surface characteristics of the Earth, such as latitude, percent karst, and landscape rugosity (Topographic Position Index) may have significant effects on the evolution of a cave-adapted fauna (Christman et al. 2016)

It is currently unclear what factors have led to the evolution of any morphological or genetic diversity within Darlingtonea kentuckensis. Darlingtonea kentuckensis has a broader than average distribution compared to most terrestrial Eastern North American troglobionts based on our review (Philips et al. unpublished). Both Valentine (1952) and Barr (1985) noted some morphological diversity among populations of D. kentuckensis. For example, Valentine noted subtle differences including a slightly more convex body form, slightly wider elytra, and more rounded elytral humeral angles (in populations on either side of the Cumberland River), but concluded there was not enough support for subspecific designation. In contrast, the population from Big Saltpeter Cave in Rockcastle County by the Rockcastle River was thought to be distinct enough to warrant the subspecific name D. k. lexingtoni Valentine. Morphologically, this taxon diagnosis was based on a slightly paler body color, very slightly narrower pronotum, flatter elytral disc, and claimed differences in the male genitalia that included subtle differences in the apex of the median lobe and one lobe of the internal sac (see Valentine 1952, Plate IV).

Barr (1985) speculated that D. kentuckensis includes at least seven subspecies or races isolated by landscape barriers. Kane et al. (1992) sampled ten D. kentuckensis populations from across the MP-II for a study of allozyme diversity. Polymorphism in nine of the eleven electrophoretic markers examined combined with the lack of variation within populations and high FST across loci suggested long-term isolation.

The exceptional species diversity in North American cave trechines (Peck 1998) makes this lineage valuable to understanding the speciation processes in troglobitic insects and other terrestrial cave organisms. Since populations of Darlingtonea occur across a broad geographic range relative to other troglobiotic taxa while belonging to a single morphologically, geographically, and genetically distinct lineage, D. kentuckensis is a convenient model for comparing observed patterns of genetic variation against those predicted by a climate-mediated process of cave colonization.

Purpose and hypotheses

If important barriers to dispersal for cave trechines in the MP-II region exist, hierarchical tests of population genetic structure should reveal a general pattern of low diversity within and high diversity among clusters of genetically similar populations. Specific geographic barriers between these genetic clusters that may be responsible for population structure can then be hypothesized and should make geographic sense without being purely attributable to the influence of isolation and genetic divergence by distance. Patterns may also reveal the presence of cryptic species or subspecies.

The Kentucky and upper Cumberland rivers represent the two primary watersheds in the MP-II. Further, the divide between the watersheds of the Kentucky and Rockcastle rivers in northern Jackson County (Barr 1985) and the upper Cumberland River in southern Pulaski County (Barr 1985, Lewis and Lewis 2005) may represent two additional major barriers to gene flow. These drainage barriers, along with an additional geological/historical barrier isolating genetically distinct groups of populations in northern and southern Pulaski County (Kane et al. 1992), may effectively divide the sampled range of Darlingtonea into four faunal regions (Table 1 and Fig. 4): on the north side (Faunal Region 1) or the south side (Faunal Region 2) of the Kentucky-Rockcastle drainage divide and north (Faunal Region 3) and south (Faunal Region 4) of the Cumberland River. Populations hypothesized by Barr (1985) from a potential fifth faunal region east of the Big South Fork of the Cumberland River were not sampled in this study. “Structure hypothesis I” tested herein predicts that sampled populations fall into four genetically distinct clusters that are geographically consistent with the hypothesis of reduced gene flow among these four major regions subdivided by major river systems.

Caves also fall into smaller, “minor” watersheds (Table 1) that could define components of population genetic structure at a finer resolution, especially if Barr’s (1985) hypothesis about the role of smaller, meandering streams in promoting cave beetle dispersal is valid. Samples from the 27 localities (each from an individual collecting event) in the final data set were assigned to watersheds based on both absolute proximity to second- and third-order streams and qualitative topographic information. Under “structure hypothesis II”, populations are expected to fall into ten genetically distinct clusters, with a pattern of genetic structure that is geographically consistent with reduced gene flow among these ten minor watersheds.

Table 1.

List of Darlingtonea populations included in a study of mitochondrial haplotypes, including population (taxon) reference codes, locality information, collection dates, sample size, faunal region, local watershed and GenBank accession codes. Faunal region 1.

Taxon Code Cave County Collection Date N Faunal Region Local Watershed/Code (River Drainage) GenBank accession number
BLO Blowing Wayne 1-Mar-2014 3 4 Otter Creek/OT (CR) MN880837, MN880838, MN880839
CLF Clifford Pearson Estill 14-Aug-2014 2 1 Station Camp Creek/SC (KR) MN880814, MN880815
CLI Climax Rockcastle 31-Jul-2014 4 2 Roundstone Creek/RO (RR) MN880810, MN880811, MN880812, MN880813
FLE Fletcher Spring Rockcastle 15-Mar-2014 3 2 Skegg Creek/SK (RR) MN880827, MN880828, MN880829
GSP Great Saltpeter Rockcastle 15-Aug-2014 4 2 Roundstone Creek/RO (RR) MN880817, MN880818, MN880819, MN880820
HIC Hicksey Jackson 14-Aug-2014 4 1 Station Camp Creek/SC (KR) MN880806, MN880807, MN880808, MN880809
HIS Hisel Jackson 1-Aug-2014 1 1 Station Camp Creek/SC (KR) MN880805
HRT Hurt Wayne 12-Jul-2014 4 4 Beaver Creek/BE (CR) MN880846, MN880847, MN880848, MN880849
JES Jesse Wayne 28-Sep-2013 4 4 Otter Creek/OT (CR) MN880836, MN880840, MN880844, MN880845
JGR John Griffin Jackson 31-Jul-2014 4 2 Horse Lick Creek/HL (RR) MN880801, MN880802, MN880803, MN880804
KOG Koger Wayne 28-Sep-2013 1 4 Beaver Creek/BE (CR) MN880850
LAI Lainhart #1 Jackson 1-Aug-2014 4 1 Station Camp Creek/SC (KR) MN880798, MN880799, MN880800, MN880816
LAK Lakes Jackson 31-Jul-2014 3 2 Horse Lick Creek/HL (RR) MN880792, MN880796, MN880797
MOR Morning Hole Jackson 14-Aug-2014 2 1 Station Camp Creek/SC (KR) MN880794, MN880795
MUL Mullins Spring Rockcastle 15-Mar-2014 2 2 Roundstone Creek/RO (RR) MN880821, MN880822
PHC Pine Hill Rockcastle 15-Mar-2014 3 2 Roundstone Creek/RO (RR) MN880830, MN880831
PIN Piney Grove Pulaski 20-Oct-2013 3 3 Pitman Creek/PI (CR) MN880855, MN880856, MN880857
POU Pourover Pulaski 20-Oct-2013 4 3 Buck Creek/BU (CR) MN880858, MN880859, MN880860, MN880861
RCH Richardson’s Pulaski 20-Oct-2013 4 3 Pitman Creek/PI (CR) MN880866, MN880867, MN880868, MN880869
ROA Roadside Pulaski 4-Jul-2012 1 3 Pitman Creek/PI (CR) MN880862
SAV Savage (Copperas Saltpeter) Clinton 28-Sep-2013 2 4 Spring Creek/SP (CR) MN880834, MN880835
SOR Sinks of Roundstone Rockcastle 15-Aug-2014 2 2 Roundstone Creek/RO (RR) MN880832, MN880833
SRI Sinks and Rises Jackson 1-Aug-2014 3 2 Horse Lick Creek/HL (RR) MN880790, MN880791, MN880793
STA Stab Pulaski 4-Jul-2012 4 3 Buck Creek/BU (CR) MN880851, MN880852, MN880853, MN880854
STL Steele Hollow McCreary 12-Jul-2014 3 4 Little South Fork/LS (CR) MN880841, MN880842, MN880843
TEA Teamers Rockcastle 15-Aug-2014 4 2 Roundstone Creek/RO (RR) MN880823, MN880824, MN880825, MN880826
WIND Wind Pulaski 4-Jul-2012 4 4 Pitman Creek/PI (CR) MN880863, MN880864, MN880865, MN880870



Collecting localities (Figs 1, 3) were prioritized based upon a technical report compiled by Harker and Barr (1979) for the Kentucky State Nature Preserves Commission that listed caves where the target taxon could be sampled. Inclusion of several additional localities that would have benefited this study was not possible due to cave access restrictions imposed in recent decades by landowners for the prevention of vandalism or by conservation authorities for the protection of the two federally endangered Myotis bat species. Appropriate measures were taken as recommended by the most recent national White Nose Syndrome decontamination protocol (v.06.25.2012) to help slow the spread of Geomyces destructans Blehert & Gargas (also known as Pseudogymnoascus destructans (Blehert & Gargas) Minnis & D.L. Lindner) the introduced fungal pathogen which has led to recent population declines in many species of North American bats.

Beetle specimens were collected by hand into 95% ethanol and placed at -20 °C for short-term storage within 48 hours of collection. Ethanol was changed after processing (individuals from each locality were sorted by genus and inventoried) and whole specimens from each location were stored together in 95% EtOH at -80 °C. Table 1 summarizes collecting information and group membership relative to each hypothesis.


Depending on the number of specimens available, up to four Darlingtonea individuals per cave were sequenced (for a total of 81 specimens) to capture a sample of within-population mitochondrial cytochrome oxidase subunit I (COI) haplotype diversity (Table 1). Gut material, if visible, was removed in order to avoid amplification of foreign DNA from prey or other organisms. Whole specimens were ground inside 1.5 ml tubes using sterile plastic pestles and incubated in a solution of CTL buffer and proteinase K for 18–24 hours at 40 °C. Total genomic DNA was extracted from whole specimens using an E.Z.N.A. Insect DNA kit from Omega Bio-Tek. Nucleic acid concentration and purity was quantified using a NanoDrop 2000 spectrophotometer. Extractions were stored post-purification at -80 °C for long-term DNA preservation.

An ~850 bp COI target region was amplified from genomic DNA using the primer pair “Pat” and “Jerry” (Simon et al. 1994). Thermal cycling conditions for polymerase chain reaction (PCR) followed those specified by the manufacturer of TaKaRa Ex Taq, which was used for all PCR reactions. Primer annealing temperatures were optimized qualitatively by visualizing PCR products from a temperature gradient on an agarose gel to maximize yield and limit nonspecific binding. A QIAquick Gel Extraction Kit from Qiagen was used to purify most PCR products before sequencing. DNA template samples were prepared for sequencing in the forward direction using a BigDye Terminator v3.1 Cycle Sequencing Kit (Applied Biosystems). Automated cycle sequencing was performed using an ABI 3130 Genetic Analyzer (Applied Biosystems) at Western Kentucky University.

Sequences were aligned using CLUSTALW (Larkin et al. 2007) using the default settings (gap open cost of 15 and a gap extend cost of 6.66), although no gaps were present. Sequences were then edited manually in Geneious version R7 (, Kearse et al. 2012) according to the following rules: IUPAC ambiguous bases were inserted where peaks in the chromatogram overlapped, making base calls questionable. The ends of sequence reads were trimmed when peaks became indistinct or read quality (%HQ) consistently fell below 20 percent (this was common among reads, especially at the 3’ end, since sequencing was performed in only one direction). Reads were translated and screened for signs of pseudogene amplification, including mid-sequence stop codons and frameshifts. Each offending read was manually inspected: in cases where the correct base was obvious upon inspection of the chromatogram, the sequence was corrected and included; in cases where the correct base was unclear, the sequence was omitted and sequencing was re-attempted for that specimen. All sequences were trimmed evenly to 413 bp to eliminate the considerable variation in sequence length that resulted from quality trimming while maximizing the number of operational taxonomic units (OTUs) included.


Partial COI sequences were collapsed into haplotypes using the online tool FaBox (Villesen 2007). Thirty-eight sites (~9%) were variable of 413 total bases in the fragment. Twenty-eight unique COI haplotypes were identified among a total of 81 individuals from 27 caves. Genetic structure among and within sampled populations was evaluated for each geographic partitioning scheme (i.e., hypothesis of structure): (I) across four faunal regions divided by the two major barriers in MP-II, and (II) across 10 minor river drainages to which caves were assigned based on proximity to second- and third-order streams and qualitative topographic information.

Arlequin 3.5 (Excoffier and Lischer 2010) was used to perform hierarchical analyses of molecular variance (AMOVA) for structure hypotheses I and II. Analysis of molecular variance estimates the percentage of genetic variation captured by different pre-defined hierarchical partitions (e.g. among all regions, among caves within each region, and among all caves). From these statistics, fixation indices (F-statistics) were calculated.

FST estimates the degree of differentiation among subpopulations within the total population. The closer FST is to 1, the greater the extent of allelic fixation or identity within populations (Holsinger and Weir 2009). FSC estimates the differentiation among populations within the groups to which they are assigned. The closer FSC is to 1, the more heterogeneity within groups. FCT estimates differentiation among those groups of populations. The closer FCT is to 1, the more divergent the groups are from each other. If strong population genetic structure exists at the group scale being analyzed (i.e., faunal regions), FCT should be high relative to FSC.

Distance matrices and network connections among COI haplotypes were also calculated in Arlequin. Fixation indices (Weir and Cockerham 1984) were calculated from observed diversity within and among populations at each level of geographic structure, and were compared (α = 0.05) to a null resampling distribution of variance components generated from 10,000 permutations in Arlequin.

An unrooted split network based on a NeighborNet algorithm was generated in SplitsTree (Huson and Bryant 2006) to identify distinct genetic clusters from all 81 COI sequences without regard to their relationships. These clusters (identified a posteriori, in contrast to the a priori geographic regions and watersheds in hypotheses I and II) defined the groups for which molecular variance was analyzed for a third structure hypothesis (III).

Network connections among haplotypes were gathered directly from Arlequin output data, and a minimum spanning network of COI haplotypes was constructed using the program HapStar (Teacher and Griffiths 2011). The resulting network was edited in Adobe Illustrator to reflect frequencies of individual haplotypes and their regional associations according to each hypothesis. Mantel tests of association between full matrices and partial submatrices of genetic and geographic distances were performed in R using the package ade4 (Chessel et al. 2004) to detect potential effects of isolation by distance. Mantel tests are commonly performed in studies of population genetics to evaluate the strength of association between genetic and geographic distance (e.g. Diniz-Filho et al. 2013). A high correlation can indicate that some of the population structure observed can be attributed to variation in allele frequencies over geographic distance, which is expected to some degree even in panmictic populations. If a large percentage of genetic variation can be explained by geographic distance, it is difficult to say how much of the observed diversity can be attributed to the particular isolating mechanisms proposed and how much is a consequence of isolation by physical distance (IBD). The population pairwise FST matrix was generated in Arlequin, and the geographic distance matrix was generated from a list of decimal degree coordinates using Geographic Distance Matrix Generator v.1.2.3 (Ersts 2015), an online open source tool provided by the Center for Biodiversity and Conservation, American Museum of Natural History.

Due to the nearly identical external morphology in adults, male genitalia was also examined in a specimen from each cave sampled to see if any differences could be found and if so, to see if there was any correlation between groups discovered via the genetic analysis.


Successful PCR amplification was found to be less reliable for older samples (some as old as five years), despite storage at -80 °C in 95% or stronger ethanol. Despite careful optimization of thermal cycling conditions, agarose gel purification of PCR products was found to considerably improve sequence read quality and was performed for most samples included in the final data set.

The distribution of cave collection sites and proportions of haplotypes from 27 populations are shown in Figs 3 and 4 respectively. Frequencies of COI haplotypes and their proportions in each major faunal region (I) or minor watershed (II) are shown in Fig. 5. A minimum spanning network of COI haplotypes is color coded for each structure hypothesis in Fig. 6A–C). A network of the 81 COI sequences (Fig. 6D) reveals five genetically distinct clusters of structure hypothesis III.

The analysis of molecular variance (AMOVA), from which F-statistics (FST, FCT, and FSC) were calculated to describe nucleotide sequence diversity at hierarchical levels, within and among groups from each hypothesis of structure are summarized in Table 2. The first two hypotheses were based upon a priori geographical hypotheses: (I) the location of caves sampled relative to two zoogeographic barriers proposed by Barr (1985) to be biologically important in MP-II, and (II) the ten minor watersheds to which sampled caves were classified based on assumptions about hydrology gathered from topographic maps (see Table 1).

AMOVA for the a posteriori structure hypothesis III, based on five distinct genetic clusters from a neighbor-joining network of COI sequences produced the greatest difference between FCT and FSC among all three analyses. In other words, when nucleotide diversity is partitioned among hierarchical levels, variance in nucleotide diversity is maximized among groups and minimized within groups. The northernmost 15 sampled populations make up three genetic clusters within an approximately ten-kilometer physical radius of one another. In this arrangement, no haplotypes are shared between the three groups, and the clusters contradict both a priori hypotheses about the locations of important major and minor water barriers to gene flow, especially in the northern part of the MP-II. Mantel tests of group submatrices found population pairwise FST to be independent of geographic distance within each cluster. Among all 15 of these populations, only a maximum of 14% of the observed variation can be explained by geographic distance.

Examination of male genitalia generally showed only slight differences among cave localities examined (Fig. 7). The median lobes were of a consistent shape as were the parameres and internal sac morphology with one notable exception. In specimen 14, the paramere expansion is absent and the internal sac appears to have a different shape within the median lobe. This morphology was found only in Hicksey Cave (abbreviated HIC in all Figs) located in the northern part of the distribution. One should note that paramere expansion is more visible in those specimens that have darker cuticle and hence individuals can appear more different than they actually are due to superficial color differences. No support for a distinct genitalic morphology of D. kentuckensis lexingtoni was observed and the three caves sampled with this subspecies (Great Saltpeter, Teamers, and Mullins Spring Caves) are no more distinct than some of the populations from sets of caves or even single caves such as individuals from Pourover Cave.

Figure 4. 

Distribution of cave collection sites and proportions of haplotypes from 27 populations of Darlingtonea kentuckensis in eastern Kentucky, USA. Circle area corresponds to number of individuals sampled per locality. Different colors indicate different haplotypes; similarity in hue qualitatively indicates sequence similarity. KR: Kentucky River; RR: Rockcastle River; CR: Cumberland River; MVF: Mount Vernon Fault; DD = drainage divide between Kentucky and Rockcastle rivers.

Table 2.

AMOVA statistics, fixation indices, and results of hypothesis tests for structure hypotheses I (four faunal regions), II (ten watersheds), and III (five genetic clusters).

Source of variation Degrees of freedom Sum of squares Variance components Percentage of variation
Among groups 3 149.425 2.19461 (Va) 56.30
Among populations within groups 23 108.884 1.44888 (Vb) 37.17
Within populations 58 14.750 0.25431 (Vc) 6.52
Total 84 273.059 3.89780 100
Fixation Indices: I
FSC 0.85069 Vb and FSC : P(random > observed) = 0.00000***
FST 0.93476 Vc and FST : P(random < observed) = 0.00000***
FCT 0.56304 Va and FCT : P(random > observed) = 0.00000***
Among groups 9 196.197 2.16311 (Va) 60.85
Among populations within groups 17 62.112 1.13762 (Vb) 32.00
Within populations 58 14.750 0.25431 (Vc) 7.15
Total 84 273.059 3.55503 100
Fixation Indices: II
FSC 0.81730 Vb and FSC : P(random > observed) = 0.00000***
FST 0.92846 Vc and FST : P(random < observed) = 0.00000***
FCT 0.60846 Va and FCT : P(random > observed) = 0.00000***
Among groups 4 221.073 3.27840 (Va) 81.93
Among populations within groups 22 37.236 0.46852 (Vb) 11.71
Within populations 58 14.750 0.25431 (Vc) 6.36
Total 84 273.059 4.00124 100
Fixation Indices: III
FSC 0.64818 Vb and FSC : P(random > observed) = 0.00000***
FST 0.93644 Vc and FST : P(random < observed) = 0.00000***
FCT 0.81935 Va and FCT : P(random > observed) = 0.00000***
Figure 5. 

Frequencies of COI haplotypes and their proportions, color coded for each hypothesis of structure; circle area corresponds to number of individuals assigned to each group. Overlain transparent dots show collecting localities. A Four faunal regions of hypothesis I (fifth region unsampled in this study: see discussion and Barr 1985, Kane et al. 1992) B ten minor watersheds of hypothesis II C five genetic clusters of hypothesis III.

Figure 6. 

A–C Minimum spanning networks of COI haplotypes, color-coded for each hypothesis of structure. A Four faunal regions of hypothesis I B ten watersheds of hypothesis II C five genetic clusters of hypothesis III D A split network of 85 COI sequences revealing the five genetically distinct clusters of hypothesis III.

Figure 7. 

Representative male genitalia from 17 of the sampled caves: 1 Wells Cave; 2 Pine Hill Cave; 3–5 Wind Cave; 6 Richardson’s Cave; 7, 8 Lainhart #1 Cave; 9, 10 and 15, 16 Pourover Cave; 11, 12 John Griffin Cave; 13 Climax Cave; 14 Hicksey Cave; 17, 18 Stab Cave; 19, 20 Piney Grove Cave; 21, 22 Dykes Bridge Cave; 23 Great Saltpeter Cave; 24 Teamers Cave; 25 Mullins Spring Cave; 26 Jesse Cave; 27 Steel Hollow Cave. Note that Wells and Dykes Bridge Caves were not included in the genetic study.


FST measures allelic identity within populations, or among-population variation. Across partitioning schemes, FST values close to one indicate that individuals within populations are more similar to each other than to individuals in other populations, corroborating the idea that in general, cave populations in this study are isolated from one another. Structure hypotheses I and II were developed based on a priori information about the locations of cave collection sites relative to (I) two hypothesized major geographic barriers to gene flow or (II) ten watersheds of higher-order streams. Results of AMOVA for evaluating structure hypotheses I and II indicated that for both hypotheses, the majority of total variation (56–61%) is accounted for by variation among the groups defined under each hypothesis. These results support both structure hypotheses I and II over a null hypothesis of panmixia. Due to the similarity of results for both structure hypotheses I and II and because they are not mutually exclusive, neither can be concluded to better represent geographic structure of genetic diversity among the populations sampled. Hence both the major rivers and even some of the smaller watersheds may be geographic barriers to gene flow. High estimates of FSC relative to FCT (Table 2), as well as shared haplotypes among groups in the northern MP-II indicates that neither hypothesis provides the most optimal scheme for partitioning the observed genetic diversity. The lack of robust support for a partitioning scheme based on small watersheds is not necessarily evidence against the influence of climate cycles on the process of lineage diversification. Many caves do not “belong” to a single watershed, but rather may connect or fall between two or more. This factor, along with the uncertainty surrounding cave connectivity via small passages accessible only by small taxa like these beetles, can make it difficult to truly know the possible connectivity of some caves to one watershed over another. Additionally, it is possible that the shape, size, and the pathway of the watersheds in this area changed throughout the recent Pleistocene and earlier. Hence the separation of populations by hypothesized barriers between caves assigned to different watersheds may have resulted from actual watershed barriers, intra-karst heterogeneity, and or climate cycles at various times that in turn helped drive or prevent cave colonization.

Structure hypothesis III was developed based on the five genetic clusters resulting from a split network. The boundaries for the five population clusters in this hypothesis were determined solely by clustering based on genetic distances among sequences, independently of any a priori geographic information. AMOVA statistics for structure hypothesis III (Table 2) indicate that for each hypothesis of structure, among-group variation accounts for a higher percentage of the total variation than within-group variation. These results support all three hypotheses as better models for structured diversity compared with a null model of panmixia. However, variation among groups (genetic clusters) in hypothesis III accounts for much more of the total variation (82%) than either hypothesis I or II (56% and 61%, respectively). Further, only in structure hypothesis III does diversity among groups (FCT = 0.82) exceed diversity within groups (FSC = 0.65). Unlike hypotheses I and II, no haplotypes are shared between the five clusters. Lastly, these five genetic clusters form natural, geographically proximate groupings. Hence the evidence supports hypothesis III as the most representative model for the geographic structure of genetic diversity among sampled populations, and especially for those in the northern MP-II part of the distribution.

If geographic distance is strongly positively correlated with genetic distance, gaps in sampling (rather than specific geographic features acting as barriers to gene flow) could be responsible for at least some of the observed clustering of populations. Results of partial Mantel tests (Table 3) indicate up to 18% of the total observed genetic variation across all 27 populations can be attributed solely to the influence of geographic distance. Across the 15 northern populations (three of the five genetic clusters), IBD could explain up to 14% of the total variation. However, low Mantel correlations for population subsets corresponding to each of these three clusters suggests that the genetic structure observed in this region (Rockcastle, Jackson, and Estill counties) is most likely due to actual barriers to gene flow and not simply isolation by distance.

Table 3.

Results of Mantel tests (10000 permutations) of association between geographic distance and population pairwise FST within and among groups from hypotheses I and III, containing the same 15 northern MP-II populations partitioned in different ways.

Hypothesis (group #) Populations included % variation explained by geographic distance Pobs>sim(α=0.05)
III (1) CLF, HIC, LAI, MOR, HIS, SRI, JGR, LAK, CLI <1 0.468
III (2) FLE, PHC, SOR <1 0.6637
III (5) TEA, MUL, GSP <1 0.673
I (1) CLF, HIC, LAI, MOR, HIS 7 0.761
I (2) SRI, JGR, LAK, CLI, MUL, GSP, TEA, SOR, PHC, FLE 19 0.0035
all northern MP-II CLF, HIC, LAI, MOR, HIS, SRI, JGR, LAK, CLI, FLE, PHC, SOR, TEA, MUL, GSP 14 0.0033
all 27 populations 18 0.0001

Barr (1985) recognized that the fragmented geology of Rockcastle County, Kentucky may account for the morphological (and genetic) variability in the region, which is topographically complex and dissected with many rivers and streams. The five clusters (including two completely outside Rockcastle County) could represent distinct lineages important in considering the ecology and evolution of Darlingtonea, but divergence times and particular geographic or geologic features consistent with the apparent locations of most putative isolating barriers have not yet been investigated systematically; only the Mount Vernon fault has been well studied.

The Mount Vernon fault (Fig. 4) runs through a cave-rich area of Rockcastle County, Kentucky. Based on its position in the otherwise relatively less faulted MP-II compared to other karst formations (KGS 2017), it may serve as a stratigraphic barrier isolating one of the three northern clusters (D. kentuckensis lexingtoni populations) from the other two (Fig. 5C red colored pie #1 and Fig. 6C). The relatively cave-poor divide between the Kentucky and Rockcastle river drainages (KGS 2017), hypothesized by Barr (1985) to represent an important stratigraphic barrier, is not supported in this study given that populations of the northernmost genetic cluster fall on both sides of the barrier. The influence of the three-way fluvial barrier proposed by Barr (1985), formed by the confluence of the Cumberland River and its Big South Fork, is not explicitly supported but cannot be ruled out due to lack of breadth and spatial resolution in population sampling. Examination of geographically proximate populations in each sector of this “river triangle” (Barr 1985, see fig.1b in Kane et al. 1992) would help to clarify its role as an isolating barrier. Though our study did not explicitly test the effect of meander frequency (Barr 1985) on terrestrial troglobiont dispersal potential, the distributional patterns observed (Fig. 4) do not conflict with the hypothesis that smaller, meandering waterways are less likely or even unlikely to act as dispersal barriers compared to large rivers.

The sampling scheme of our study makes it difficult to extricate signal due to population structure from that due to IBD for the two genetic clusters on either side of the Cumberland River, which are strongly clustered spatially (Fig. 4). An ideal scheme would evenly sample many population pairs on either side of and at increasing distances from each proposed barrier. Under this sampling regime, results of partial Mantel tests within and among groups separated by each proposed barrier could be used to detect population structure amid underlying “noise” from IBD. Isolation between groups across fluvial barriers with different calculated meander frequencies could also be formally compared. Such a systematic sampling method would be challenging for this group of organisms however, as caves are unevenly distributed across the landscape and access restrictions further reduce the number of available cave sampling localities.

Overall, the limits of neither major nor minor watersheds alone adequately model the observed distribution of genetic diversity across sampled populations of D. kentuckensis. Geographic distance and landscape features, both stratigraphic and fluvial, appear to have each contributed to this distribution. Determination of the boundaries of cryptic species or subspecies, inference of their pattern of relatedness, and identification of predictive characteristics of isolating barriers will require further sampling of additional populations and more complete and/or additional molecular loci.


Based on CO1 data alone, there is a wide range of divergence values between taxa that can be defined as separate species on their own evolutionary trajectory from other lineages (Hebert et al. 2003). No formal taxonomic changes are proposed herein as a result of this study, as full or nascent species could be represented by all, some, or none of the five genetic clusters discovered among twenty-seven sampled populations of D. kentuckensis, depending upon the species definition favored. Both genetic and some morphological evidence supports the hypothesis that D. kentuckensis consists of isolated populations that could be considered as separate cryptic species or perhaps subspecies. Hebert et al. (2003) gives an average sequence divergence of 11.2% between species of beetles within the same genus, but divergence ranges from below 1% to 16–32% depending upon the paired taxa examined. Genetic divergence between each of the five populations studied herein differ by ~1.3%, a percentage that is within the range of CO1 sequence divergence between species, although it is certainly on the low side. Regardless, populations within the range of the subspecies D. kentuckensis lexingtoni do form a genetically distinct cluster that is especially supported by this study; additionally all three northernmost clusters are geographically proximate but genetically distinct, with little evidence that isolation by distance is an influence on the pattern of genetic structure. The observed strong correlation between pairwise FST and geographic distance among the two southern populations may either be an artifact of sampling deficiency that overlooks intermediate haplotypes or a reflection of a real historical sequence of colonization events. Therefore, these results can be viewed as a starting point for continued investigation, using additional molecular markers and denser sampling, of the historical phylogeography and species limits in this group and other related taxa.


This research was supported by funds from TKP and JRJ, in addition to internal funding from Western Kentucky University through a Graduate Student Research Grant and a Research and Creative Activities Program (RCAP) grant. Naomi Rowland provided valuable advice and assistance with molecular work. We thank Dr. Matthew Niemiller for his general guidance and valuable comments on the manuscript, Dr. Karen Ober for her contribution of specimens, data, and advice, and Elise Valkanas for her phylogenetic studies on eastern North American cave beetles (supported by an NSF- REU grant) which inspired and contributed to this project. We thank Den and Sheila Roenfeldt for supporting a travel fellowship which allowed OFB to present this research and gather feedback from colleagues at the International Society for Subterranean Biology’s 2016 conference. Our appreciation to Arnaud Faille and one anonymous reviewer whose insightful comments improved the manuscript. We thank the many landowners who granted permission for us to sample caves on private property, as well as Jim Currens, Julian Lewis, Ben Miller, Jason Polk, Lee Florea, Kurt Helf, The Rockcastle Karst Conservancy, the Green River Grotto, John Andersland, Rob Neidlinger, and the Kentucky Speleological Society.


  • Barr TC (1969) Evolution of the (Coleoptera) Carabidae in the Southern Appalachians. In: Holt PC, Roane MK, Parker BC (Eds) The distributional history of the biota of the Southern Appalachians. Virginia Polytechnic Institute and State University, Blacksburg, 67–92.
  • Barr TC (1971) Trechoblemus in North America, with a key to North American genera of Trechinae (Coleoptera: Carabidae). Psyche 78: 140–149.
  • Barr TC (1973) Refugees of the ice age. Natural History 82 5): 26–35, 72–73.
  • Barr TC (1979) The taxonomy, distribution, and affinities of Neaphaenops, with notes on associated species of Pseudanophthalmus (Coleoptera, Carabidae). American Museum novitates: no. 2682: 1–20.
  • Barr TC (1981) Pseudanophthalmus from Appalachian caves (Coleoptera: Carabidae): the engelhardti complex. Brimleyana 5: 37–94.
  • Barr TC (1985) Pattern and process in speciation of trechine beetles in eastern North America (Coleoptera: Carabidae: Trechinae). In: Ball GE (Ed.) Taxonomy, Phylogeny and Zoogeography of Beetles and Ants. Dr W. Junk, Dordrecht, 350–407.
  • Barr TC (2004) A classification and checklist of the genus Pseudanophthalmus Jeannel (Coleoptera: Carabidae: Trechinae). Virginia Museum of Natural History Special Publication 11: 1–52.
  • Chessel D, Dufour AB, Thioulouse J (2004) The ade4 package – I: One-table methods. R News 4: 5–10.
  • Christman MC, Doctor DH, Niemiller ML, Weary DJ, Young JA, Zigler KS, et al. (2016) Predicting the Occurrence of Cave-Inhabiting Fauna Based on Features of the Earth Surface Environment. PLoS ONE 11(8): e0160408.
  • Currens JC (2002) Kentucky is karst country! What you should know about sinkholes and springs. Information Circular 4, Series XII, University of Kentucky, Kentucky Geological Survey, 35 pp.
  • Derkarabetian S, Steinmann DB, Hedin M (2010) Repeated and time-correlated morphological convergence in cave-dwelling harvestmen (Opiliones, Laniatores) from montane Western North America. PLoS One 5(5): e10388.
  • Diniz-Filho JA, Soares TN, Lima JS, Dobrovolski R, Landeiro VL, de Campos Telles MP, Rangel TF, Bini LM (2013) Mantel test in population genetics. Genetics and Molecular Biology 36(4): 475–85.
  • Ersts PJ (2015) Geographic Distance Matrix Generator (version 1.2.3). American Museum of Natural History, Center for Biodiversity and Conservation.
  • Faille A, Rene Tänzler R, Toussaint EFA (2015) On the way to speciation: Shedding light on the karstic phylogeography of the micro-endemic cave beetle Aphaenops cerberus in the Pyrenees. Journal of Heredity 106(6): 692–699.
  • Goudarzi F, Hemami M, Rancilhac L, et al. (2019) Geographic separation and genetic differentiation of populations are not coupled with niche differentiation in threatened Kaiser’s spotted newt (Neurergus kaiseri). Scientific Reports 9: 6239.
  • Hanski I (1999) Metapopulation Ecology. Oxford University Press.
  • Harker DF, Barr TC (1979) Caves and associated fauna of eastern Kentucky. Eastern Kentucky Coal Field: Preliminary Investigations of natural features and cultural resources. Vol. 3. Kentucky Nature Preserves Commission, Frankfort, 130 pp.
  • Hebert PD, Ratnasingham S, de Waard JR (2003) Barcoding animal life: cytochrome c oxidase subunit 1 divergences among closely related species. Proceedings of the Royal Society of London. Series B: Biological Sciences, 270(Suppl 1): S96–S99.
  • Hedin M, Thomas SM (2010) Molecular systematics of eastern North American Phalangodidae (Arachnida: Opiliones: Laniatores), demonstrating convergent morphological evolution in caves. Molecular Phylogenetics and Evolution 54: 107–121.
  • Holsinger JR (1988) Troglobites: The Evolution of Cave-Dwelling Organisms. American Scientist 76(2): 146–153.
  • Holsinger KE, Weir BS (2009) Genetics in geographically structured populations: defining, estimating and interpreting F(ST). Nature Reviews: Genetics 10: 639–650.
  • Katz AD, Taylor SJ, Davis MA (2018) At the confluence of vicariance and dispersal: Phylogeography of cavernicolous springtails (Collembola: Arrhopalitidae, Tomoceridae) codistributed across a geologically complex karst landscape in Illinois and Missouri. Ecology and Evolution 2018: 10306–10325.
  • Kearse M, Moir R, Wilson A, Stones-Havas S, Cheung M, Sturrock S, Buxton S, Cooper A, Markowitz S, Duran C, Thierer T, Ashton B, Meintjes P, Drummond A (2012) Geneious Basic: An integrated and extendable desktop software platform for the organization and analysis of sequence data. Bioinformatics 28(12): 1647–1649.
  • Kentucky Geological Survey (2017) Kentucky Geologic Map Information Service, University of Kentucky.
  • Larkin MA, Blackshields G, Brown NP, Chenna R, McGettigan PA, McWilliam H, Valentin F, Wallace IM, Wilm A, Lopez R, Thompson JD, Gibson TJ, Higgins DG (2007) Clustal W and Clustal X version 2.0. Bioinformatics 23(21): 2947–2948.
  • Lewis JJ, Lewis SL (2005) Cave fauna study for the Interstate 66 EIS (Somerset to London, Kentucky). Proceedings of the 2005 National Cave and Karst Management Symposium: 15–20.
  • Maddison DR, Kanda K, Boyd OF, Faille A, Porch N, Erwin TL, Roig-Juñent S (2019) Phylogeny of the beetle supertribe Trechitae (Coleoptera: Carabidae): Unexpected clades, isolated lineages, and morphological convergence. Molecular Phylogenetics and Evolution 132: 151–176.
  • Niemiller ML, Fitzpatrick BM, Miller BT (2008) Recent divergence with gene flow in Tennessee cave salamanders (Plethodontidae: Gyrinophilus) inferred from gene genealogies. Molecular Ecology 17(9): 2258–2275.
  • Niemiller ML, Near TJ, Fitzpatrick BM (2012) Delimiting species using multilocus data: diagnosing cryptic diversity in the southern cavefish Typhlichthys subterraneus (Teleostei: Amblyopsidae). Evolution 66, 846–866.
  • Niemiller ML, McCandless JR, Reynolds RG, Caddle J, Tillquist CR, Near TJ, Pearson WD, Fitzpatrick BM (2013) Effects of climatic and geological processes during the Pleistocene on the evolutionary history of the northern cavefish, Amblyopsis spelaea (Teleostei: Amblyopsidae). Evolution 67: 1011–1025.
  • Peck SB (1998) A summary of diversity and distribution of the obligate cave-inhabiting faunas of the United States and Canada. Journal of Cave and Karst Studies 60(1): 18–26.
  • Simon C, Frati F, Beckenbach A, Crespi B, Liu H, Flook P (1994) Evolution, weighting, and phylogenetic utility of mitochondrial gene sequences and a compilation of conserved polymerase chain reaction primers. Annals of the Entomological Society of America 87(6): 651–701.
  • Valentine JM (1952) New genera of anophthalmid beetles from cumberland caves (Carabidae, Trechini). Geological Survey of Alabama, Museum Paper 34, 41 pp.
  • Wiens JJ, Chippindale PT, Hillis DM (2003) When are phylogenetic analyses misled by convergence? A case study in Texas cave salamanders. Systematic Biology 52: 501–514.
  • White WB, Culver DC, Pipan T (2019) Encyclopedia of Caves, 3rd Edition, Academic Press, 1250 pp.