Morphometrics and phylogeography of the cave-obligate land snail Helicodiscus barri (Gastropoda, Stylommatophora, Helicodiscidae)
expand article infoNicholas S. Gladstone, Matthew L. Niemiller§, Evelyn B. Pieper, Katherine E. Dooley§, Michael L. McKinney
‡ University of Tennessee, Knoxville, United States of America
§ The University of Alabama in Huntsville, Huntsville, United States of America
Open Access


Molecular studies have recently led to the detection of many cryptic species complexes within morphologically ambiguous species formerly undescribed by the scientific community. Organisms such as land snails are at a particularly higher risk of species misidentification and misinterpretation, in that gastropod systematics are based almost entirely on external shell morphology. Subterranean ecosystems are associated with especially high degrees of cryptic speciation, largely owing to the abiotic similarities of these systems. In this study, we attempt to diagnose the potential cryptic diversity in the troglobitic land snail Helicodiscus barri. Land snails are generally associated with having low vagility, and as such this species’ broad, mosaic distribution indicates the misdiagnosis of this organism as a single species. We analyze both mitochondrial (16S, CO1) and nuclear (28S, H3) genetic data for 23 populations. Phylogeny for H. barri was reconstructed using both maximum-likelihood and Bayesian approaches to assess relationships among populations, and two species delimitation methods (mPTP and ABGD) were used to detect the presence of unique molecular operational taxonomic units (MOTUs). Species delimitation results revealed seven and sixteen MOTUs respectively, suggesting the presence of several cryptic lineages within H. barri. To assess how external shell morphology corresponds with patterns of genetic and environmental variation, two morphometric approaches were used incorporating 115 shells from 31 populations. Both morphometric approaches reveal a significant environmental influence on shell morphology, and one approach showed the significance of MOTU groups. We discuss the delimitation and morphometric results and additionally provide discussion on the taxonomic and conservation implications of this study.


Helicodiscidae, subterranean ecology, morphometrics, MOTUs, cryptic species


Caves provide a model system for studying the evolutionary processes and historical factors related to biogeography and speciation (Juan et al. 2010). Cave systems, characterized by geographic isolation and relatively simple biological communities, often are viewed as analogous to oceanic islands (Culver and Pipan 2009, Snowman et al. 2010). Strong selective pressures and the isolation of subterranean ecosystems can result in morphological stasis among otherwise genetically distinct species, largely due to the parallel or convergent evolution of these lineages (Lefébure et al. 2006, Finston et al. 2007, Niemiller et al. 2012). Further, many troglobites (i.e., terrestrial cave-obligates) exhibit broad, mosaic distribution patterns which, in conjunction with morphological stasis, often confound traditional approaches of delimitating species boundaries (Jochum et al. 2015). Consequently, troglobites are ideal models to address fundamental questions in ecology and evolution and provide a platform to approach a more modernized integration of taxonomic methods.

An increasing number of studies has examined population genetic and phylogeographic hypotheses of subterranean fauna (e.g., Moulds et al. 2007, Snowman et al. 2010, Weckstein et al. 2016), which have greatly increased our understanding of colonization history, speciation, dispersal, and biogeography of troglobitic taxa (Juan et al. 2010). Additional phylogeographic studies have uncovered considerable levels of cryptic diversity in subterranean species (Finston et al. 2007, Juan and Emerson 2010, Niemiller et al. 2012). Due to increasing advances in imaging technology, studies that incorporate morphometric analyses often complement such molecular findings (Jochum et al. 2015, Armbruster et al. 2016, Burress et al. 2017, Inäbnit et al. 2019). The misidentification of species can hinder assessments of biodiversity and conservation of cryptic species. Therefore, an integrative taxonomic evaluation of troglobitic taxa is needed to fully assess species richness within these systems, and to better inform their respective evolutionary histories. Moreover, cryptic species complexes may be comprised of groups already at significant risk of extinction (Niemiller et al. 2013).

Land snails (Phylum Mollusca, Class Gastropoda) are a species-rich group, with over 24,000 currently recognized species and over 35,000 species thought to exist globally (Barker 2001, Lydeard et al. 2004). The land snail fauna of eastern North America is exceptionally diverse, with over 500 documented species (Hubricht 1985, Nekola 2014). However, this likely represents an underestimate of total species richness in this region. Larger species are often associated with mesic forest ecosystems with high levels of moisture, leaf litter, and calcium (Goodfriend 1986, Pearce and Örstan 2006). Yet, land snails utilize a variety of microhabitats often neglected in sampling efforts within these areas (Cameron and Pokryszko 2005). Further, land snails occur at high density in karst-rich landscapes, and subterranean habitats are particularly under-sampled within the region (Clements et al. 2008, Niemiller and Zigler 2013).

Nearly 75% of all land snails in eastern North America are considered terrestrial micromolluscs (< 5 mm) and comprise a significant portion of all land snail diversity (Nekola 2005, Liew et al. 2008). Many of these species tend to be particularly undersampled and often require the collection of soil and leaf litter samples to discover them (Liew et al. 2008, Nekola and Coles 2010, Durkan et al. 2013). Regions hypothesized to have higher levels of snail biodiversity have had varying and potentially insufficient sampling effort, with many species remaining undescribed (Dourson 2007, Douglas et al. 2014, Dinkins and Dinkins 2018). Moreover, there is a paucity of studies examining intraspecific morphological variation in micromolluscs, obscuring accurate geographic ranges for these species (Nekola and Coles 2010). Thus, because most land snail species are delimited based on conchology (i.e., shell variation), a high incidence of misidentification of minute species occurs in many natural history collections (Hubricht 1985, Nekola and Coles 2010). The continued misidentification of species can have significant impacts on biodiversity assessments and conservation management (Bickford et al. 2007).

Strictly employing morphological data to delimit extant species in the genomic era is often met with criticism (Hermsen and Hendricks 2008, Duminil and Di Michele 2009, Carstens et al. 2013). An integrative taxonomy, i.e., a combination of morphological, ecological, and genetic data when considering phylogenetic relationships, is necessary to facilitate proper interpretations of biological patterns (Dayrat 2005, Weigand et al. 2012, Inäbnit et al. 2019). For gastropods, there are few discrete shell characters that can be used in phylogenetic hypotheses, and conchology is highly variable in response to environmental factors and other selective pressures (Goodfriend 1986, Smith and Hendricks 2013). However, morphometric analyses can contribute to species hypotheses when combined with genetic data (Hermsen and Hendricks 2008, Miller 2016, Inäbnit et al. 2019). Moreover, applying morphometric analyses can inform the causal mechanisms for shape variation between gastropod populations (Vergara et al. 2017).

Terrestrial micromolluscs of the genus Helicodiscus Morse, 1864 are found throughout the eastern United States (Hubricht 1985). This genus is known for its unique conchological sculpture, often exhibiting depigmented soft bodies and prominent spiraling striae on the shells of both surface and subterranean species. Many of these species are calciphiles, and two species – H. barri Hubricht, 1962 and H. notius specus Hubricht, 1962 – have even adopted a cave-obligate existence (Hubricht 1962). The distributions of these troglobites span both the Interior Low Plateau (ILP) and the Appalachians karst regions, covering multiple physiographic provinces within their ranges. The latter species is only known from six caves in Kentucky, Tennessee, and Virginia, whereas the former is known from 49 caves in Tennessee, Alabama, and Georgia. Two additional Helicodiscus species that were previously thought to be troglobitic – H. hadenoecus Hubricht, 1962 and H. punctatellus Morrison, 1942 – have been discovered at surface localities widely disjunct from their otherwise subterranean distribution (Coney et al. 1982; Hotopp et al. 2013). These distribution patterns suggest the potential for cryptic diversity among subterranean taxa within this genus. Morphological stasis is highly prevelant in troglobites despite significant genetic divergence, and, therefore, the mosaic distributions of these snails warrant investigation (Juan and Emerson 2010, Weigand et al. 2012).

