Phylogeography and population genetics of a headwater-stream adapted crayfish, Cambarus pristinus (Decapoda: Cambaridae), from the Cumberland Plateau in Tennessee

Assessments of genetic diversity for imperiled species can provide a baseline for determining the relative impacts of contemporary anthropogenic threats (e.g., habitat fragmentation) on population connectivity and identify historical factors contributing to population structure. We conducted a population genetics and phylogeographic assessment of the imperiled Pristine Crayfish (Cambarus pristinus) sampled throughout its range encompassing two morphologically distinct forms. Pristine Crayfish exhibit a disjunct distribution throughout lower order tributaries suggesting they are headwater-adapted species. The two morphologically distinct forms of the Pristine Crayfish are found in the upper Caney Fork (nominal Caney Fork form) and the Big Brush Creek and Cane Creek systems (Sequatchie form). We used 19 microsatellite loci and the cytochrome oxidase subunit I (COI) gene to assess population connectivity and genetic diversity of the Pristine Crayfish. Haplotypes recovered from the COI gene revealed that historic connectivity was maintained within each form of the Pristine Crayfish. However, the divergence between forms was higher (2.3%) than within forms (< 2.0%), suggesting each form is on an independent evolutionary trajectory, supporting recognition of the Sequatchie form as a distinct taxon. Microsatellite analyses for the Caney Fork form recovered a high degree of population isolation and support for six genetically isolated population. In addition, genetic diversity metrics per population and for the Caney Fork form were low suggesting that the Caney Fork form is at an increased risk of extinction under anthropogenic disturbances. We suggest that each form receive continued listing protection and conservation resources and that the Sequatchie form be treated as a unique taxon.


