Morphometrics and Phylogeography of the Cave-Obligate Land Snail <i>Helicodiscus barri</i> (Gastropoda, Stylommatophora, Helicodiscidae)

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

The numbers associated with each unique color corresponds to the associated mPTP MOTUs found in Table 3

Undocumented Biodiversity within Subterranean Ecosystems
Caves provide a model system in which to study the evolutionary processes and historical factors related to biogeography and speciation . These systems, characterized by geographic isolation and relatively simple biological communities, are often viewed as analogous to oceanic islands . Strong selective pressures and the isolation of subterranean ecosystems can result in morphological stasis among otherwise genetically distinct species, largely due to the parallel evolution of these lineages (Niemiller et al. 2012). Further, many troglobites (i.e., terrestrial cave-obligates) exhibit broad, mosaic distribution patterns which, in conjunction with this morphological stasis, often confound traditional means of delimitating species boundaries (Jochum et al. 2015). Thus, 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.
Despite the relatively inhospitable conditions of subterranean environments (when compared to the surface), a taxonomically diverse fauna has been documented within the caves of the United States. Over 1,200 species have been documented, with many still awaiting formalized species descriptions (Hobbs 2012. The Interior Low Plateau (ILP) and Appalachians karst regions are particularly species rich areas within the United States, exhibiting the first and second highest levels of North American troglobitic diversity, respectively (Culver et al. 2003). Karst regions are largely considered biodiversity hotspots (Myers et al. 2000), with high levels of endemism associated with rock outcrops and caves (Clements et al. 2006). New molecular tools (e.g., environmental DNA sampling methods) are now being utilized to continue documenting biodiversity within these subterranean habitats , and increased awareness of the fragility of these systems is promoting cave conservation and management .
Phylogeography is the study of historical processes that have impacted the modern geographic distribution of species' populations by utilizing genetic data. Phylogeographic studies can provide a greater understanding of the importance of the hydrological and geological barriers that contribute to species diversification and how they shape biogeographic patterns (Avise 2000). The relative roles of vicariance and dispersal within subterranean ecosystems has been of increasing interest in light of continued advances in molecular biology and genomic methods (Porter et al. 2007;). An increasing number of studies has investigated population genetic and phylogeographic hypotheses of subterranean fauna (e.g., Moulds et al. 2007;Snowman et al. 2010;Weckstein et al. 2016). Although conclusions have sometimes differed, these studies have greatly increased our understanding of colonization history, speciation, dispersal, and biogeography of troglobitic taxa . Regarding the conservation and management of such organisms, additional phylogeographic studies have uncovered considerable levels of cryptic species diversity (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 (Armbruster et al. 2016;Burress et al. 2017;Jochum et al. 2015). 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 ).

Influence of Geological and Climatic Processes on Subterranean Biodiversity
The Interior Low Plateau (ILP) and Appalachians karst regions in North America are considerably biodiverse and cave-rich areas. Tennessee alone has over 10,000 documented caves, representing nearly 20% of all caves within the United States (Niemiller and Zigler 2013).
Cave systems in eastern Tennessee are presumably more isolated due to the significant structural fragmentation of the Appalachians Ridge and Valley limestone by synclinal ridges of clastics, such as sandstone and shale. Conversely, limestones are exposed over larger areas within the Mississippian Plateaus in central Tennessee along the Highland Rim (Barr 1968). The discontinuity of Appalachians karst is largely supported by the high level of biological endemism, with many species being known from only a single cave (Christman et al. 2005;Niemiller and Zigler 2013). Troglobitic taxa within the ILP, on the other hand, have much greater dispersal potential, and have notably broader geographic distributions (Christman and Culver 2001;Niemiller and Zigler 2013).
Substantial climatic fluctuations during the Pleistocene have significantly influenced modern interpretations of biodiversity and the geographic distributions of species within North America (Webb III 1992). Many organisms were driven to extinction by unfavorable climatic conditions, and others began transitioning to new niches which subsequently favored increased diversification (Hewitt 1996). Subterranean ecosystems are often considered important climatic refugia due to the relative stability of these systems. Many troglobitic taxa groups such as carabid beetles (Barr 1969), fishes , crayfishes (Buhay 2007), and arachnids (Bryson et al. 2014) are thought to be Pleistocene relicts (i.e., relict lineages of surface species driven to extinction following the Pleistocene), colonizing climatically stable subterranean habitats during periods of glacial advancement or recession. This aligns with the "Pleistocene-effect" model proposed by Holsinger (1988).

Morphological Disparity within Terrestrial Micromolluscs
Land snails (Phylum Mollusca, Class Gastropoda) are a significantly species-rich animal group.
Over 24,000 species are currently recognized, and over 35,000 species are thought to exist globally (Barker 2001;Lydeard et al. 2004). The land snail fauna within eastern North America is exceptionally diverse, with over 500 species documented (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 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 under-sampled 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). Areas hypothesized to have comparatively higher levels of snail biodiversity have had varying and presumably insufficient collection efforts, with many cryptic species remaining unaccounted (Dourson 2007;Douglas et al. 2014;Dinkins and Dinkins 2018).
This has led to a lack in the documentation of intraspecific morphological variation of micromolluscs, while also obscuring accurate geographic ranges for these species (Nekola and Coles 2010). Thus, as most land snail species are delimited on the basis of conchology (i.e., shell morphology), many collections sustain a high incidence of misidentification of minute species (Hubricht 1985;Nekola and Coles 2010), while also ignoring patterns of genetic variation. The continued misidentification of species can have significant impacts on biodiversity assessments and subsequent conservation management (Bickford et al. 2007).
Strictly utilizing morphological data to delimit extant species in the post-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). 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). Moreover, utilizing morphometric analyses can inform the causal mechanisms for shape variation between gastropod populations (Vergara et al. 2017).
Terrestrial micromolluscs of the genus Helicodiscus are found throughout the eastern United States (Hubricht 1985). This genus is markedly known for its unique conchological sculpture, often exhibiting depigmented soft bodies and prominent spiraling striae on the shells of both surface and subterranean species (Figure 1). Many of these species are calciphiles, and two species -H. barri Hubricht, 1962 andH. notius specus Hubricht, 1962 -have even adopted a cave-obligate existence (Hubricht 1962). The distributions of these troglobites span both the 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 and H. punctatellus -have recently 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 within this genus among subterranean taxa. Further, 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).