Here, we conduct the first study examining morphological variation and phylogeography of the cave-obligate land snail Helicodiscus barri. Recent cave bioinventory efforts within the ILP and Appalachians karst regions have yielded several additional specimens of this species for comparison across multiple physiographic provinces. The disjunct, mosaic distribution pattern of H. barri in conjunction with a lack of clear morphological variation is consistent with a high potential for cryptic diversity, as observed in other subterranean taxa (Snowman et al. 2010, Loria et al. 2011, Niemiller et al. 2012, Inäbnit et al. 2019). We examined museum accessions of H. barri while sampling caves within the states of Tennessee, Alabama, and Georgia for additional specimens. Using phylogeographic approaches, we (1) assessed patterns of genetic variation of H. barri; (2) employed two species delimitation methods (ABGD and mPTP) to infer the presence of cryptic lineages; and (3) tested if current species hypotheses based on conchology correspond with patterns of genetic variation. Further, we assessed morphological variation between H. barri populations using traditional morphometrics (TM) and landmark-based geometric morphometrics (GM).


Specimen collection

Shell specimens were collected from 31 populations of cave-dwelling Helicodiscus barri from the dark zone of caves within both the ILP and Appalachian karst regions in Tennessee and Alabama (115 total individual specimens collected). Each survey typically involved two to four researchers (maximum 12), with a search effort of two to 36 person-hours per cave visit. In total, 74 caves were visited from 13 March 2013 to 19 June 2018 by NSG, totaling ca. 300 person-hours. Snail specimens were preserved in 100% ethanol and identified using published keys and species descriptions (Pilsbry 1948, Hubricht 1962, Dourson 2010), as well as examination by taxonomic specialists (Dan and Judy Dourson). Specimens from ten additional populations were provided by the Field Museum of Natural History (FMNH), Florida Museum of Natural History (FLMNH), and the Auburn Museum of Natural History (AUM). In total, 154 shells were examined (see Table 1). The geographic distribution of populations utilized within this study can be found in Figure 1.

Helicodiscus barri populations incorporated in this study, including 17 new populations. Cave names, Tennessee Cave Survey (TCS) cave number, county and state are provided, as well as information regarding which populations were considered in morphometric and genetic analyses.

Sample Cave TCS No. County State References n Morphology Genetic
NSG-DI3 Bowman Cave TDI3 Dickson TN This study 2 X X
NSG-KN112 Brents Cave TKN112 Knox TN This study 3 X X
AUM28348 Bull Run Cave TDA4 Davidson TN Hubricht 1964 2 X X
MLN 14-054.12; NSG-JK3 Carter Cave TJK3 Jackson TN Gladstone et al. 2018 5 X X
NSG-VB547 Cave Between the Caves TVB547 Van Buren TN Lewis 2005 6 X
NSG-RN5 Cave Creek Cave TRN5 Roane TN This study 1 X
MLN 14-007 Christmas Cave TDK72 DeKalb TN Gladstone et al. 2018 1 X X
MLN 13-056 Clarksville Lake Cave TMY11 Montgomery TN Gladstone et al. 2018 2 X X
FMNH239117 Collier Cave ALD100 Lauderdale AL Peck 1989 4 X
FMNH239122; NSG-DI6 Columbia Caverns TDI6 Dickson TN Hubricht 1962 7 X X
NSG-KN50 Conner Creek Cave TKN50 Knox TN This study 5 X X
FMNH239121 Culbertson Cave TUN22 Union TN Hubricht 1985 1 X
NSG-AN5 Demarcus Cave TAN5 Anderson TN This study 3 X
MLN 14-015.3 Dry Cave TFR9 Franklin TN Gladstone et al. 2018 2 X
AUM27534-T2 Frazier Hollow Cave DK11 DeKalb TN This study 1 X X
NSG-DI27 East Fork Cave TDI27 Dickson TN This study 2 X
AUM28173 Hering Cave Madison AL This study 1 X X
NSG-FR14 Keith Cave TFR14 Franklin TN Lewis 2005 10 X X
UF 405128 Lady Finger Bluff Trail Perry TN Gladstone et al. 2018 1 X
WC13-165 Lovelady Cave THM56 Hamilton TN This study 1
NSG-MM10 McCorkle Cave TMM10 McMinn TN This study 1
NSG-VB9 McCoy Cave TVB9 Van Buren TN This study 2 X X
AUM27855 New Salem Cave Nr1 TSM10 Smith TN This study 2 X X
MLN 15-007.9 Oaks Cave TUN5 Union TN Gladstone et al. 2018 1 X X
FMNH305126; NSG-AN12 Offut Cave TAN12 Anderson TN Hubricht 1985 8 X X
MLN 15-006.19; NSG-CM8 Panther Cave No. 1 TCM8 Campbell TN Gladstone et al. 2018 7 X X
FMNH239120 Parkers Cave GKH119 Chattooga GA Holsinger and Peck (1971) 2 X
NSG-KN108 Pedigo Cave Nr. 2 TKN108 Knox TN This study 1 X
NSG-AN6 Robert Smith Cave TAN6 Anderson TN This study 2 X
MLN 13-000 Rockhouse Cave ALM312 Limestone AL Gladstone et al. 2018 1 X
AUM27652 Rogers Hollow Cave TUN23 Union TN This study 2 X
FMNH239118 Shelta Cave AMD4 Madison AL Peck 1989 3 X
NSG-OV440; GC1 Slippery Slit Cave TOV440 Overton TN Lewis 2005 3 X X
KSZ15-313 Smartt Farm Cave GWK124 Walker GA This study 1
NSG-VB657 Swamp River Cave TVB657 Van Buren TN This study 1 X
MLN-16.0228 Weavers Cave TAN22 Anderson TN Gladstone et al. 2018 2 X X
NSG-KN80 Wilke Waller Cave TKN80 Knox TN This study 1
Figure 1. 

Geographic distribution of Helicodiscus barri from this study in relation to karst adapted from Weary and Doctor (2014). Triangles represent cave populations.

DNA extraction, amplification, and sequencing

Genomic DNA was obtained from soft tissue of each live specimen collected. The shells of smaller individuals were removed prior to DNA extraction. Tissue was removed from larger shells by breaking a small opening into the abapertural side of the shell or the shell base, so that the shell was not completely destroyed and remained identifiable. Each DNA extraction was performed using the Qiagen® DNeasy Blood and Tissue kit following the manufacturer’s protocol (Qiagen Sciences, Louisville, KY). Polymerase chain reaction (PCR) was used to amplify fragments of the mitochondrial (mt) 16S ribosomal RNA locus using the primer pair 16Sa/16Sb (Palumbi et al. 1991), mt cytochrome oxidase subunit 1 (COI) locus using the primer pair LCOI490/ HCO2198 (Folmer et al. 1994), nuclear 28S ribosomal RNA locus using the primer pair 28Sna1/28Sna2 (Kano et al. 2002), and nuclear histone 3 (H3) locus using the primer pair H3F/H3R (Colgan et al. 2000). PCR products were purified using ExoSAP-IT (Affymetrix) and sequenced in both directions with BigDye chemistry at Eurofins MWG Operon (Louisville, KY).

Genetic analyses

Forward and reverse sequences were assembled into contigs and edited in Sequencher v.5.1 (Gene Codes Corporation, Ann Arbor, MI). Alignments were modified by the manual trimming of the 3’ and 5’ primer ends. Ambiguous base calls and double peaks within heterozygotes were assessed visually with the chromatograms. Sequences were then aligned using MUSCLE under default parameters implemented in MEGA X v.10.0.5 (Kumar et al. 2018). All sequence data generated from this study was accessioned into GenBank (see Suppl. material 1). PartitionFinder v.2.1.1 (Lanfear et al. 2012) was used to determine the best model of sequence evolution for each partition based on the Bayesian information criterion (BIC). A general time-reversible model of sequence evolution with corrections for a discrete gamma distribution and a proportion of invariant sites (GTR+Γ+I) was chosen for 16S. The Hasegawa et al. (1985) model (HKY) with corrections for a discrete gamma distribution was chosen for the first and second codon positions of both CO1 and H3 as well as for 28S. A symmetrical model with corrections for a discrete gamma distribution (SYM+ Γ) was chosen for the third codon position of CO1 and H3 (Zharkikh 1994). Due to uneven coverage of genetic data across specimens, three unique H. barri datasets were assessed: CO1, mtDNA (CO1 + 16S), and mtDNA + nDNA (CO1 + 16S + 28S + H3). Discus rotundatus was used as an outgroup for all phylogenetic analyses. Summary statistics of the H. barri molecular dataset including haplotype and nucleotide diversity, number of segregating sites, haplotypes, and mutations were calculated in DnaSP v.6.12.01 (Librado and Rozas 2009). Uncorrected p-distances within and between cave populations were used as a metric of genetic divergence and calculated in MEGA X v.10.0.5 (Kumar et al. 2018). A haplotype network for all specimens for which genetic data were available was created in SplitsTree v.4.14.8 (Huson and Bryant 2005) using the NeighborNet network method with uncorrected p-distances.

