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.

Figure 1. 

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

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).

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 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 (http://wwwabi.snv.jusieu.fr/public/abgd/abgdweb.html) 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 (http://mptp.h-its.org) for the maximum 100,000 MCMC generations, with 25% of samples (25,000) conservatively discarded as burn-in.

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.5.16.0.0. 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 http://life.bio.sunysb.edu/morph/). 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.

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%.

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.

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.

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.

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.

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.

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).

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).