Investigation Goals
In this project, I conduct the first study of the morphological variability and phylogeography of the cave-obligate land snail Helicodiscus barri. Recent cave bioinventory efforts within the ILP and Appalachians karst regions have allowed for the increased collection of the troglobitic species Helicodiscus barri. Other cave-dwelling Helicodiscus species have rarely been discovered (Hubricht 1962), and as such H. barri serves as the focal point of this investigation.
This species is known from caves within both the Appalachians and ILP karst regions, across multiple established 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. I examined museum accessions for H. barri while additionally sampling caves within the states of Tennessee, Alabama, and Georgia. Using phylogeographic approaches, I (1) assessed patterns of genetic variation (i.e., lasting gene flow, relative genetic similarity) of H. barri and employed two species delimitation methods (ABGD and mPTP) to infer the presence of cryptic lineages, and (2) tested if current species hypotheses based on conchology correspond with patterns of genetic variation. Further, I assessed morphological variation between H. barri populations using traditional morphometrics (TM) and landmarkbased geometric morphometrics (GM). MATERIALS AND METHODS

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 personhours. Snail specimens were preserved in 100% ethanol and identified using published keys, species descriptions, and examination by taxonomic specialists (Pilsbry 1948;Hubricht 1962;Dourson 2010;Daniel Dourson, personal ID). 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 all, 154 shells were examined ( Table 1). The geographic distribution of populations utilized within this study can be found in

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), the mt cytochrome oxidase subunit 1 (COI) locus using the primer pair LCOI490/ HCO2198 (Folmer et al. 1994), the nuclear 28S ribosomal RNA locus using the primer pair 28Sna1/28Sna2 (Kano et al. 2002), and the 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 (Huntsville, AL).