Phylogenetic analyses and species delimitation

Phylogenetic trees were inferred utilizing both a Maximum-Likelihood (ML) and Bayesian-inference (BI) approach. ML analyses were conducted using RAxML v.8.0 (Stamatakis 2014) as implemented through the T-REX web server (Boc et al. 2012). A consensus tree was generated from the CO1, mtDNA, and mtDNA + nDNA datasets using rapid bootstraps for 100,000 replicates under a GTR+Γ+I model of evolution. The BI analyses were conducted in MrBayes v.3.2.6 (Ronquist et al. 2012) using a random start tree with three heated and one cold chain (default temperature of 0.1). This was run twice for 50,000,000 generations and sampled every 1,000 generations under the models of evolution determined by PartitionFinder. The first 25% of samples (12,500,000) were discarded as burn-in. Convergence of runs was assessed utilizing Tracer v. 1.4 (Rambaut and Drummond 2007).

The generation of molecular barcodes is often utilized in species delimitation in understudied groups (or in this case, those that are morphologically ambiguous; Pons et al. 2006; Rubinoff 2006; Weigand et al. 2012, 2014). As such, two species delimitation approaches were subsequently used in the identification of Molecular Operational Taxonomic Units (MOTUs; Floyd et al. 2002): 1) Automatic Barcode Gap Recovery (ABGD; Puillandre et al. 2012), and 2) Multi-rate Poisson Tree Processes (mPTP; Kapli et al. 2017). ABGD partitions samples into candidate species based on a statistically inferred barcode gap. The barcoding gap is defined as a notable disparity between pairwise genetic distances, presumably between intraspecific and interspecific distances. This process is applied recursively to newly obtained groupings of sequences, to assess the possibility of internal division. This method was employed on the CO1 dataset excluding the outgroup (n = 24) via the ABGD web server ( using the Kimura two-parameter (K2P; Kimura 1980) model with a standard X (relative gap width) = 1.5.

The initial development of PTP models assumed one exponential distribution for speciation events and one for all coalescent events (Zhang et al. 2013). In contrast, the mPTP approach fits speciation events for each candidate species to a unique exponential distribution, greatly improving the quality of results (Kapli et al. 2017). This method requires a rooted phylogenetic tree and partitions samples into candidate species based upon the number of substitutions under assumed Poisson processes. Intraspecific substitution rates should be notably smaller than interspecific rates. This method does not require an ultrametric tree, which is ideal given little reliable fossil data for Helicodiscidae and the variability of molecular clock rates in Stylommatophoran gastropods (Thomaz et al. 1996, Chiba 1999, Van Riel et al. 2005). A rooted tree was generated for the CO1 dataset using the methods previously outlined for RAxML under the models of evolution determined by PartitionFinder. Analysis was carried out on the mPTP webserver ( for the maximum 100,000 MCMC generations, with 25% of samples (25,000) conservatively discarded as burn-in.

Morphometric analyses

Specimens were photographed using a Canon 6D digital SLR camera mounted on the Macropod PRO Micro Kit (Macroscopic Solutions, Tolland, CT). Each shell was photographed using a Canon MP-E 65mm f/2.8 1–5× macro lens in three views: ventral, dorsal and apertural. MacroMagnification settings were extracted from the images using ExifTool v. Images were imported to Adobe Photoshop CS5 Extended v.12.1 and were subsequently scaled. Reproductive anatomy was not evaluated, given that most specimens were taken from museum collections with no soft tissue available. Two morphometrics methods were employed: geometric morphometrics (GM) and traditional morphometrics (TM).

GM techniques allow for the quantification and assessment of morphological variation. Biologically-meaningful landmarks (LMs) and semilandmarks (SLMs) of the Helicodiscus shells were digitized using tpsDig2 v. 2.32 (Rohlf 2015, available at Nine homologous static LM were placed across each specimen. LM 1 is Type 1, characterizing the discrete juxtaposition of the homologous shell structure. LM 4 and LM 5 are Type 2, characterizing geometric maxima of curvature. All remaining LMs were Type 3, characterizing more than one region of each shell (Bookstein 1997). These nine LMs were combined with two manually traced curves of three equidistant SLMs anchored on LM 4–5 and LM 5–6 (see Figure 2A). Appending tps curves to SLM was achieved using tpsUtil v. 1.76 (Rohlf 2015). This results in a total of nine fixed and six semi-landmarks.

Figure 2. 

A Landmark scheme for geomorphometric analyses. Red circles represented landmarks (LM), blue circles represent semi-landmarks (SLM). B Shell measurements utilized for the traditional morphometric (TM) analyses.

To eliminate variation due to orientation of the shell or size, a Procrustes superimposition was performed using the geomorph package v.2.0. in RStudio v 1.1.456 with R v. 3.5 (Adams and Otarola-Castillo 2013, Adams et al. 2014, RC Team 2014). These data were subjected to principal component analysis (PCA) to evaluate the distribution of populations in morphospace. Several alternative landmark schemes were tested but provided no notable differences in PCA results. The most conservative approach was employed to reduce the number of variables introduced into the downstream analyses. Correlation coefficients (PC loadings) of individual variables were assessed visually to determine which specific variables are significant to each PC as to interpret what shell characteristics account for variability of the dataset. To quantify error associated with landmark placement and shell placement during photography, a set of replicate images with digitized landmarks was used to calculate the disparity using the morphol.disparity function in geomorph (Adams and Otarola-Castillo 2013, Adams et al. 2014).

Specimens were grouped by MOTUs to identify detectable differences in shell variation in concordance with genetic variation. Reduced datasets of those specimens with obtained genetic data were used for these groups. However, in the absence of genetic data, all specimens for which data are available were grouped by the physiographic province associated with the collection locality (i.e., Cumberland Plateau, Eastern Highland Rim, Valley and Ridge, and Western Highland Rim). These regions possess unique environmental characteristics (e.g., soil physiochemistry, rock type, vegetation; Fenneman 1917) and were utilized as broad-scale categories to test for the effect of environmental variation on conchology. Before assessing the significance of these groups in explaining morphological variation, the first and second PCs were subjected to a test of spatial autocorrelation (SAC) to prevent the increase of Type 1 errors introduced to the analyses (Perez et al. 2010). SAC was determined using the Moran’s I statistic (Sokal and Oden 1978) and was found to be non-significant (PC 1=0.2184, PC 2=0.0832). These groups were subjected to Procrustes Analysis of Variance (ANOVA) following a randomized residual permutation procedure (RRPP) for 10,000 iterations.

Seven unique shell measurements from Burch (1962) were utilized in the TM approach (see Figure 2B): Shell width (SW), shell height (SH), aperture width (AW), aperture height (AH), body whorl height (BW), penultimate whorl height (PW), and angle of apex (AA). These shell characteristics are often utilized in morphometric analyses and are readily utilized in land snail species identification guides (Pearce and Örstan 2006, Dourson 2010). Scaled data were converted to a Euclidean distance matrix and subject to Permutational Multivariate Analysis of Variance (PERMANOVA). This test was performed in the vegan package in R for 10,000 permutations (Oksanen 2018). P-values extracted from pairwise comparisons were corrected using a Bonferroni test.


Genetic analyses

Molecular sequence data were obtained from 32 specimens of 23 populations. Tissue samples were scarce, and the saturation of the land snail soft body with mucopolysaccharides inhibited the success of standard extraction procedures and subsequent sequencing. Thus, we were unable to obtain full genetic coverage (i.e., all four target genes sequenced) for all specimens. Summary statistics generated for the genetic data is presented in Table 2. The mtDNA dataset used for the gene tree estimation was unambiguously aligned (1316 base pairs; bp). A concatenated alignment of all specimens in which all four genes were amplified was also unambiguously aligned and assessed (n = 16; 3040 bp). The CO1 dataset was also assessed independently, as it was later utilized for the downstream species delimitation approaches. The CO1 dataset was unambiguously aligned (n=24; 704 bp). No shared haplotypes were observed between cave populations at CO1 (see Table 3), even in cases where caves were less than 15 m apart from one another (e.g., Demarcus Cave (AN5) and Robert Smith Cave (AN6) in Anderson County, Tennessee). The generated haplotype network strongly resembles MOTU delimitation results (see Figure 3), with the two most diverse MOTUs identified possessing five haplotypes each. Mean uncorrected p-distances between cave populations at CO1 was 16.14% (range 2.6–23.2%), indicating significant geographic isolation. For the concatenated genetic dataset, mean uncorrected p-distances between cave populations was 6% (range 1.3–10.3%). Due to the rarity of this species, there were only four instances of obtaining sequences of more than one individual per cave (Columbia Caverns (DI6), Keith Cave (FR14), Offut Cave (AN12), Panther Cave No. 1 (CM8)). Of these, two populations exhibited two haplotypes at CO1. Intrapopulation variation of these four populations was low, with a mean CO1 uncorrected p-distance of 1.48±0.3%.

Summary statistics generated for all four genes assessed (mtDNA: CO1, 16S; nDNA: 28S, H3).

n bp h Hd Π S Eta
CO1 24 668 18 0.957 ± 0.031 0.14346 ± 0.012 164 245
16S 28 605 20 0.968 ± 0.019 0.09019 ± 0.011 71 99
28S 29 1305 6 0.424 ± 0.111 0.00327 ± 0.001 20 20
H3 25 337 6 0.427 ± 0.122 0.00322 ± 0.001 7 7

Species delimitation results from both ABGD and mPTP analyses. Haplotype diversity, specimen ID, state, karst region, and physiographic province also included.

mPTP ABGD Haplotype ID Specimen ID State Karst Physiographic
PG1 AG1 H1, H2 NSG-CM8-T1, NSG-CM8-T2, MLN-15-006.19 TN APP VR
PG2 AG2 H3 NSG-AN12-T1, NSG-AN12-T2, NSG-AN12-T3, NSG-AN12-T4, NSG-AN12-T5 TN APP VR
PG3 AG3, AG4, AG5, AG6, AG7 H4, H5, H6, H7, H8 NSG-JK3-T1, NSG-VB9-LL, GC1, NSG-DI3-LL, NSG-FR14-T1, NSG-FR14-T2 TN ILP CP, EHR, WHR
PG5 AG9 H10 MLN-14-0153 TN ILP CP, EHR
PG6 AG10, AG11, AG12 H11, H12, H13 AUM27534-T2, NSG-KN50-NP, NSG-AN6-NP TN APP, ILP VR, WHR
PG7 AG13, AG14, AG15, AG16 H14, H15, H16, H17, H18 MLN-16.0228, NSG-KN112-T1, NSG-DI6-T1, NSG-DI6-T2, AUM28173 AL, TN APP, ILP CP, VR, WHR
Figure 3. 

Haplotype network generated using the NeighborNet network method with uncorrected p-distances with the CO1 dataset. Species delimitation results are depicted using major color groups for the mPTP results, and subcolor groups for the ABGD results.

Phylogenetic analyses

Both ML and BI approaches resulted in highly similar tree topologies for each unique concatenated genetic dataset (CO1, mtDNA, mtDNA + nDNA). The outstanding difference between the ML and BI phylograms generated from the mtDNA dataset was a resolution of polytomy from the ML approach in the Bowman Cave (DI3), Carter Cave (JK3), Keith Cave (FR14), McCoy Cave (VB9), and Slippery Slit Cave (OV440) clade. The mtDNA + nDNA phylograms also differed with the Hering Cave (AMD6) population representing a monotypic group in the ML approach, and grouping with the Brent’s Cave (KN112), Columbia Caverns (DI6), and Weavers Cave (AN22) clade as it does in all other phylograms assessed. Additionally, there were several notable distinctions in the CO1 phylograms produced between BI and ML approaches (see Suppl. material 2). Despite these differences, only the representative phylogenies utilizing the BI approach for the CO1, mtDNA, and mtDNA + nDNA datasets are shown (Figures 4, 5). All other trees are placed within Suppl. materials 2, 3.

Figure 4. 

Phylogenetic tree of the CO1 dataset (808 bp) using the BI methodology. Posterior probabilities generated from the analysis are shown for each clade with the top numbers. Confidence values given from the bootstrapped ML method are shown for each clade with the bottom numbers. The ‘x’ symbols indicate varying topology between the BI and ML analyses. ML trees are reported in the Appendix for cross-reference. Species delimitation results are depicted using major color groups for the mPTP results, and subcolor groups for the ABGD results.

Figure 5. 

Phylogenetic trees of the concatenated mtDNA (CO1 + 16S; 1316 bp) and the full mtDNA + nDNA (CO1 + 16S + 28S + H3; 3040 bp) datasets. Posterior probabilities generated from the analyses are shown for each clade with the top numbers. Confidence values given from the bootstrapped ML method are shown for each clade with the bottom numbers. The ‘x’ symbols indicate varying topology between the BI and ML analyses. ML trees are reported in the Appendix for cross-reference. Species delimitation results are depicted using major color groups for the mPTP results, and subcolor groups for the ABGD results.

Due to an inability to amplify all genes per specimen, some specimens are not represented in all phylogenies. However, among the representatives included in all three datasets, there is a consistent topology. The only differences between the mtDNA tree and the mtDNA + nDNA BI trees are 1.) the resolution of polytomy and varied topology in the Bowman Cave (DI3), Carter Cave (JK3), Keith Cave (FR14), McCoy Cave (VB9), Slippery Slit Cave (OV440) clade, and 2.) the relative placement of the Frazier Hollow Cave (DK11), Robert Smith Cave (AN6), and Conner’s Creek Cave (KN50) clade. All other clades remain consistent. Bootstrap support for the ML approach were notably lower at deeper nodes in each phylogram, and the same occurred with posterior probabilities generated from the BI approach. Node posterior probabilities and confidence values increased overall after the addition of the nDNA data. Comparison of both mtDNA and nDNA phylograms show the existence of at least seven monophyletic clades across the Appalachians and ILP karst regions (Figure 5). The monotypic Dry Cave (FR9) and Demarcus Cave (AN5) samples seem to be considerably divergent from other groups. While the former is known from the southern extent of the Eastern Highland Rim, the latter monotypic clade is in immediate proximity to Robert Smith Cave (less than 15 m) yet both are significantly delineated in the CO1 phylogram and the subsequent delimitation approaches.

Species delimitation

The ABGD method generated two partition strategies. At prior intraspecific divergence (P) values between 0.0010 and 0.0215, sixteen MOTUs were recognized in initial and recursive partitions. Both partition schemes remained stable at these values until reaching congruency at P = 0.0359, grouping all populations together into a single MOTU. The barcode gap was discovered at 0.14–0.16 K2P distance. The PTP results generated seven MOTUs for both single and multi-coalescent rate models (see Suppl. material 4). Both delimitation approaches show highly similar MOTU designations, with most identified groups being known from individual caves (Figure 4). There were four cases of both delimitations methods producing the same results (PG1, PG2, PG4, PG5). Three groups of five, four, and three MOTUs generated by ABGD were consolidated into three MOTUs generated by mPTP (PG3, PG6, PG7), respectively. The consolidated PG3 MOTU group is largely clustered within the Eastern Highland Rim (AG3, AG4, AG5, AG7), with only one disjunct representative being found in a fragmented karst formation on the eastern extent of the Western Highland Rim (AG6). The PG6 and PG7 MOTU groups exhibit an irregular geographic structure, with both possessing representatives from both karst regions (Appalachians and ILP). Further, the ABGD results suggest two MOTU groups (AG15, AG16) within a single cave population at Columbia Caverns, with AG16 comprising this cave on the eastern extent of the Western Highland Rim and another in the southern extent of the Cumberland Plateau in the state of Alabama.

Morphometric analyses

In total, 65 specimens were incorporated into both the GM and TM datasets from 28 cave populations. The disparity test used to indicate possible error introduced from shell and landmark placement (2.28%) was negligible. PC plots for each grouping are displayed in Figure 6A–D. For the GM PCA, the first three principal components account for 69.52% of the total variance. PC 1 (31.24%) was interpreted as the curvature of the shell, with higher PC scores exhibiting a higher angle of apex and larger shell height. The X-coordinates of LM 3, LM 2, and LM 6 all had the highest PC loadings associated with PC 1 (0.418, 0.312, 0.243 respectively). PC 2 (25.18%) was interpreted as the size of the secondary body whorl in relation to the aperture, with higher PC scores exhibiting significantly wider secondary body whorls with an annular apertural structure. The Y-coordinate of LM 8 and the X-coordinates of LM 9 and LM 1 had the highest PC loadings associated with PC 2 (0.281, 0.189, 0.159 respectively). PC 3 (13.10%) was interpreted as the size of the aperture, with higher PC scores exhibiting larger apertures and higher shell width. The X-coordinates of LM 6 and SLM 15 and the Y-coordinate of LM 4 had the highest PC loadings associated with PC 3 (0.354, 0.333, 0.301 respectively). A smaller morphometric dataset (n = 39) was assessed for those individuals for which molecular data was available. Only the MOTU groups from the mPTP analysis were considered, as these were the larger groups. The first three principal components for this smaller dataset account for 74.90% of the total variance (PC 1 = 31.03%; PC 2 = 28.76%; PC 3 = 15.11%).

Figure 6. 

PCA results from both geometric morphometric (left) and traditional morphometric (right) analyses. A, B Total morphometric dataset (n=65) grouped by physiographic province. C, D Morphometric dataset with complimentary molecular data (n=39) grouped by MOTUs from the mPTP analysis.

For the TM PCA, the first three principal components accounted for 79.36% of the total variance. PC 1 (66.03%) was interpreted as the overall size of the shell, with high PC scores exhibiting larger shell height and shell width (PC loadings = 0.4553, 0.4370 respectively). PC 2 (13.23%) was interpreted as the height of the shell, with higher PC scores exhibiting much larger penultimate whorls and shell height (PC loadings = 0.6336, 0.6011 respectively). PC 3 (9.85%) was interpreted as the curvature of the shell, with higher PC scores having higher angles of apex and smaller shell width (PC loadings = 0.6815). For the smaller mPTP dataset, the first three principal components accounted for 85.98% of the total variance. Procrustes ANOVA and PERMANOVA tested the influence of environmental variation (i.e., respective physiographic province) on external shell morphology, indicating significance for both morphometric approaches. MOTU groups did not significantly explain shell variation with the GM approach, but it was significant for the TM approach.


Many molecular studies of troglobitic taxa have revealed previously unknown cryptic lineages in North America (Buhay and Crandall 2009, Snowman et al. 2010, Niemiller et al. 2012, Weckstein et al. 2016). Troglobites are hypothesized to have fewer opportunities for dispersal than obligately-subterranean aquatic species (i.e., stygobites), due to limited connectivity of terrestrial subterranean passages (Culver et al. 2009). This may promote isolation and short-range endemism in troglobites (Culver et al. 2009, Niemiller and Zigler 2013). No phylogeographic study of troglobitic snails has been conducted in North America, and all other molecular studies of troglobites in the Appalachians and Interior Low Plateau have focused on organisms with comparatively higher vagility and dispersal potential (e.g., Buhay et al. 2007, Niemiller et al. 2008, Niemiller et al. 2012, Snowman et al. 2010, Loria et al. 2011). Using both a multilocus molecular and a morphometrics approach, we investigated genetic diversity within H. barri to identify potential cryptic populations within the species’ range, and to further determine whether external shell morphology was a useful indicator of differing patterns of genetic variation.

Genetic diversity of Helicodiscus barri

Despite limited sampling success of this rare species, our study revealed high genetic diversity in H. barri. Haplotypic diversity is strongly dictated by individual caves, and there appears to be little to no dispersal between cave systems regardless of proximity. Mitochondrial genetic divergence among H. barri populations is significantly higher (16.14%) compared to other troglobitic invertebrate taxa studied in the region (e.g. 3.1% for Nesticus spiders (Snowman et al. 2010); 0.06% for Tetracion millipedes (Loria et al. 2011); 2.1% for Ptomaphagus beetles (Leray et al. 2019)), suggesting that the low vagility of land snails accentuates the isolation caused by subsurface habitat fragmentation. Rates of mitochondrial gene evolution for land snails vary considerably, with estimates of 1.6–12.9% per million years for ribosomal genes and 2.8–13% for CO1 (Thomaz et al. 1996, Chiba 1999, Van Riel et al. 2005). Further, land snails often exhibit high levels of intraspecific genetic divergence and population structure (Guiller et al. 1994, Davison et al. 2009, Perez et al. 2014). An estimated 1.6% divergence per million years has been a proposed standard for other gastropods (Liu and Hershler 2007, Murphy et al. 2012, Harris et al. 2013). With this conservative estimate, CO1 sequence divergence suggests the average timing of isolation between H. barri populations is 10.1 million years, and up to 14.5 million years. In this scenario, not only do these results indicate the evolutionary independence of these cave populations, they suggest that the subterranean colonization of this species predates Pleistocene glaciation. The Climatic Relict Hypothesis suggests that environmental stress (such as the often implicated Holsinger (1988) “Pleistocene-effect” model) drives the colonization of organisms into subterranean environments (Leys et al. 2003, Culver and Pipan 2009). Using this gastropod CO1 molecular clock, the Climate-relict hypothesis is not supported. Rather, a scenario in which a geographically widespread proto-troglobitic (i.e., troglophilic) species colonized different subterranean systems independent of obvious environmental stress is favored instead.

Mitochondrial divergence estimates of ≥10% per million years are, however, often associated with terrestrial gastropods on island systems (Chiba 1996, Thacker and Hadfield 2000, Van Riel et al. 2005), which are highly comparable to cave systems due to the isolation of subterranean environments and the discontinuity of these habitats across karst landscapes (Culver 1970, Snowman et al. 2010). In this latter scenario with a high rate of mitochondrial gene evolution (10% per million years), average timing of isolation is 1.6 million years. This would suggest a climatically-driven subterranean colonization during the mid-Pleistocene, failing to reject the Climatic Relict Hypothesis. Thus, it is difficult to differentiate between varying biogeographic scenarios without the application of a Helicodiscus-specific molecular clock model. However, the development of an accurate molecular clock is problematic due to a notable gap in fossil material for the genus in North America. There are several occurrences of fossil material in surface and cave habitats across central and eastern North America during the Pleistocene (Rinker 1949, Wetmore 1962, Slaughter 1966, Schultz and Cheatum 1970, Guilday et al. 1977, Dalquist and Stangl 1984, Eshelman and Hager 1984), and one fossil record in central North America in the upper Miocene (Liggert 1997). Moreover, dating based on biogeographic barrier formation is also problematic, as timing estimates of cave formation across the distribution of H. barri are highly variable. The formation of some caves in the eastern Appalachians and Cumberland Plateau have been estimated to occur in the late Pliocene to middle Pleistocene (Davies 1953, Anthony and Granger 2004, White 2009), while other caves along the Tennessee River Valley and the Highland Rim have been estimated to form in the late Mesozoic to the early Tertiary (Moneymaker 1948, Barr 1961). Therefore, assessing the timing of colonization is currently beyond the scope of this study.

Delimitation analyses revealed up to sixteen unique MOTUs within H. barri, largely organized by geographic and geological similarity (see Figure 7). Most MOTUs belong to similar rock groups, arranged largely in association with each respective physiographic province. There were two unique cases of MOTUs being distributed across both the Appalachians and ILP karst, each exhibiting irregular geographic structure. PG7 is distributed across four cave populations from the northeastern Valley and Ridge, the southernmost contact zone of the Cumberland Plateau and the Eastern Highland Rim, and the westernmost extent of the Western Highland Rim. PG6 is distributed across three cave populations in the eastern Central Basin and the northeastern Valley and Ridge. Further, the ABGD results reveal two distinct MOTUs (AG15, AG16) within a single cave population at Columbia Caverns (DI6). AG15 is comprised of a single individual from Columbia Caverns, whereas AG16 is comprised of one individual from Columbia Caverns and another from the Hering Cave population in northern Alabama. This pattern may be the product of multiple cave colonization events in Columbia Caverns, or perhaps this demonstrates a case of sympatric speciation because of niche partitioning (e.g., Cooper et al. 2002, Niemiller et al. 2008), as these individuals were found in two separate areas of this large cave system. However, due to a low sample size, the aforementioned limited fossil data, and the uncertainty in estimating biogeographic barrier formation, it is difficult to determine the evolutionary history of this species and the geologic context whereby these unique MOTU groups may have developed.

Figure 7. 

Geographic distribution of MOTUs generated from the mPTP delimitation method in relation to karst adapted from Weary and Doctor (2014). Triangles represent cave populations. The numbers associated with each unique color corresponds to the associated mPTP MOTUs found in Table 3.

Utility of shell morphometrics in species delimitation of cryptic terrestrial micromolluscs

There has been much debate regarding the use of gastropod shell morphology in phylogenetic analyses (Emberton 1995, Wagner 2001, Uit de Weerd et al. 2004, Smith and Hendricks 2013, Miller 2016). Shell variation, while informative at lower taxonomic resolutions (e.g., Smith and Hendricks 2013), may not be useful in accurate delimitation of cryptic lineages, owing to the high responsiveness of shell structure to environmental factors and commonality of local adaptations in land snails (Goodfriend 1986, Fiorentino et al. 2008, Stankowski 2011, Razkin et al. 2017). Moreover, though many subterranean taxa (including Helicodiscus barri) exhibit disjunct, fragmented distributions, the ecological similarity of subterranean environments can lead to the protraction of morphological distinguishability between distinct genetic lineages (Losos and Mahler 2010, Eme et al. 2018, Inäbnit et al. 2019). Terrestrial micromolluscs pose additional difficulty in morphological delimitation due to their small size and similarities in external shell morphology, and molecular approaches have been favored (e.g., Weigand et al. 2012). Recent study of troglobitic Zospeum snails show that external shell morphology shows high variability both within and between cave populations, further obscuring the taxonomic identity of these cryptic groups without molecular data and intensive study of internal shell structure and soft tissue histology (Jochum et al. 2015).

Results herein indicate geographic variation of shell morphology as shown by the distinction of physiographic province groups, although intensive study of habitat variation was not performed. Both GM and TM methods resulted in significant differences among physiographic provinces. These findings further suggest an environmental influence on overall external shell morphology, agreeing with previous studies (Goodfriend 1986, Fiorentino et al. 2008, Vergara et al. 2017). The smaller MOTU dataset, comparatively, exhibited large morphological overlap between MOTU groups. However, results were significant for distinction from the TM groupings (see Table 4). This significance is most likely the result of population-based similarity in shell size, rather than the respective MOTU, of which many consist of multiple populations. Further, the GM groups were not significantly different for the MOTU dataset. The small sample size is a potential drawback to the utilization of these morphometric approaches. Terrestrial micromolluscs are notoriously difficult to sample (Boag 1982, Durkan et al. 2013), and sampling in cave environments significantly increases this difficulty. Moreover, as many of these populations remain understudied, morphometric methods may negatively impact populations subject to high amounts of collection and disturbance. This said, application of molecular barcodes may be most useful in the identification of these terrestrial micromolluscs (Weigand et al. 2011, 2014).

Results from both TM and GM analyses. Asterisk (*) denotes significant p-values.

Group Degrees of Freedom Sums of Squares R2 F P
Procrustes ANOVA
MOTUs (mPTP) 4 0.00763 0.14037 1.388 0.1311
Physiographic Province 3 0.01291 0.13503 3.1742 2.00E-04*
Permutational MANOVA
MOTUs (mPTP) 4 80.084 0.30107 3.6614 3.00E-04*
Physiographic Province 3 69.84 0.15589 3.7553 0.0021*

Taxonomic and conservation implications

The discovery of cryptic evolutionary lineages within H. barri has significant conservation implications. Recent reassessment of the conservation status of H. barri listed this species as Vulnerable (G3) under NatureServe criteria and Least Concern (LC) under the IUCN Red List criteria (Gladstone et al. 2018). Though our study suggests that this species is more geographically wide-spread than previously known, the distribution of individual MOTUs is greatly reduced, sometimes being restricted to a single cave. However, this species’ presence in both karst regions despite separation by a considerable amount of non-karst strata, and the discovery of a single specimen from surface habitat (see Table 1) suggests that it may not be limited to cave systems. Rather, like other Helicodiscus species, H. barri could be highly calciphilic, dwelling in rock talus piles or potentially interstitial habitats (Gladstone et al. 2018, Dr. Jeff Nekola, personal comm.). Few studies have investigated the significance of epikarst and other subsurface habitats to troglobitic fauna in North America (Culver et al. 2012), and a more intensive sampling effort may be necessary to assess the importance of these habitats in facilitating the dispersal of such snail fauna.

This study offers an important first step in outlining the presence of cryptic lineages within H. barri. However, many aspects of this species’ ecology and life history remain unknown, and the subsequent assessment of distinguishing ecology or habitat requirements for these cryptic groups is essential for their conservation and management. As with other recently discovered cryptic species, additional study of MOTU distribution, ecology, and conservation status are all necessary (Niemiller et al. 2013, Schlesinger et al. 2018).


We thank Annette Engel, Evin Carter, Charles Stephen, Denise Kendall Niemiller, Kirk Zigler, Lindsey Hayter, and Nathaniel Mann for assistance with collections, and Dan and Judy Dourson for assistance with specimen identification and helpful advice. This project was supported by the University of Tennessee, Knoxville, the National Speleological Society, and the Conchologists’ of America Organization.


  • Adams DC, Otarola-Castillo E (2013) geomorph: an R package for the collection and analysis of geometric morphometric shape data. Methods in Ecology and Evolution 4: 393–399.
  • Adams DC, Otarola-Castillo E, Sherratt E (2014) geomorph: Software for geometric morphometric analyses. R package version 2.0. packages/geomorph/index.html
  • Anthony DM, Granger DE (2004) A Late Tertiary Origin for Multilevel Caves Along the Western Escarpment of the Cumberland Plateau, Tennessee and Kentucky, Established by Cosmogenic super (26) Al and super (10) Be. Journal of Cave and Karst Studies 66(2): 46–55.
  • Armbruster JW, Niemiller ML, Hart PB (2016) Morphological Evolution of the Cave-, Spring-, and Swampfishes of the Amblyopsidae (Percopsiformes). Copeia 104(3): 763–777.
  • Barker GM (2001) Gastropods on land: phylogeny, diversity, and adaptive morphology. In: GM Barker (Ed.) The Biology of Terrestrial Molluscs. CABI Publishing, Wallingford, New Zealand, 1–146.
  • Barr TC (1961) Caves of Tennessee. Bulletin 64. Tennessee Division of Geology, Nashville, TN.
  • Bickford D, Lohman DJ, Sodhi NS, Ng PKL, Meier R, Winker K, Ingram KK, Das I (2007) Cryptic species as a window on diversity and conservation. Trends in Ecology & Evolution 22(3): 148–155.
  • Boag DA (1982) Overcoming sampling bias in studies of terrestrial gastropods. Canadian Journal of Zoology 60: 1289–1292.
  • Boc A, Diallo AB, Makarenkov V (2012) T-REX: a web server for inferring, validating and visualizing phylogenetic trees and networks. Nucleic Acids Research 40(1): 573–579.
  • Bookstein FL (1997) Morphometric tools for landmark data: geometry and biology. Cambridge University Press, Cambridge, UK.
  • Buhay JE, Moni G, Mann N, Crandall KA (2007) Molecular taxonomy in the dark: evolutionary history, phylogeography, and diversity of cave crayfish in the subgenus Aviticambarus, genus Cambarus. Molecular Phylogenetics and Evolution 42(2): 435–448.
  • Buhay JE, Crandall KA (2009) Taxonomic revision of cave crayfish in the genus Cambarus subgenus Aviticambarus (Decapoda: Cambaridae) with descriptions of two new species, C. speleocoopi and C. laconensis, endemic to Alabama, USA. Journal of Crustacean Biology 29: 121–134.
  • Burch JB (1962) How to Know the Eastern Land Snails. In: William C (Ed.) Brown Company Publishers, Dubuque, IA: 214 pp.
  • Burress PBH, Burress ED, Armbruster JW (2017) Body shape variation within the Southern Cavefish, Typhlichthys subterraneus (Percopsiformes: Amblyopsidae). Zoomorphology: 1–13.
  • Cameron RAD, Pokryszko BM (2005) Estimating the species richness and composition of land mollusc communities: Problems, consequences and practical advice. Journal of Conchology 38(5): 529–548.
  • Clements R, Ng PKL, Lu X, Ambu S, Schilthuizen M, Bradshaw CJA (2008) Using biogeographical patterns of endemic land snails to improve conservation planning for limestone karsts. Biological Conservation 141(11): 2751–2764.
  • Colgan DJ, Ponder WF, Eggler PE (2000) Gastropod evolutionary rates and phylogenetic relationships assessed using partial 28S rDNA and histone H3 sequences. Zoologica Scripta 29(1): 29–63.
  • Coney CC, Tarpley WA, Bohannan R (1982) Ecological studies of land snails in the Hiawassee River Basin of Tennessee, U.S.A. Malacological Review 15: 69–106.
  • Cooper SJB, Hinze S, Leys R, Watts CHS, Humphreys WF (2002) Islands under the desert: molecular systematics and evolutionary origins of stygobitic water beetles (Coleoptera: Dytiscidae) from central Western Australia. Invertebrate Systematics 16(4): 589–590.
  • Culver DC, Pipan T (2009) The biology of caves and other subterranean habitats. Second edition. Oxford University Press, Oxford, UK.
  • Dalquist WW, Stangl FB (1984) Late Pleistocene and early Recent mammals from Fowlkes Cave, southern Culbertson County, Texas. Carnegie Museum of Natural History Special Publication 8: 432–455.
  • Davies WE (1953) Geology of Pennsylvania caves. National Speleological Society Bulletin 15: 3–9.
  • Dinkins BJ, Dinkins GR (2018) An Inventory of the Land Snails and Slugs (Gastropoda: Caenogastropoda and Pulmonata) of Knox County, Tennessee. American Malacological Bulletin 36(1): 1–22.
  • Douglas DA, Dourson DC, Caldwell RS (2014) The Land Snails of White Oak Sinks, Great Smoky Mountain National Park, Tennessee. Southeastern Naturalist 13(1): 166–175.
  • Dourson DC (2010) Kentucky’s Land Snails and Their Ecological Communities. Goatslug Publications, Bakersville, NC.
  • Durkan TH, Yeung NW, Meyer WM, Hayes KA, Cowie RH (2013) Evaluating the efficacy of land snail survey techniques in Hawaii: implications for conservation throughout the Pacific. Biodiversity and Conservation 22(13–14): 3223–3232.
  • Emberton KC (1995) When shells do not tell: 145 million years of evolution in North America’s polygyrid land snails, with a revision and conservation priorities. Malacologia 37: 69–110.
  • Eme D, Zagmajster M, Delić T, Fišer C, Flot JF, Konecny‐Dupré L, Pálsson S, Stoch F, Zakšek V, Douady CJ, Malard F (2018) Do cryptic species matter in macroecology? Sequencing European groundwater crustaceans yields smaller ranges but does not challenge biodiversity determinants. Ecography 41(2): 424–436.
  • Eshelman R, Hager M (1984) Two Irvingtonian (medial Pleistocene) vertebrate faunas from north-central Kansas. Contributions in Quaternary Vertebrate Paleontology: a volume in memorial to John E. Guilday, Special Publication of Carnegie Museum of Natural History 8: 384–404.
  • Fenneman NM (1917) Physiographic subdivision of the United States. Proceedings of the National Academy of Sciences of the United States of America 3(1): 17–22.
  • Finston TL, Johnson MS, Humphreys WF, Eberhard SM, Halse SA (2007) Cryptic speciation in two widespread subterranean amphipod genera reflects historical drainage patterns in an ancient landscape. Molecular Evolution 16: 355–365.
  • Fiorentino V, Salomone N, Manganelli G, Giusti F (2008) Phylogeography and morphological variability in land snails: the Sicilian Marmorana (Pulmonata, Helicidae). Biological Journal of the Linnean Society 94: 809–823.
  • Folmer O, Black M, Hoeh W, Lutz R, Vrijenhoek R (1994) DNA primers for amplification of mitochondrial cytochrome c oxidase subunit I from diverse metazoan invertebrates. Molecular Marine Biology and Biotechnology 3: 294–299.
  • Gladstone NS, Carter ET, McKinney ML, Niemiller ML (2018) Status and Distribution of the Cave-Obligate Land Snails in the Appalachians and Interior Low Plateau of the Eastern United States. American Malacological Bulletin 36(1): 62–78.
  • Guilday JE, Parmalee PW, Hamilton HW (1977) The Clark’s Cave bone deposit and the Pleistocene paleoecology of the central Appalachian Mountains of Virginia. Bulletin of Carnegie Museum of Natural History 2: 1–87.
  • Guiller A, Madec L, Daguzan J (1994) Geographical patterns of genetic differentiation in the land snail Helix aspersa Müller (Gastropoda: Pulmonata). Journal of Molluscan Studies 60: 205–221.
  • Harris JD, Ferreira AF, De Frias Martins AM (2013) High levels of mitochondrial DNA diversity within oxychilid land snails (subgenus Drouetia Gude, 1911) from Sáo Miguel island, Azores. Journal of Molluscan Studies 79(2): 177–182.
  • Hasegawa M, Kishino H, Yano T (1985) Dating of the human-ape splitting by a molecular clock of mitochondrial DNA. Journal of Molecular Evolution 22(2): 160–174.
  • Hermsen EJ, Hendricks JR (2008) W(h)ither fossils? Studying morphological character evolution in the age of molecular sequences. Annals of the Missouri Botanical Garden 95(1): 72–100.
  • Holsinger JR (1988) Troglobites: the evolution of cave-dwelling organisms. American Scientist 76(2): 146–153.
  • Hotopp KP, Pearce T, Nekola JC, Slapcinsky J, Dourson DC, Winslow M, Kimber G, Watson B (2013) Land Snails and Slugs of the Mid-Atlantic and Northeastern United States. Carnegie Museum of Natural History, Pittsburgh, PA, USA.
  • Inäbnit T, Jochum A, Kampschulte M, Martels G, Ruthensteiner B, Slapnik R, Nesselhauf C, Neubert E (2019) An integrative taxonomic study reveals carychiid microsnails of the troglobitic genus Zospeum in the Eastern and Dinaric Alps (Gastropoda, Ellobioidea, Carychiinae). Organisms Diversity & Evolution: 1–43.
  • Jochum A, Slapnik R, Klussmann-Kolb A, Pall-Gergely B, Kampschulte M, Martels G, Vrabec M, Nesselhauf C, Weigand AM (2015) Groping through the black box of variability: An integrative taxonomic and nomenclatural re-evaluation of Zospeum isselianum Pollonera, 1887 and allied species using new imaging technology (Nano-CT, SEM), conchological, historical and molecular data (Ellobioidea, Carychiidae). Subterranean Biology 16: 123–165.
  • Juan C, Emerson BC (2010) Evolution underground: shedding light on the diversification of subterranean insects. Journal of Biology 9(3): 17.
  • Kano Y, Chiba S, Kase T (2002) Major adaptive radiation in neritopsine gastropods estimated from 28S rRNA sequences and fossil records. Proceedings of the Royal Society of London B: Biological Sciences 269(1508): 2457–2465.
  • Kapli P, Lutteropp S, Zhang J, Kobert K, Pavlidis P, Stamatakis A, Flouri T (2017) Multi-rate Poisson tree processes for single-locus species delimitation under maximum likelihood and Markov chain Monte Carlo. Bioinformatics 33(11): 1630–1638.
  • Kimura M (1980) A simple method for estimating evolutionary rates of base substitutions through comparative studies of nucleotide sequences. Journal of Molecular Evolution 16: 111–120.
  • Kumar S, Stecher G, Li M, Knyaz C, Tamura K (2018) MEGA X: Molecular Evolutionary Genetics Analysis across Computing Platforms. Molecular Biology and Evolution 35(6): 1547–1549.
  • Lanfear R, Calcott B, Ho SYW, Guindon S (2012) PartitionFinder: combined selection of partitioning schemes and substitution models for phylogenetic analyses. Molecular Biology and Evolution 29(6): 1695–1701.
  • Lefébure T, Douady CJ, Gouy M, Trontelj P, Briolay J, Gibert J (2006) Phylogeography of a subterranean amphipod reveals cryptic diversity and dynamic evolution in extreme environments. Molecular Ecology 15(7): 1797–1806.
  • Leray VL, Caravas J, Friedrich M, Zigler KS (2019) Mitochondrial sequence data indicate “Vicariance by Erosion” as a mechanism of species diversification in North American Ptomaphagus (Coleoptera, Leiodidae, Cholevinae) cave beetles. Subterranean Biology 29: 35–57.
  • Leys R, Watts CHS, Cooper SJB, Humphreys WF (2003) Evolution of subterranean diving beetles (Coleoptera: Dytiscidae Hydroporini, Bidessini) in the arid zone of Australia. Evolution 57(12): 2819–2834.
  • Liggert GA (1997) The beckerdie local biota (early Hemphillian) and the first Tertiary occurrence of a crocodilian from Kansas. Transactions of the Kansas Academy of Sciences 100(3–4): 101–108.
  • Loria SE, Zigler KS, Lewis JL (2011) Molecular phylogeography of the troglobiotic millipede Tetracion Hoffman, 1956 (Diplopoda, Callipodida, Abacionidae). International Journal of Myriapodology 5: 35–48.
  • Losos JB, Mahler LD (2010) Adaptive radiation: the interaction of ecological opportunity, adaptation, and speciation. In: Bell MA et al. (Eds) Evolution since Darwin: the first 150 years. Sinauer Associates, Sunderland, MA, 381–420.
  • McGuire JA, Linkem CW, Koo MS, Hutchison DW, Lappin KA, Orange DI, Lemos‐Espinal J, Riddle BR, Jaeger JR (2007) Mitochondrial introgression and incomplete lineage sorting through space and time: phylogenetics of crotaphytid lizards. Evolution: International Journal of Organic Evolution 61(12): 2879–2897.
  • Miller JP (2016) Geometric morphometric analysis of the shell of Cerion mumia (Pulmonata: Cerionidae) and related species. Folia Malacologica 24(4): 239–250.
  • Morrison JPE (1942) Preliminary report on mollusks found in the shell mounds of the Pickwick landing basin in the Tennessee River Valley. Bulletin of the Bureau of American Ethnology 129: 337–392.
  • Morse ES (1864) Observations on the terrestrial Pulmonifera of Maine, including a catalogue of all the species of terrestrial and fluviatile Mollusca known to inhabit the state. Journal of the (Portland) Maine Society of Natural History 1: 1–63.
  • 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(3): 846–866.
  • Niemiller ML, Fitzpatrick BM, Shah P, Schmitz L, Near TJ (2013) Evidence for repeated loss of selective constraint in rhodopsin of amblyopsid cavefishes (Teleostei: Amblyopsidae). Evolution 67: 732–748.
  • Palumbi S, Martin A, Romano S, McMillan WO, Stice L, Grabowski G (1991) The Simple Fool’s Guide to PCR. Version 2.0. University of Hawaii, Honolulu.
  • Pearce TA, Örstan A (2006) Terrestial Gastropoda. In: Sturm CF, Pearce TA, Valdés A (Eds) The Mollusks: A Guide to Their Study, Collection, and Preservation. American Malacological Society, USA, 261–285.
  • Perez KE, Defreitas N, Slapcinsky J, Minton RL, Anderson FE, Pearce TA (2014) Molecular phylogeny, evolution of shell shape, and DNA barcoding in Polygyridae (Gastropoda: Pulmonata), an endemic North American clade of land snails. American Malacological Bulletin 32(1): 1–31.
  • Perez SI, Diniz‐Filho JAF, Bernal V, Gonzalez PN (2010) Spatial regression techniques for inter‐population data: studying the relationships between morphological and environmental variation. Journal of Evolutionary Biology 23(2): 237–248.
  • Pilsbry HA (1948) Land Mollusca of North America (north of Mexico). Volume II, Part II. The Academy of Natural Sciences of Philadelphia. Philadelphia, PA.
  • Pons J, Barraclough TG, Gomez-Zurita J, Cardoso A, Duran DP, Hazell S, Kamoun S, Sumlin WD, Vogler AP (2006) Sequence-based species delimitation for the DNA taxonomy of undescribed insects. Systematic Biology 55(4): 595–609.
  • Oksanen J, Blanchet GF, Friendly M, Kindt R, Legendre P, McGlinn D, Minchin PR, O’Hara RB, Simpson GL, Solymos P, Stevens MHH, Szoecs E, Wagner H (2018) vegan: Community Ecology Package. R package version 2.5-1.
  • R Core Team (2014) R: A language and environment for statistical computing. R Foundation for Statistical Computing, Vienna, Austria. URL
  • Razkin O, Gomez-Moliner BJ, Vardinoyannis K, Martinez-Orti A, Madeira MJ (2017) Species delimitation for cryptic species complexes: case study of Pyramidula (Gastropoda, Pulmonata). Zoologica Scripta 46: 55–72.
  • Rinker GC (1949) Tremarctotherium from the Pleistocene of Meade County, Kansas. Contributions from the Museum of Paleontology, University of Michigan 7(6): 107–112.
  • Rohlf FJ (2015) The tps series of software. Hystrix, the Italian Journal of Mammalogy 26: 1–4.
  • Ronquist F, Teslenko M, Van Der Mark P, Ayres DL, Darling A, Höhna S, Larget B, Liu L, Suchard MA, Huelsenbeck JP (2012) MrBayes 3.2: efficient Bayesian phylogenetic inference and model choice across a large model space. Systematic Biology 61(3): 539–542.
  • Schlesinger MD, Feinberg JA, Nazdrowicz NH, Kleopfer JD, Beane JC, Bunnell JF, Burger J, Corey E, Gipe K, Jaycox JW, Kiviat E (2018) Follow-up ecological studies for cryptic species discoveries: Decrypting the leopard frogs of the eastern US. PloS One 13(11): e0205805.
  • Schultz GE, Cheatum EP (1970) Bison occidentalis and associated invertebrates from the late Wisconsin of Randall County, Texas. Journal of Paleontology 44(5): 836–850.
  • Slaughter BH (1966) The Moore Pit local fauna; Pleistocene of Texas. Journal of Paleontology 40(1): 78–91.
  • Smith UE, Hendricks JR (2013) Geometric morphometric character suites as phylogenetic data: extracting phylogenetic signal from gastropod shells. Systematic Biology 62(3): 366–385.
  • Snowman CV, Zigler KS, Hedin M (2010) Caves as islands: mitochondrial phylogeography of the cave-obligate spider species Nesticus barri (Araneae: Nesticidae). Journal of Arachnology 38: 49–56.
  • Stankowski S (2011) Extreme, continuous variation in an island snail: local diversification and association of shell form with the current environment. Biological Journal of the Linnean Society 104: 756–769.
  • Thacker RW, Hadfield MG (2000) Mitochondrial phylogeny of extant Hawaiian tree snails (Achatinellinae). Molecular Phylogenetics and Evolution 16(2): 263–270.
  • Thomaz D, Guiller A, Clarke B (1996) Extreme divergence of mitochondrial DNA within species of pulmonated land snails. Proceedings of the Royal Society of London Series B – Biological Sciences 263: 363–368.
  • Uit de Weerd D, Piel WH, Gittenberger E (2004) Widespread polyphyly among Alopiinae land snail genera: when phylogeny mirrors biogeography more closely than morphology. Molecular Phylogenetics and Evolution 33: 533–548.
  • Van Riel P, Jordaens K, Van Houtte N, Martins AMF, Verhagen R, Backeljau T (2005) Molecular systematics of the endemic Leptaxini (Gastropoda: Pulmonata) on the Azores islands. Molecular Phylogenetics and Evolution 37: 132–143.
  • Weckstein JD, Johnson KP, Murdoch JD, Krejca JK, Tayika DM, Veni G, Reddell JR, Taylor SJ (2016) Comparative phylogeography of two codistributed subgenera of cave crickets (Orthoptera: Rhaphidophoridae: Ceuthophilus spp.). Journal of Biogeography 43: 1450–1463.
  • Weigand AM, Jochum A, Pfenninger M, Steinke D, Klussmann-Kolb A (2011) A new approach to an old conundrum – DNA barcoding sheds new light on phenotypic plasticity and morphological stasis in microsnails (Gastropoda, Pulmonata, Carychiidae). Molecular Ecology Resources 11(2): 255–265.
  • Weigand AM, Götze MC, Jochum A (2012) Outdated but established?! Conchologically driven species delineations in microgastropods (Carychiidae, Carychium). Organisms Diversity & Evolution 12(4): 377–386.
  • Weigand AM, Jochum A, Klussmann-Kolb A (2014) DNA barcoding cleans house through the Carychiidae (Eupulmonata, Ellobioidea). American Malacological Bulletin 32(2): 236–245.
  • Wetmore A (1962) Notes on fossil and subfossil birds. Smithsonian Miscellaneous Collections 145(2): 1–17.
  • White WB (2009) The evolution of Appalachian fluviokarst: competition between stream erosion, cave development, surface denudation, and tectonic uplift. Journal of Cave and Karst Studies 71(3): 159–167.
  • Zharkikh A (1994) Estimation of evolutionary distances between nucleotide sequences. Journal of Molecular Evolution 39(3): 315–329.