Introduction
Headwater streams play an important role in maintaining biodiversity of riverine systems. Although specific environmental attributes can vary across different landscape types, headwater streams exhibit environmental and ecological characteristics (e.g., gradient, in-stream habitat, water temperature, nutrient inputs) that differ from larger streams within the same river system (Richardson and Danehy 2007). Accordingly, these streams support a unique and diverse community of headwater-adapted species within the riverscape (Vannote et al. 1980). Additionally, these headwater streams are connected linearly or dendritically by a hierarchical arrangement of lotic habitats of increasing stream size. Hydrological and ecological barriers to species dispersal, such as shifts in stream size, gradient or physiography, have promoted isolation and speciation in headwater communities and contributed to patterns of population structure within species (Dudgeon et al. 2006, Strayer 2006, Strayer and Dudgeon 2010. For example, headwater-adapted species may experience reduced dispersal potential across larger rivers due to increases in predation or absence of cover habitat within these larger river habitats (Fraser et al. 1995;. The effect of natural instream features on population connectivity is often coupled with intrinsic biological characteristics of a species, such as dispersal ability and reproductive strategy (Turner et al. 1996;Turner and Trexler 1998;Hollingsworth and Near 2009).
Obligate riverine organisms are generally constrained to the linear or dendritic patterns of the streams they inhabit and often have patchy distributions in optimal habitats, due to the hierarchical, heterogenous nature of the riverscape (Vannote et al. 1980;Frissell et al. 1986;Fagan 2002;Tonkin et al. 2018). Nevertheless, many species can persist with sufficient recruitment between populations despite this limited dispersal (Brown and Kodric-Brown 1977). However, human-mediated habitat degradation and fragmentation will artificially increase population isolation and reduce recruitment. Increased isolation can reduce gene flow, which negatively affects population persistence and recolonization potential (Ward et al. 1994;Alp et al. 2012;Pilger et al. 2015Pilger et al. , 2017. Habitat fragmentation also leads to declines in population size, contributing to declines in genetic diversity through inbreeding and genetic drift. Reductions in genetic diversity may reduce the adaptive potential of impacted taxa and impair their ability to endure future stochastic events, possibly resulting in population or species extinction (Frankham 2005;Wright et al. 2008;Frankham et al. 2014).
Habitat degradation and fragmentation represent major threats to biodiversity and have been shown to drive species imperilment and extinction (Rands et al. 2010, Haddad et al. 2015. Headwater species are particularly susceptible to both and can be impacted in two distinct ways (Ward et al. 1994;Alp et al. 2012;Pilger et al. 2015Pilger et al. , 2017. First, these headwater habitats are closely associated with the surrounding terrestrial environment, so alterations within the drainage area can directly impact the in-stream environment ( Vannote et al. 1980;England and Rosemond 2004). Second, in-stream habitat alterations can create barriers to dispersal isolating previously connected populations (Sterling et al. 2012; Barnett et al. 2019;Blanton et al. 2019;Fluker et al. 2019).
Cambarus pristinus, the Pristine Crayfish, is a benthicadapted, headwater species endemic to the Cumberland Plateau of Tennessee (Withers and McCoy 2005;Rohrbach and Withers 2006;Williams et al. 2009;Johansen 2018). It has a patchy distribution within cobble-sized substrate of moderate pools within lower order stream reaches (< 4th) and tributaries of the Upper Caney Fork River (Cumberland River drainage) and Big Brush Creek (Sequatchie River drainage) watersheds ( Fig. 1a: Withers and McCoy 2005;Rohrbach and Withers 2006;Johansen et al. 2016;Johansen 2018). Current taxonomic work based on morphology suggests populations within Big Brush Creek (Sequatchie River drainage) and Cane Creek (Cumberland River drainage) watersheds are a separate species, referred to as the Sequatchie form, based on morphological differentiation from the northern, nominal form that resides in the Caney Fork and Bee Creek (Cumberland River system) watersheds (Fig. 1a). C. pristinus is generally rare at occupied sites (< 15 individuals/100 m 2 ) but is occasionally encountered at higher densities (42-59 individuals/100 m 2 ; Johansen 2018). Because small populations and habitat specialization are correlated to elevated extinction risk due to reduced genetic diversity and increased genetic drift (Purvis et al. 2000;Frankham 2005;Mace et al. 2008), C. pristinus is currently ranked as threatened by the Tennessee Wildlife Resources Agency (TWRA) and is under review for federal listing by the United State Fish and Wildlife Service (USFWS).
Major threats identified for C. pristinus are characteristic of rare species and include a limited range, habitat specialization, and small local densities (Taylor et al. 2007;Johansen et al. 2016). Anthropogenic disturbances such as surface mining and logging practices occur throughout the Cumberland Plateau and degrade aquatic habitat by increasing pollutants and changing water chemistry (i.e., increasing conductivity; McGrath et al. 2004;Schorr et al. 2006Schorr et al. , 2013Muncy et al. 2014;Gangloff et al. 2015) and are among leading threats to crayfish diversity (Taylor et al. 2007;Allert et al. 2013;Richman et al. 2015). Occupancy for the Caney Fork form has been positively associated with streams with reduced conductivity, moderate depths, and a high degree of cobble substrates (Johansen 2018) suggesting that disturbances that alter these parameters will have a detrimental effect on C. pristinus population persistence (McGrath et al. 2004;Schorr et al. 2006). Recent surveys have documented population loss for both forms after bridge construction and culvert alterations. In addition, the Sequatchie form has not been detected at several localities within the Big Brush Creek watershed where silviculture plantations dominate the terrestrial landscape of the watershed (Withers and McCoy 2005;Rohrbach and Withers 2006;Johansen et al 2016;Johansen 2018). These recent population losses and habitat fragmentation have likely increased the patchy distribution of C. pristinus and reduced their dispersal ability henceforth reducing gene flow and increasing the risk of further population loss and species extinction.
Several natural instream barriers exist within the range of C. pristinus that serve as potential dispersal barriers for benthic aquatic species. These barriers may have historically impacted population connectivity for C. pristinus resulting in long-standing population structure. The drainage divide between the Cumberland River and Sequatchie River systems is a potential barrier for the Sequatchie form of C. pristinus. Within the Caney Fork form, many occupied tributaries join the main stem Caney Fork River after it flows off the Western Escarpment of the Cumberland Plateau (Fig. 1a). The Western Escarpment is characterized by stark drops in elevation where many streams have, through erosion, confined, fast-flowing gulfs or ravines that may act as long-standing barriers or filters to gene flow among the tributary populations of the Caney Fork form (Bouchard 1975).
Finally, the preference of C. pristinus to occupy 2nd-4th order streams indicate that larger streams may serve as a barrier or filter to dispersal due to shifts in stream hydrology (Hollingsworth and Near 2009;Hurry et al. 2015;Schmidt and Schaefer 2018).
Because historical and contemporary factors may impact population connectivity, we examined multiple patterns of genetic structure within C. pristinus. Our first objective was to determine if there was support for long-standing isolation among populations and forms of C. pristinus due to perceived geological or geofluvial instream barriers. We assessed phylogeographic structure using the mitochondrial DNA (mtDNA) cytochrome oxidase subunit I (COI) gene. The mtDNA COI gene is commonly used to assess divergence among and within closely related crayfish species (Taylor and Knouft 2006;Glon et al. 2018Glon et al. , 2019Bloom et al. 2019). For our second objective, we used microsatellite loci to identify contemporary patterns of genetic diversity and population structure in C. pristinus. Microsatellites are often used for population genetics assessments (Pritchard et al. 2000;Selkoe and Toonen 2006;Abdul-Muneer 2014;Allendorf 2017;Micheletti and Storfer 2017;Blanton et al. 2019) and will provide a baseline for understanding genetic diversity in C. pristinus.

Sample collection
We surveyed 13 of 30 historical localities spanning the range of C. pristinus including populations representing both the nominal Caney Fork and Sequatchie forms (Table 1; Fig. 1a;Hobbs 1965;Withers and McCoy 2005;Rohrbach and Withers 2006;Johansen 2018). An additional four sites were sampled to increase collections of the Sequatchie form (Table 1). Most genetic samples were collected by removing one chela that was preserved in 95% ethanol. Whole specimens were Fig. 1 a Localities sampled for the two morphological forms of C. pristinus. Circles identify localities for the Caney Fork Form and diamonds denote those for the Sequatchie Form. Symbols filled with color denote localities where C. pristinus was encountered and colors corresponds to those used to depict phylogeographic relationships in (b, c). White-filled symbols denote localities sampled where C. pristinus was not detected. The white star is where the Caney Fork River becomes > 4th order. The Cumberland Plateau ecoregion is shaded gray. The thick black line represents the Western Escarpment of the Cumberland Plateau. Numbers denote site ID and correspond to those in Table 1: 1-Caney Fork, 2-Pokepatch Creek, 3-West Fork, 4-Meadow Creek, 5-Little Laurel Creek,  The two haplotype networks recovered from mitochondrial COI sequences of C. pristinus. Circles represent haplotypes and the size of each circle is proportional to the number of individuals with a given haplotype; connecting lines represent mutations among haplotypes. Colors correspond to localities depicted in Fig. 1a and open circles represent unsampled haplotypes. c Maximum-likelihood 50% majority rule consensus tree from the unique mitochondrial COI haplotypes recovered for C. pristinus. Bootstrap support values are found at each node and average pairwise sequence divergence estimates between clades or taxa are in parentheses collected if chelae were reduced/absent or if an individual's total carapace length (TCL) was < 16 mm to ensure enough tissue for genetic analyses.

DNA extraction
We extracted DNA from chela or abdominal tissue using a Qiagen DNEasy Blood and Tissue kit. We removed exoskeleton debris from all tissue prior to extraction because it can impair DNA quality and downstream PCR amplification of target loci. Two elutions were used to generate a 180 µL and 80 µL volume of extracted DNA. DNA was quantified on a NanoDrop ND-1000 Spectrophotometer (Thermo Fisher Scientific). Samples over 100 ng/µL were diluted to a concentration of 30 ng/µL for microsatellite locus amplification. Any concentrations under 100 ng/µL were used undiluted in microsatellite locus amplification. No dilutions were created for mitochondrial gene sequencing.

Mitochondrial DNA sequencing and analyses
Using previously published primers (Folmer et al. 1994), we amplified and sequenced the mtDNA COI gene for 51 individuals. Five individuals were examined per site, except for site 10 where only a single individual was collected (Table 1; Fig. 1a). Polymerase chain reactions (PCR) consisted of a 25 µL reaction volume with 2.00 µL of individual DNA, 18.75 mM MgCl 2 (New England Biolabs), 2.50 µL of Standard 10X buffer (New England Biolabs), 10 pM primer LCO1490 (5ʹ-GGT CAA CAA ATC ATA AAG ATA TTG G-3ʹ), 10 pM primer HCO2198 (5ʹ-TAA ACT TCA GGG TGA CCA AAA AAT CA-3ʹ), 1.25 U/mL Taq Polymerase (New England Biolabs), and 17.00 µL of PCR pure water. Thermocycler conditions included 35 cycles of the following: 30 s at 94 °C, 30 s at 50 °C, and 90 s at 72 °C, which occurred after an initial denaturation step of 3 min at 94 °C. These steps were followed by a final extension of 5 min at 72 °C. Sanger sequencing was performed by Yale's DNA Analysis Facility and resulting sequences were aligned and edited using CodonCode Aligner v8.0.1 (Soft Genetics).
To test for pseudogene presence in our sequence data, we followed recommendations of Buhay (2009), Song et al. (2008 and Sinclair et al. (2004) and performed visual inspections of sequence chromatograms for double peaks, examined sequence divergence between individuals for values 6% or higher, and performed BLAST searches to see if any of our sequences aligned with distantly relates species. We did not detect signatures of pseudogenes in our sequences and retained all sequences for further analyses.
We examined phylogeographic relationships among the resulting COI haplotypes by constructing a haplotype network in TCS v1.2.1 using a statistical parsimony analysis (Templeton et al. 1992) with a 95% connection limit. Final haplotype network graphics were created in TCSBeautifier  (Kumar et al. 2016). We also constructed a phylogenetic tree using MEGA v 7.0.26 (Kumar et al. 2016) using a maximum-likelihood analysis with 1000 bootstraps. The Hasegawa-Kishino-Yano with gamma distribution (HKY + G) substitution model (Hasegawa et al. 1985) was used as the parameter model of sequence evolution and was determined by the Akaike information criterion in MEGA's model selection test consisting of 24 possible substitution models (Nei and Kumar 2000). The final graphic was created in FigTree v 1.4.4 (Rambaut 2012).

Microsatellite genotyping
Hereditec (Langsing, NY, USA) developed species-specific microsatellite primers, which identified over 3000 potential loci after possible duplicate loci were removed. From this set we examined 82 loci for successful amplification and variation in the Caney Fork form of C. pristinus. The Sequatchie form was excluded from microsatellite analyses due to low sample sizes. Of the 82 loci examined, 42 were successfully amplified but only 19 loci were polymorphic and scoreable (Table S2) and retained for this study. PCR reactions consisted of 10.00 µL total volume including 1.00 µL of 10X concentrate standard Taq reaction buffer -Mg free (New England Biolabs), 30 mM MgCl 2 (New England Biolabs), 2 nM dNTPs (New England Biolabs), 2.5 pM forward primer, 5 pM reverse primer, 0.5 U Taq Polymerase, 0.1 µM M13-labeled dye (Applied Biosystems, Inc.), 1.00 µL DNA, and 5.65 µL PCR water. Thermocycler conditions consisted of 35 cycles of the following: 30 s at 94 °C, 30 s at the primer specific annealing temperature, and 90 s at 72 °C, which occurred after an initial denaturation step of 1 min at 94 °C. These steps were followed by a final extension of 5 min at 72 °C and held at 12 °C until the PCR product could be retrieved. Multiplexed PCR products were genotyped by the University of Florida Interdisciplinary Center for Biotechnology Research Genotyping and Gene Expression Core using the Liz 600 standard and an AB3730 sequencer. The resulting genotypes were scored in GeneMarker v1.6 (SoftGenetics), with manual confirmation and edits as needed.
Linkage disequilibrium (LD) for all locus pairs and departures from Hardy-Weinberg equilibrium (HWE) per locus were tested with GENEPOP v1.1.2 (Rousset 2008) in R. Parameters for HWE and LD were evaluated with the default Markov chain parameter settings of 10,000 dememorization steps, 1000 batches, 10,000 iterations per batch and p-values were adjusted following Bonferroni corrections to reduce Type I error (Rice 1989).

Sibling relationship inference
Samples of closely related individuals can bias allele frequency distributions (Allendorf and Phelps 1981). We used COLONY v2.0.6.7 (Jones and Wang 2010) to infer sibling relationships (sib-ships) between individuals within each sampled site. Sib-ships (SS) were inferred using a medium run-length with random mating and both male and female polygamy permitted. Sixteen individuals were removed from further analyses due to multiple SS inferences (Table 2).

Spatial genetic structure
Bayesian analyses implemented in STRU CTU RE v2.3.4 (Pritchard et al. 2000) were used to examine population structure for the Caney Fork form. Analyses with no a priori population assignment and using sampling sites as loc priors consisted of the following parameters: a generalized admixture model, correlated allele frequencies, 10 iterations for each value of K (K = 1-10), and 20,000 burn-in Markov chain Monte Carlo (MCMC) steps followed by 200,000 MCMC steps. The most likely number of population clusters was determined by using the mean log-likelihood (Ln[Pr(X|K)]) (Pritchard et al. 2000), ∆K (Evanno et al. 2005) and estimates from the Puechmaille (Puechmaille 2016) methods in the program StructureSelector (Li and Liu 1 3 2018). Because fine-scale population structure can be hidden within population clusters identified at larger spatial scales using the ∆K method (Vähä et al. 2007;Janes et al. 2017), we ran additional STRU CTU RE analyses independently for each population cluster identified in our analyses using sampling sites as loc priors and keeping all other parameters constant (Hubisz et al. 2009). The program CLUMPAK v1.1 pipeline (Kopelman et al. 2015) ran through StructureSelector (Li and Liu 2018) was used to summarize and present the output of independent runs of each K. Discriminant analysis of principle components (DAPC) can describe between-subpopulation variation while minimizing within-subpopulation variation noise. DAPC was conducted in Adegenet v2.1.1 in R (Jombart 2008). The optimal number of clusters (K) was determined using the k-means procedure with the function find.cluster and inferred from the Bayesian Information Criterion (BIC; Fig.  S1). Parameters were established using 10,000 iterations and 100 randomly chosen starting centroids to allow the algorithm to converge. Five principal components (PCs) were retained as determined by the function optim.a.score with ten simulations.
Pairwise F ST and Jost's D values were used to test for a signature of isolation-by-distance (IBD) using two pairwise geographic distance matrices measured in Google Earth. One matrix consisted of pairwise log-transformed riverine distances (km) and the second matrix consisted of pairwise log-transformed Euclidean distance (km) to account for potential overland movements (Oliveira et al. 2015;Thomas et al. 2018). In addition, significance was tested between geographic and genetic distances using a Mantel test with 10,000 randomizations of the data using IBD v1.52 (Bohonak 2002).
To highlight areas of increased/decreased gene flow within the range of C. pristinus, we ran the program EEMS (Estimated Effective Migration Surfaces; Petkova et al. 2016). EEMS defines deviations of genetic distance expected under a pure IBD model within a linear system, such as rivers. EEMS is also more robust than other analyses like STRU CTU RE or principal component analyses to uneven sample sizes (Petkova et al. 2016). Following methods described by Petkova et al. (2016), we adjusted parameters to achieve a 20-30% acceptance rate to improve MCMC convergence. We ran analyses using 500 and 700 demes with 10,000,000 MCMC steps and an initial burn-in of 1,000,000 steps for each analysis. Model convergence was assessed by looking at log-probability plots. Finally, all analyses were combined into a single figure using REEMSPLOTS in R (Petkova et al. 2016).

Genetic diversity estimates
Genetic diversity metrics for each population cluster recovered from our genetic structure analyses were calculated in the R package PopGenReport v3.0.4. Metrics calculated for each cluster included mean number of alleles per locus (N a ), percent of polymorphic loci, allelic richness (AR), private allelic richness (PAR), and observed (H o ) and expected (H e ) heterozygosity. The inbreeding coefficient (F IS ) was determined using GENETIX v4.05 (Belkhir et al. 2004) with 10,000 permutations to test for significance. All metrics calculated, except for PAR, were also estimated for the Caney Fork form overall.
Effective population size (N e ) was estimated using the LD method in NeESTIMATOR v2.01 (Waples and Do 2010;Do et al. 2014) to account for a single-year dataset and the SS method implemented in COLONY (Jones and Wang 2010;Wang et al. 2016). We included the SS method due to its robustness when dealing with sample size variability (Jones and Wang 2010). Effective population size was estimated per recovered population cluster and for the Caney Fork form with the upper and lower bounds determined by 95% CI and alleles with a frequency of < 0.02 excluded to prevent upward bias of N e (Waples 2006).
We used BOTTLENECK v1.2.02 (Piry et al. 1999) to test for recent population declines in each population cluster and for the Caney Fork form overall. We ran BOTTLENECK using a two-phase model (TPM) with 0, 10, and 20% multistep mutation rates (Zachariah Peery et al. 2012) with 36% variation (Di Rienzo et al. 1994) and tested for a significant excess of heterozygosity using a Wilcoxon sign-rank test.

Sample collections
We obtained tissue samples from ten of 13 sampled historical sites and one non-historical site. A total of 226 individuals were collected with at least 20 individuals per site except for sites 5, 7, 10, and 11, which had fewer than 20 individuals (Table 1). All sites were included in our mitochondrial analyses to assess phylogeographic relationships and only sites containing the Caney Fork form were used for our microsatellite analyses.

Mitochondrial DNA
We recovered 12 unique haplotypes from 51 individuals of C. pristinus sequenced for the mtDNA COI gene (Fig. 1b, Table S1). Eight haplotypes were recovered for the 40 Caney Fork form individuals and four haplotypes were recovered for the 11 Sequatchie form individuals sequenced. The 95% connection limit for haplotypes in our statistical parsimony analysis was 11 mutations. Average pairwise distance between sites ranged from 0.00 to 0.024768. We recovered two separate haplotype networks that exceeded the connection limit. One network consisted of all Caney Fork form haplotypes and the second included all Sequatchie form haplotypes (Fig. 1b). One common haplotype was shared among all Caney Fork form sites, except site 7 (Fig. 1b). A common haplotype was shared among two of the three sites for the Sequatchie form (Fig. 1b). Haplotypes within the Caney Fork form differed by one to three mutations and those within the Sequatchie by one to four mutations (Fig. 1b).
Similarly, two distinct clades, each representing a form, were recovered with 87% bootstrap support from our maximum likelihood analysis (Fig. 1c). There was low divergence and no definable geographic structure within forms (Fig. S2). Average pairwise distance using uncorrected p-distances between sites ranged from 0.00 to 2.5% (Table S3). Sequence divergence estimates from uncorrected p-distances were greater between forms at 2.3% than within forms (Fig. 1c). Sequence divergence within the Caney Fork form was lowest at 0.12%, and sequence divergence within the Sequatchie form was similarly low at 0.18% (Fig. S2).

Microsatellite marker validation
We successfully genotyped 165 individuals for 16 or more of the 19 loci examined (53/3135 missing genotypes; 1.8%). We found no evidence of scoring error, stutter, or allelic dropout for any locus. Evidence of null alleles was found in five loci but was inconsistent across sites. There was no evidence of LD or departures from HWE across loci. Therefore, we retained all 19 loci for subsequent analyses.

Sibling relationship inference
Sixteen SS pairs from five of the eight sites were identified in COLONY (Table 2). We removed one individual from each pair to avoid biasing population genetic results (Wang 2018). This reduced the number of individuals retained for further analyses to 149 (Table 3).

Spatial genetic structure
All STRU CTU RE analyses recovered six distinct population clusters under the mean log-likelihood and Puechmaille methods, except the mean log-likelihood analysis with no a priori populations designated, which recovered nine clusters. Our hierarchical analyses using the ∆K method were congruent with the six distinct population clusters recovered using the mean log-likelihood and Puechmaille methods. The six clusters identified included: sites 1, 2, and 5, recovered as a single cluster (Cluster A; Fig. 2c); and all other sites were each a distinct population cluster. Little to no admixture was observed among the six clusters (Clusters B-F; Fig. 2c).
Six distinct clusters, using five principal components (PCs), also were inferred from the DAPC analysis (Fig. 2b, Table 3 Genetic diversity measures based on 19 microsatellite loci for six distinct population clusters examined for the Caney Fork form of C. pristinus Cluster ID refers to each distinct population cluster recovered from our population structure analyses (Fig. 2) S1a-c). These six clusters were similar to the six distinct population clusters recovered from our STRU CTU RE analyses except that some individuals from Cluster A were assigned to other clusters in our DAPC; this discrepancy is likely due to the limits of the membership assignment algorithm in DAPC (Fig. S1d, e). Pairwise F ST (range: 0.002-0.450) and Jost's D values (range: 0.002-0.416) between sites were similar across several comparisons and the majority were significant after Bonferroni correction. Site 7 (Fig. 1a) consistently had the highest values under both metrics indicating that site 7 is the most genetically differentiated site (Table 4). The few sites where pairwise comparisons were not significant were between sites found within Cluster A recovered from our STRU CTU RE and DAPC analyses providing support for connectivity among sites within Cluster A and isolation of all other sites.
Our EEMS analysis highlighted several areas of increased and decreased gene flow (Fig. 2a) that are expected under IBD. Sites 6, 7, and 8 showed moderately decreased levels of gene flow while sites 3 and 4 showed slightly reduced to neutral levels of gene flow. Sites 1, 2, and 5 were recovered with low to moderately increased gene flow. Sites 1, 2, and 5 make up Cluster A from our STRU CTU RE and DAPC and suggests that a migration corridor that readily allows dispersal exists in the northern extent of the Caney Fork River watershed (Fig. 2a) despite the shift in stream hydrology as stream order increases. Sites 3, 4, 6, 7, and 8 correspond to population clusters C, B, D, F, and E respectively and highlight that the shift in stream topography along the Western Escarpment line reduces dispersal.

Genetic diversity estimates
A total of 113 alleles were amplified across all loci with an average of 5.78 alleles per locus across all clusters (Range: 2-12 alleles). Allelic diversity was low overall but relatively uniform among clusters, with N a ranging from 1.95 to 4.74  Fig. 1 and Table 1. b DAPC results showing first two discriminate functions. Colors represent the six distinct population clusters recovered from STRU CTU RE analysis and shapes indicate the original locality assignment for each individual. c STRU CTU RE analysis results using loc priors to assess population structure. Black boxes represent sites and colored bars within each box represent population of ancestry for each individual. Numbers denote site ID and letters denote population cluster assignment alleles and AR ranging from 1.83 to 3.79 alleles (Table 3). Observed (H o ) and expected heterozygosity (H e ) were similar across all clusters (range: 0.34-0.51 H o ; 0.38-0.54 H e ) except for cluster F (Ho = 0.18; He = 0.22; Table 3), which had relatively lower heterozygosity. Private allelic richness (PAR) was < 0.10 except for clusters E and F. Three clusters (A, B, and C) did not have F IS values that differed significantly from zero. The remaining clusters (D, E, and F) had significant positive F IS values indicating the presence of inbreeding for those clusters (Table 3). However, the 95% confident intervals spanned zero reducing our confidence in the significant F IS values recovered (Colegrave and Ruxton 2003). A signature of a bottleneck event was detected for cluster D with a 20% two phase model (TPM). There was no evidence for deviations from HWE after Bonferroni correction within any cluster (Table 3). N e using the LD method varied among clusters (39.2-∞) but was low overall for the Caney Fork form with an estimate of 26.3 individuals (95% CI 23.9-29.0). N e using the SS inference method was low for all clusters (23-49) and was low for the species overall with an estimate of 92 individuals (95% CI 68-121; Table 5).

Discussion
Hierarchical shifts in stream hydrology, changes in geology, and drainage boundaries can exist in lotic system as natural barriers to dispersal (Hughes et al. 1995;Fetzner Jr. and Crandall 2003). In addition, benthic, headwateradapted species naturally exhibit hierarchical population structure due to inherent biological factors such as low dispersal potential and habitat specificity (Fluker et al. 2014;Hurry et al. 2015;Schmidt and Schaefer 2018). Headwater species may stave off local extinction processes through dispersal, but human-mediated habitat fragmentation can disrupt dispersal and further isolate populations, reducing their long-term viability (Ward et al. 1994;Alp et al. 2012;Paz-Vinas et al. 2015). Table 4 Pairwise measures of genetic and geographic distances among sampled sites and used for the isolation-by-distance (IBD) analysis of the Caney Fork form of C. pristinus based on 19 microsatellite loci (Fig. 3) Site ID corresponds to those used in Fig. 1a  For conservation management, it is important to disentangle naturally occurring population structure from structure resulting from recent disturbances (Zellmer and Knowles 2009;Davis et al. 2014;Epps and Keyghobadi 2015). We expected to find signatures of long-standing isolation within C. pristinus attributed to reduced dispersal at three spatial scales: between the major drainage divides separating the Sequatchie River and Caney Fork River, among tributary systems of the Caney Fork River due to shifts in geology and topography associated with the physiographic break of the Western Escarpment of the Cumberland Plateau and among headwater tributaries due to shifts in hydrology associated with larger order streams. Although, we identified patterns of historic isolation between the Caney Fork and Sequatchie forms of C. pristinus, the observed structure was not related to any of our hypothesized barriers. Furthermore, we found no evidence for long-standing geographic structure within either of the forms of C. pristinus as indicated by short branch lengths and a common haplotype shared among the majority of sites within each morphological form. This suggests that within each form (or Clade) some level of gene flow was maintained historically among populations.
In contrast, the results from our microsatellite analyses identified contemporary patterns of population structure within the Caney Fork form. Based on our EEMs analysis, some of the observed structure and areas of reduced gene flow were related to the physiographic shift from the Cumberland Plateau to the Western Escarpment and, to a lesser degree, anincrease in stream size. This pattern was not reflected in our mitochondrial dataset, so we are unable to determine if the separation of the three upper Caney Fork clusters by the Western Escarpment from downstream clusters is due to long-standing isolation or related to contemporary patterns of human-mediated habitat alteration, or a combination of both. However, population structure observed within tributary systems not separated by these physiographic breaks, such as within clusters in the upper Caney Fork River, and the overall low levels of genetic diversity recovered across population clusters indicate recent declines and increased population isolation, likely associated with human habitat alterations that have negatively impacted the species. Low allelic richness, small N e , evidence of population bottlenecks, and signatures of inbreeding observed across population clusters and for the Caney Fork form as a whole, are indicative of a declining species at an elevated risk of extinction (Spielman et al. 2004).

Phylogeography of C. pristinus
Previous studies have shown that natural barriers to dispersal exist in lotic systems (e.g., hierarchical shifts in stream hydrology, changes in geology, and drainage boundaries) and that these barriers contribute to long-standing isolation in headwater-adapted aquatic organisms (Hughes et al. 1995;Fetzner Jr. and Crandall 2003;Nguyen et al. 2004; Table 4. Comparisons using log-transformed river distances are in black and those using log-transformed Euclidean distances are in gray (River-F ST : R 2 = 009; p = 0.13, Jost's D: R 2 = 0.14; p = 0.08; Euclidean-F ST : R 2 = 0.16; p = 0.05, Jost's D: R 2 = 0.16; p = 0.05). Jost's D values are indicated with circles and F ST values are indicated with diamonds Table 5 Effective population size (N e ) estimates for the LD and SS methods based on 19 microsatellite loci and population census size (N c ) estimates expressed as number of individuals/100 m from Johansen et al. (2016) for the Caney Fork form of C. pristinus A "-" indicates that the metric wasn't reported. Cluster IDs represent both site ID used in Table 1, denoted as a number and cluster ID's representing distinct population clusters recovered from our population structure analyses denoted in bold Hollingsworth Jr. and Near 2009;Kanno et al. 2011;Lamphere and Blum 2011;Hurry et al. 2015). However, the degree of population structure in headwater-adapted species varies with some like Etheostoma basilare (Corrugated Darter) exhibiting strong patterns of population structure attributed to historic vicariance (Hollingsworth Jr. and Near 2009;Clay et al. 2020), while others like Euastacus bispinosus (Glenelg Spiny Crayfish) and Etheostoma sagitta spilotum (Kentucky Arrow Darter) maintain connectivity across river systems and physiographic shifts within their range (Miller et al. 2014;Blanton et al. 2019). The latter suggests that some headwater-adapted species may be less constrained by inherent dispersal potential or factors such as tributary density, stream size, drainage shape, or distance (Turner and Trexler 1998;Schmidt and Schaefer 2018).
Based on the habitat specificity exhibited by C. pristinus, we predicted that larger river habitats and geographic breaks would serve as filters or barriers to dispersal, contributing to reduced gene flow among populations separated by these features. We recovered two genetically distinct and divergent clades representing the two morphological forms of C. pristinus, suggesting long-standing isolation of these forms. However, the geographic distribution of these two forms does not correspond to existing geographic breaks examined, given that the Sequatchie form spans the Cumberland-Tennessee River drainage divide. It is unlikely that the current distribution is explained by long distant dispersal due to the extensive distance (including most of the Tennessee and Cumberland River mainstems) separating extant populations and the lack of records of either form in streams of > 4th order in size that they would need to traverse. More likely, an ancestral population present on the Cumberland Plateau experienced a vicariant event resulting in the isolation of the two forms either in their present locations or between the Sequatchie and Caney Fork rivers; in the latter case, the isolation event may have been followed by range expansion through underground, overland, or headwater transfer of the Sequatchie form into its present range in the upper Caney Fork River system (Finlay et al. 2006;Berendzen et al. 2008;Niemiller and Zigler 2013). Estimates of divergence times from closely related taxa endemic to the Cumberland Plateau suggests the two forms of C. pristinus diverged sometime during or after the Pleistocene (Crandall et al. 2015). Climate oscillations and hydrological changes occurring during that time have been linked to the speciation of several lineages of freshwater taxa within the region (Thornbury 1965;Near et al. 2001;Berendzen et al. 2003Berendzen et al. , 2008Kozak et al. 2006). A study estimating divergence times for the two lineages would improve our ability to link the observed patterns to specific geological events and to explore alternative scenarios that may have led to divergence between the two forms.
Overall haplotype diversity within each morphological form was low, and two sites, 7 (Caney Fork form) and 10 (Sequatchie form), exhibited a single haplotype unique to each site. This may indicate a recent reduction in population size for each form and limited gene flow to sites 7 and 10 (Schrimpf et al. 2011). Alternatively, it is also possible that the low levels of haplotype diversity reflect the species' ecological or evolutionary history. Similarly, low haplotype diversity is observed in the broadly distributed C. parvoculus and C. jezerinaci (Mountain Midget crayfish and Spiny Scale crayfish, respectively) of the Cumberland Plateau (Thoma and Fetzner Jr. 2008). Low haplotype diversity has been observed in species with small N e and among species that have undergone demographic stochasticity or a historic range reduction (Bucklin and Wiebe 1998;Matocq and Villablanca 2001), which are observed features of C. pristinus.
In contrast to our mtDNA results that implied population connectivity across the range of the Caney Fork form, our contemporary levels of population structure inferred from our microsatellite data indicate that the Caney Fork form is comprised of six genetically differentiated populations. Patterns of isolation by distance (IBD) due to the hierarchical, linear nature of streams are expected for freshwater organisms (Wright 1943;Frissell et al. 1986 Fetzner Jr. andCrandall 2003;Sexton et al. 2013). However, we did not recover a significant relationship between genetic and river kilometer distances for the Caney Fork form (Fig. 3). We did recover a significant relationship between genetic distance and Euclidean distance but it only explained a small portion of the genetic differentiation observed. The lack of an IBD pattern using river distance and improbable occurrence of overland dispersal, despite a significant IBD association using Euclidean distance, suggests that something other than geographic distance contributed to population structure for the Caney Fork form. In headwater-adapted species with limited dispersal potential, the expected linear relationship between genetic and geographic distance may only occur at smaller geographic scales and become nonlinear at greater geographic scales. In such cases, habitat patch arrangement may be a more important predictor of genetic differentiation than geographic distance alone (van Strien et al. 2015). Under this scenario, a non-linear IBD relationship may occur at larger distances or spatial scales (Schmidt and Scahefer 2018). For example, IBD may be important in explaining population structure at finer spatial scales of the Caney Fork form, such as for population structure observed between Clusters A, B, and C, but landscape variables (i.e., Western Escarpment, shifts in hydrology of the mainstem Caney Fork River) may be better explanatory factors for population structure among all the Caney Fork form clusters.
As noted, our microsatellite data indicate a potential relationship between the Western Escarpment and reduced levels 1 3 of gene flow among populations separated by this feature. We also found fine-scale population structure among populations that crossed the Caney Fork River, suggesting reduced gene flow occurs as a stream size becomes > 4th order. These results suggest the shift from the Cumberland Plateau to the Western Escarpment may act as a natural barrier or filter to gene flow, and to a lesser extent, larger order rivers may also limit gene flow among populations within the Caney Fork form. However, the uplift of the Cumberland Plateau and subsequent incision of the Caney Fork tributaries into the Plateau likely occurred during the Pliocene and into Pleistocene (Sasowsky et al. 1995;Anthony and Granger 2007). If these features limit dispersal among populations, it is interesting that we did not see evidence of long-standing isolation around these features (Western Escarpment and > 4th order mainstem habitats) in our mtDNA results.
Although we cannot rule out that some degree of the observed contemporary population structure reflects reduced dispersal due to long-standing, natural instream features (i.e., the Western Escarpment of the Cumberland Plateau, large mainstem habitats), it is likely that recent metapopulation dynamics of C. pristinus have been disrupted due to human-mediated habitat degradation that has contributed to observed declines in occurrence, population size, and abundance (Johansen et al. 2016). Small populations, such as the small census estimates for C. pristinus (Johansen et al. 2016), are more sensitive to genetic drift, which in the absence of gene flow, can reduce genetic diversity leading to the loss of beneficial alleles and the fixation of detrimental alleles, and ultimately to the loss of adaptive potential (Frankham 1995a, b). Point estimates of N e for each population and the Caney Fork form overall were below 100. The 100/1000 rule states that an N e of 100 is required to prevent inbreeding depression and genetic drift while an N e of 1000 is necessary to maintain evolutionary potential. This rule provides a useful benchmark to assess the role population size has in the current conservation status of a species (Frankham et al. 2014). Several populations of the Caney Fork form, particularly those in the southern portion of the range, show signs of susceptibility to the negative effects of genetic drift and inbreeding via low N e estimates and high degrees of isolation.
Inbreeding is another potential threat to the long-term viability of isolated populations. Inbreeding exacerbates the effects of genetic drift and can interact with a small population size to push a species into a vortex that can rapidly push a species to extinction (Keller and Waller 2002;Frankham 2005;Wright et al. 2008). We did detect genetic signatures of inbreeding in several of our population clusters and for the Caney Fork form as a whole. This is concerning due to the small population size and genetic isolation we observed in the Caney Fork form.
Additionally, we detected a bottleneck in Puncheoncamp Creek (cluster D), suggesting a recent drastic decline in population size, which could reduce genetic diversity through the random loss of alleles. Although we do not have direct evidence of what caused this decline, historic mining practices within the headwaters of Puncheoncamp Creek (cluster D) and subsequent conversion to silviculture practices are a likely cause of population decline (Moore 1985;Withers and McCoy 2005;Schorr et al. 2006). Mining and silviculture activities negatively impact crayfish populations by increasing sedimentation and metal pollutants and altering water chemistry and riparian buffers (Allert et al. 2012(Allert et al. , 2013Helms et al. 2013;Loughman et al. 2015;Richman et al. 2015). The site for cluster D was different from other sites we visited for this study. The site consisted of a bedrock stream bed with few patches of large cobble substrate, a high degree of bank erosion, and a large amount of sediment deposition compared to our other sites. Anthropogenic disturbances are prevalent throughout the range of C. pristinus (Moore 1985;Withers and McCoy 2005;Schorr et al. 2006), and it is possible other clusters may have experienced bottleneck events, but our method requires a reduction of 50-80% of the effective population to detect bottleneck twothirds of the time (Hoban et al. 2013). It is also possible that other factors such as pre-bottleneck genetic diversity, bottleneck persistence, the timing of the bottleneck event, and population growth obscured genetic signals needed to test for bottleneck effects (Williamson-Natesan 2005; Zachariah Peery et al. 2012).
Signatures of isolation, inbreeding, and population bottlenecks indicate a species at an elevated risk of extinction and our assessment of genetic diversity for the Caney Fork form supports this conclusion. The loss of genetic diversity often precedes the distributional and demographic indicators of a declining species (Brook et al. 2002;Spielman et al. 2004). While not directly comparable due to variation in alleles, values for these metrics for the Caney Fork form were similar to those for other imperiled crayfishes and relatively lower than values associated with non-imperiled crayfish (Gouin et al 2006(Gouin et al , 2011Miller et al. 2014;Duncan et al. 2020; Table S4). The pattern of population structure observed in the Caney Fork form indicates gene flow is limited and is likely contributing to our low genetic diversity estimates. In concert with our mitochondrial assessment, C. pristinus appears to exhibit reduced adaptive potential when compared to non-imperiled crayfishes. However, more work is needed to explicitly test the role of contemporary land uses ongenetic diversity, adaptive potential, and population connectivity (Frankham 1995a, b;Matocq and Villablanca 2001).

Conservation implications
The morphological differences and genetic divergence observed between the two forms of C. pristinus, indicate that the forms are on independent evolutionary trajectories supporting the recognition of each form as a distinct species. We recommend that the Sequatchie form be provisionally referred to as Cambarus aff. pristinus. Due to the type locality residing within the upper Caney Fork River watershed (Hobbs 1965) and the holotype specimen representing a population of the Caney Fork form we recommend that only the Caney Fork form retains the Cambarus pristinus nomenclature. Management strategies should be developed independently for each lineage. This will be particularly important if translocation or captive propagation are utilized to supplement declining populations or re-establish gene flow between extant population (Vrijenhoek 1998).
Although our study indicates that IBD and natural instream features, primarily the shift to the Western Escarpment of the Cumberland Plateau and to a lesser degree stream size, contribute to the observed isolation observed in C. pristinus (Caney Fork form), they do not fully explain our genetic diversity results. Low genetic diversity measures combined with population loss, inbreeding, and bottleneck events suggest that contemporary anthropogenic impacts also negatively impact C. pristinus population connectivity. Additional work from a landscape genetics approach would provide a more robust test of effects of specific landscape and environmental factors that have contributed to the observed patterns of genetic diversity and isolation we recovered.
We identified patterns of low heterozygosity, low allelic diversity and inbreeding indicative of a species suffering the negative impacts of genetic stochasticity related to small population size and isolation (Brook et al. 2002). Although, C. pristinus may benefit from translocation of individuals among isolated populations, such actions would require a better understanding of the basic biology of this species (Ingvarsson 2001). For example, we do not know how contemporary population connectivity and effective population size estimates compare to historic patterns. Anecdotal evidence supports the idea that this species was never particularly abundant across its range (Hobbs 1965) and there is evidence that rare organisms possess traits that enable them to remain viable in a state of rarity (Kunin and Gaston 1993). Therefore, understanding the basic ecology (i.e., microhabitat requirements, territoriality, intraspecific competition, or sex-biased dispersal) and life-history traits will be vital to the success of conservation strategies utilizing translocation or captive propagation (Johansen 2018).
Regular genetic monitoring will be essential for successful conservation. Our estimates of N e were low and illustrate that only 10% of individuals estimated for census size contribute genes to subsequent generations. Thus, markrecapture census estimates may typically overestimate the health of the population. Additionally, populations will often exhibit genetic indicators of decline before population loss or size reduction (Brook et al. 2002;Frankenham et al. 2014); regular genetic monitoring will establish trends for the species and identify early warning signs for isolated populations and the species overall.
Although we could not assess genetic diversity and gene flow for C. aff. pristinus due to low sample size, our surveys and observations of habitat decline for this taxon highlight a need for increased management. Silviculture is prevalent throughout its range and many streams we visited had increased sedimentation and loss of interstitial spaces. Although not explicitly tested, it is likely that C. aff. pristinus has similar habitat requirements to C. pristinus, indicating silviculture may be detrimental to population persistence (Withers and McCoy 2005;Rohrbach and Withers 2006). This is further supported by our repeated failures to detect individuals at several historical sites within streams surrounded by silviculture and indicates local population loss. We propose that additional surveys in concert with a genetic assessment be conducted to determine the current status of C. aff. pristinus. In addition, we recommend surveys that use an occupancy modeling framework to determine if this taxon is extirpated from historic sites or if individuals are simply difficult to detect.
In conclusion, we recommend continued recognition of Cambarus pristinus as a species of conservation concern, with elevated risks of extinction. We propose a similar ranking for C. aff. pristinus. Both have small ranges, are patchily distributed, occur in low densities, and have been potentially extirpated from parts of their ranges. Furthermore, we demonstrate that C. pristinus shows patterns of genetic variation consistent with small, fragmented populations which would be susceptible to stochastic and long-term threats. The data we have provided here can serve as a guide to developing effective conservation strategies for both species. However, more work is needed to understand the basic biology and ecology of each distinct lineage to guide specific conservation efforts that may reverse their statuses.