Genetic Analyses
Forward and reverse sequences were assembled into contigs and edited in Sequencher v.5.1 (Gene Codes Corporation). 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 were accessioned into GenBank (see Table A1). PartitionFinder v.2.1.1 (Lanfear et al. 2012) was used to determine the best model for the 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, Kishino, and Yano (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

Phylogenetic Analyses and Species Delimitation
Phylogenetic trees were inferred utilizing both a Maximum-Likelihood (ML) and Bayesianinference (BI) approach. The 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 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 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 eleven LMs were combined with two manually traced curves of three equidistant SLMs anchored on LM 4-5 and LM 5-6 (see Figure 4A). 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. 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 et al. 2013(Adams et al. , 2014RC 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 utilized 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 data set. 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 et al. 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 (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 nonsignificant (PC 1=0.2184, PC 2=0.0832). These groups will be subject 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

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 inhibits the success of standard extraction procedures and subsequent sequencing. Thus, I was unable to obtain full genetic coverage (i.e., all four target genes sequenced) for all specimens. Summary statistics generated for the genetic data can be found 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 data were 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 and Robert Smith Cave). The generated haplotype network strongly resembles MOTU delimitation results (Figure 5)   CO1. Intrapopulation variation of these four populations was low, with a mean CO1 uncorrected p-distance of 1.48%.

Phylogenetic Analyses
Both ML and BI approaches resulted in highly similar tree topologies for each unique there were several notable distinctions in the CO1 phylograms produced between BI and ML approaches (see Figure A1). Despite these differences, only the representative phylogenies utilizing the BI approach for the CO1, mtDNA, and mtDNA + nDNA datasets are shown in

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 (see Figure 8A). 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 (see Figure 8B). The PTP results generated seven MOTUs for both single and multi-coalescent rate models. Both delimitation approaches show highly similar MOTU designations, with most identified groups being known from individual caves. 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 grouped 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

Morphometric Analyses
In total, 65 individuals 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. 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 Xcoordinates 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 Xcoordinates 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%). PC plots for each grouping are displayed in Figure 9A-D.
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 testing the influence of environmental variation (i.e., respective physiographic province) on external shell morphology indicated 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 (Table 4).

CHAPTER 4 DISCUSSION
Many molecular studies of troglobitic taxa have discovered 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 thought to have fewer opportunities for dispersal than obligately-subterranean aquatic species (i.e., stygobites), due to limited connectivity of terrestrial subterranean passages. This may promote isolation and short-range endemism in troglobites 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. 2011;Snowman et al. 2010;Loria et al. 2011). Using both a multilocus molecular and morphometrics approaches, my goal for this study was to investigate genetic diversity within H.
barri to potentially identify cryptic populations within the species' range, and to further determine whether external shell morphology was a useful indicator of differing patterns of genetic variation. 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, intraspecific population structure is known to exhibit high genetic divergence (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 pre-dates the Pleistocene glaciation.

Genetic Diversity of Helicodiscus barri
The Climatic Relict Hypothesis, which 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;. This scenario rejects this hypothesis, and instead favor a scenario in which a geographically widespread prototroglobitic (i.e., troglophilic) species colonized different subterranean systems independent of obvious environmental stress.
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 of 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 infer an accurate scenario 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;Dalquest and Stangl 1984;Eshelman and Hager 1984), and one fossil record in central North America in the upper Miocene (Liggert 1997). Moreover, dating on the basis of 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 results reveal up to sixteen unique MOTUs within H. barri, largely organized by geographic and geological similarity (see Figure 10).  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 10. 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 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 (here 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). 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 in this study. 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 did not exhibit significance with 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, utilization of molecular barcodes may be most useful in the identification of these terrestrial micromolluscs (Weigand et al. 2011(Weigand et al. , 2014.  (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.

Taxonomic and Conservation Implications
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 Schlesinger et al. 2018).    Carter, looking at the unique animals that live underground.
Immediately following the completion of his BS, Nicholas was accepted into the Master of Geology program at the University of Tennessee Knoxville, beginning his graduate career in 2017. His research continues to focus on the evolution of cave animals and their conservation.