A genomic perspective on the conservation status of the endangered Nashville crayfish (Faxonius shoupi)

The Nashville crayfish (Faxonius shoupi, Hobbs 1948) was federally listed as an endangered species in 1986 due to its limited distribution in the Mill Creek watershed; this waterway lies in the rapidly developing Nashville basin and has experienced habitat degradation due to agricultural run-off, contamination, and urban development. Recovery efforts, including dam removal and restoration of riparian zones, have improved conditions in Mill Creek and F. shoupi has increased in numbers and recolonized extirpated stream segments. However, a history of demographic bottlenecks and restricted gene flow may have negatively impacted the long-term recovery of this species. A recently discovered population of F. shoupi in a disjunct segment of the Lower Tennessee River at the Pickwick Tailwater may provide an additional source of genetic variation. Uncertainty surrounding the origins of the Pickwick population and its taxonomic relationship to F. shoupi in Mill Creek raises questions about the conservation and management implications of this population. We used mitochondrial sequencing and SNP genotyping to assess genetic variation and connectivity of F. shoupi in the Mill Creek drainage and to investigate the taxonomy and demographic history of the newly discovered population at Pickwick. We found substantial genetic variation and evidence of connectivity for samples throughout Mill Creek for both mitochondrial and genome-wide SNPs. Our results also suggest a recently severed connection between crayfish in Pickwick and Mill Creek. Unique mitochondrial haplotypes and SNP variation in the Pickwick population highlight the need for prioritizing this population in future conservation and management plans for this species.


Introduction
The Nashville Crayfish (Faxonius shoupi) is the only federally protected crayfish in the state of Tennessee. Listing of this species in 1986 was enacted in response to the limited distribution of F. shoupi in Middle Tennessee and the loss of quality habitat and ongoing development in the Nashville Basin (USFWS 1986). The known distribution of F. shoupi from the time of its discovery in 1938-1939(first identified as Cambarus propinquus sanborni, Fleming 1939 until 2019, has primarily been limited to Mill Creek (MC) and its tributaries (Davidson and Williamson counties; Fig. 1). There have also been reported occurrences of F. shoupi in Big Creek (Elk River drainage), the South Harpeth Creek, and in Richland Creek (Cumberland River drainage) (Bouchard 1974), however the Big Creek and Richland Creek records may have been errors (Miller et al. 1990). The greatest threats to the persistence of F. shoupi in MC has been poor water quality and siltation from agriculture, competition with the invasive crayfish, and rapid urbanization within Williamson and Davidson counties [U.S. Fish and Wildlife Service (USFWS) 2018]. Restoration projects have been implemented in MC in an effort to improve water quality, restore channel stability, enhance aquatic habitat, and restore riparian buffer function. Recent surveys have shown that F. shoupi can now be observed in high numbers, even in heavily developed areas; as a result, the USFWS proposed 1 3 to remove this species from the federally endangered species list (USFWS 2018). However, disruption of gene flow due to habitat disturbance and impoundments as well as loss of genetic variation from recent bottlenecks can still inhibit the long-term adaptive potential and resiliency of this species.
Repeated sampling efforts of F. shoupi populations throughout the MC drainage show that populations were negatively impacted by urbanization and development (Bouchard 1984;Barrociere 1986, USFWS 2018. Specimen numbers from surveys in 1986 were less than a third of the numbers collected from the same sites in 1969 (Hobbs et al. 1969). Also, in the Sevenmile tributary to MC, F. shoupi appears to have been replaced by an undescribed species later identified as F. durelli (Bouchard 1972). More recent surveys in MC and its tributaries indicate that F. shoupi has responded favorably to ongoing habitat restoration efforts. Recolonization surveys have shown that F. shoupi is able to quickly recolonize habitat following local extirpations (Carpenter 2002). Despite recovery in numbers, these demographic bottlenecks can still have negative fitness consequences by depleting standing genetic variation, fixing deleterious alleles, and reducing adaptive potential to environmental change.
The recent discovery of a disjunct population of F. shoupi in the Lower Tennessee River at the Pickwick Tailwater (PW) in Hardin County, Tennessee suggests that this species historically had a wider distribution. This newly discovered PW population is separated from F. shoupi at MC by nearly 700 river kilometers. Hundreds of crayfish survey efforts, spanning over five decades, have been conducted within the intervening tributary systems and have not recovered additional specimens of F. shoupi (Hobbs & Shoup 1942, O'Bara et al. 1985, Schuster et al. 2008, Couch & Schuster 2020. However, crayfish surveys are typically performed in wadeable sections of rivers, either by seining or snorkeling (Larson & Olden 2016); few surveys targeting crayfish have been performed in deeper portions of the Cumberland and Tennessee Rivers. It is possible that populations of F. shoupi occur in deeper sections of these rivers and have not been detected. The restricted distribution of F. shoupi is likely due to a lack of suitable habitat in the mainstem of the Cumberland River and Lower Tennessee River. Faxonius shoupi inhabits streams with an abundance of loose slab rock for cover habitat (Barrociere 1986). Tributaries to the ascending arm of the Tennessee River (Pickwick Tailwater/Kentucky Reservoir) lie within the Western Highland Rim and Coastal Plain provinces which exhibit gravel and sand substrates and would not be suitable for F. shoupi.
The PW population was initially identified as F. shoupi based on similarities in traditionally diagnostic morphological characters; first form male gonopod, cephalothorax morphology and chelae morphology are a close match to F. shoupi specimens from MC populations. However, the PW population does exhibit a unique color pattern that differs from color patterns in MC populations; PW color patterns were characterized by pale vermiculations on a Fig. 1 Map showing all sites sampled for F. shoupi, F. placidus, and F. durelli. Sites are color coded by species as indicated in the legend contrasting darker pigmented region of the cephalothorax, compared to the small pale speckles restricted to punctations observed on the MC specimens (Fig. 2). Mill Creek and PW sites differ in stream order and depth. Mill Creek is a second order stream with a maximum depth of less than two meters. In contrast, river width at the PW site is more than 300 m and crayfish were collected from a depth of nearly four meters. Importantly, the substrate of the Tennessee River below Pickwick dam is consistent with F. shoupi habitat in Mill Creek and is characterized by an abundance of bedrock with large loose slab rocks. 1 3 common at this location and easily collected, indicating a robust reproducing population.
The relationship of the F. shoupi PW to populations in MC has important implications for management strategies and conservation status of this species. Comparisons of gonopod and chelae morphology suggest that the newly discovered population at PW is a disjunct population of F. shoupi; discovery of an additional, viable population of F. shoupi would indicate that this species is not as vulnerable as previously thought. However, genetic evidence is needed to investigate the taxonomic relationship of the PW population to F. shoupi, as traditional morphological traits frequently contradict molecular-based species assignments (Taylor and Knouft 2006). It is also unclear if the PW population of F. shoupi is the result of a recent, anthropogenic release or if this population represents an older native population that was only recently discovered. Genetic tools can be used to test hypotheses regarding colonization history. Anthropogenic mediated introductions, such as a bait-bucket release, would be established by a small number of propagules; this founder event would leave a detectable genetic signature (Johnson et al. 2011). Specifically, population bottlenecks quickly remove rare alleles from a population resulting in an excess of heterozygosity relative to expectations for a population in mutation-drift equilibrium. Also, this bottleneck event would reduce genome-wide genetic variation relative to source populations of the species (Nei et al. 1975, Allendorf 1986). Information on the taxonomic relationship and origins of the crayfish at PW is needed to effectively incorporate this population into management strategies for this species.
We used both mitochondrial sequencing and genotypeby-sequencing datasets to investigate population structure and demographic histories of F. shoupi in MC and to investigate the history and taxonomic status of the recently discovered population in the Lower Tennessee River downstream of Pickwick Dam. Our specific objectives were to assess levels of genetic variation within populations of F. shoupi, test for genetic signatures of recent population declines, estimate levels of gene-flow throughout the MC drainage, and reconstruct the relationship between the MC and PW sites. Phylogenetic reconstructions and a Bayesian multispecies coalescent model were applied to estimate population divergence times between the two disjunct populations.

Collections
A total of 76 tissue samples were obtained for molecular analysis by TWRA from crayfish morphologically identified as F. shoupi from three sites within MC, Tennessee and 29 samples were obtained from one site in the Tennessee River at Pickwick Tailwater (Table 1, Fig. 1). Collections at PW were conducted by TWRA scuba divers operating from an anchored dive boat at depths ranging from approximately 2-4 m and a river width of 308 m. Crayfish were collected by overturning loose slab rocks and quickly grabbing specimens with gloved hands. The specimens were placed in mesh holding bags, brought to the surface where they were processed for data before being released. Sampling sites within the MC drainage were selected to span the currently known distribution of the species and included upper MC (MC-U), middle MC (type locale; MC-M) and lower MC (MC-L) ( Table 1; Fig. 1). Tissues from outgroup species were collected from the Green River in Buffalo River drainage (F. durelli, type locality), Barrens Fork in Collins River drainage (F. placidus, BF population) and Sinking Creek in the Lower Cumberland-Old Hickory Lake drainage (F. placidus, LC population, type locality) ( Table 1). Tissue samples from F. shoupi were obtained from leg sections separated at autonomizing joints (approximately 5 mm) and individuals were immediately released on site. Crayfish are known to undergo autotomy of appendages and regeneration initiates after ~ 6 days, depending on the intermolt period (Durand 1960). Complete specimens for outgroups were preserved in 70% ethanol and will be deposited as voucher

Mitochondrial sequencing
Genomic DNA was extracted from leg tissues using the EZNA Tissues DNA Kit (Omega Bio-tek) except that DNA was eluted in water in the final step. The 5' end of the mitochondrial gene cytochrome oxidase I (COI) was amplified using universal primers LCO1490 (GGT CAA CAA ATC ATA AAG ATA TTG G) and HCO2198 (TAA ACT TCA GGG TGA CCA AAA AAT CA) (Folmer et al. 1994). Conditions for polymerase chain reactions (PCR) were as follows: initial denaturation step of 5 min at 95 °C followed by 35 cycles of 15 s at 95 °C, 15 s at 54 °C, 45 s at 72 °C. This program ended with a final extension of 10 min at 72 °C. PCR products were cleaned prior to cycle sequencing reactions by exonuclease I/shrimp alkaline phosphatase (New England Biolabs) and used for bi-directional Sanger sequencing on an ABI 3730 automated sequencer (MCLAB). Sequence chromatograms were imported and visualized using SEQUENCHER version 5.2 (Gene Codes Corp.). Sequences were aligned using the ClustalW alignment algorithm (Thompson et al. 1994) as implemented in Bioedit (Hall 1999). Alignments were refined by eye and protein-coding genes were examined for stop codons to check for miscalled bases and pseudogenes.

Mitochondrial haplotype diversity
Estimates of genetic variation within sampled populations at the mitochondrial COI gene, including genotypic diversity (H GD ) and nucleotide diversity (π M ) were performed in Arlequin 3.5.2 (Excoffier et al. 2010). Mitochondrial haplotypes from all F. shoupi individuals and outgroups were trimmed to an equal length of 632 base pairs and used in a reconstruction of a minimum spanning haplotype network as performed by the package pegas (Pardis 2010) in the R-studio suite (Racine 2012). Structuring of genetic variation for mitochondrial COI haplotypes was investigated using Analysis of Molecular Variance (AMOVA) as performed by Arlequin (Excoffier et al. 1992). The hierarchical analysis investigated the proportion of molecular variance due to difference between drainages (i.e. MC populations and PW; F ST ), variation between populations within drainages (F SC ), and variation within sampled populations (F CT ).

Phylogenetic analyses of mitochondrial haplotypes
Phylogenetic reconstructions of mitochondrial haplotypes were estimated using both Bayesian and maximum likelihood (ML) optimality criteria. Both analyses assumed the HKY-Gamma model of evolution (Hasegawa et al. 1985); this was selected as the best-fit model under Bayesian Information Criteria in as performed by MEGA X (Kumar et al. 2018). A Bayesian phylogenetic reconstruction was performed using MrBayes version 3.2.7a on the CIPRES Science Gateway (Miller et al., 2010). The Markov chain Monte Carlo (MCMC) algorithm ran for 10,000,000 generations, sampling every 1000 generations. Two independent runs were performed and the resulting trees were combined after the deletion of a burnin (first 10% of trees). A majorityrule consensus tree was generated and nodal support was estimated by posterior probabilities. A ML analysis was performed using MEGA X and the tree with the highest log likelihood was retained. Nodal support was estimated by bootstrap analysis with 1000 replicates.

Genotyping by sequencing
Genotyping by sequencing (GBS) was used to identify and genotype SNPs across all four F. shoupi sites, F. placidus LC and F. placidus BF (Table 2) following the protocol described by Elshire et al. (2011). Briefly, GBS library construction used the same DNA extractions that were used for mitochondrial Sanger sequencing. Extracted DNA was standardized to 5 ng/μl and quantified using the Quantifluor dsDNA System (Promega). Standardized DNA was digested with the restriction enzyme ApeKI (New England Biolabs) and barcoded adapters were ligated to cut restriction sites. Each individual was labelled with a unique barcode (4-8 bp in length). Ligated DNA was pooled and PCR amplified using primers complementary to ligated adapters. The resulting PCR product was cleaned using the Qiaquick PCR purification kit and the distribution of PCR fragment lengths was measured using an Agilent 2100 BioAnalyzer (Agilent Technologies). The resulting library was cleaned using the AxyPrep Mag bead purification system (Agilent) and sequenced on an Illumina NexSeq 500 (Illumina Inc.) with the 75 bp, single-end read chemistry.

Genotyping and SNP filtering
Single Nucleotide Polymorphisms were identified and filtered using the de novo map STACKS pipeline (Catchen et al. 2013). In ustacks, the minimum number of raw reads needed to form a stack (allele) was set to m = 3 and the number of mismatches between alleles within individuals M = 2. A catalog of loci was constructed in cstacks where n = 2 mismatches were allowed between sample tags. The value of n was increased from the default value of n = 1 following recommendations by Paris et al. (2017) suggesting that n = M and an increase in n allows for fixed differences across populations. The populations function was used to estimate population summary statistics and to generate input datasets for downstream analyses. One source of potential bias in GBS datasets is mutations that occur at restriction enzyme cut-sites; the accumulation of substitutions in restriction sites can influence population genetic inferences by causing allele drop-out and biasing datasets to more conserved regions of the genome (Gautier et al. 2013;Silliman et al. 2021). We investigated the influence of filtering parameters for missing data on the retention of loci and variable sites using the filtering flags in the populations function. Two separate datasets were used for the final analyses that were generated using different filtering parameters. Estimates of contemporary effective population sizes and tests for bottlenecks used the p1r70 dataset, where SNPs were required to be present in 70% of individuals in a single population. This dataset was used to maximize the number of loci within populations instead of shared loci across populations. Estimates of within population genetic variation and analyses of population structure were performed using the p2r60 dataset, so that only SNPs present in at least two out of four of the F. shoupi populations were included in these analyses. Maximum heterozygosity was set at 0.60 to remove homeologous loci. Tests for Hardy-Weinberg equilibrium (HWE) were performed in Plink 2.0 (Chang 2015) and loci with p-values below 0.01 in any single populations were removed from further analyses.

Genetic variation at SNP loci
Estimates of genetic variation including mean observed heterozygosity (H o ), expected heterozygosity (He), and nucleotide diversity at variable sites (π N ) were estimated in the program populations within the STACKS workflow using the p2r60 SNP dataset. We estimated the average allelic richness (A RN ) and private allelic richness (P RN ) using a rarefaction approach to correct for unequal sample size across sites as implemented in the software ADZE (Szpiech et al. 2008).
Due to the loss of shared loci across species we did not include F. placidus LC and F. placidus BF in rarefaction analyses.

Demographic histories
Contemporary effective populations size (N e ) estimates for all populations of F. shoupi and F. placidus were performed on SNP datasets (p1r70) using the LD method (Waples and Do 2008), as implemented in NeEstimator v2.01 (Do et al. 2014). The LD method measures the strength of association of alleles at independent loci to infer the degree of drift as a function of N e . A threshold of 0.02 for the allelic frequency was chosen in order to reduce bias resulting from low-frequency variants. Confidence intervals (95%) were estimated using parametric bootstrap. Linkage disequilibrium based estimates of N e were adjusted for physical linkage using the equation by Waples et al. (2016).
We used the estimated chromosome number from the congeneric F. virilis for N e adjustments (Chr = 200, Mlinarec 2016).
Rapid declines in population size can also leave a signature on patterns of allelic variation within populations due to a loss of low frequency alleles. We tested for evidence of recent population bottlenecks from SNP datasets (p1r70) using the program BOTTLENECK (Piry et al. 1999). A test for excess heterozygosity (compared with mutationdrift equilibrium) was performed using the 1-tailed Wilcoxon signed rank test (α = 0.05), which detects historical  Cornuet and Luikart 1996). Simulation parameters were set for 1000 iterations and assumed an infinite allele model (IAM). We also tested for a mode-shift in allele frequency distributions to detect a genetic signature of more recent bottlenecks (Luikart et al. 1998

Population structure
Pairwise population differentiation at SNP loci was summarized using the p2r60 dataset using pairwise F ST ' and Φ ST as implemented in the program Populations in the STACKS workflow. F ST ', an analogue of F ST , is a standardized measure of differentiation that is independent of the amount of genetic variation that is present within populations (Meirmans 2006). The Φ ST statistic is an AMOVA-based measure of population differentiation that incorporates information regarding the molecular distance between haplotypes (Excoffier 1992). A Bayesian assignment test was implemented in the software STRU CTU RE v. 2.3.4 (Pritchard et al. 2000) on all F. placidus and F. shoupi populations (p2r60); F. durelli was not included due to small sample size. We performed ten independent runs for each value of K (1-6) assuming the admixture model. Markov chain Monte Carlo (MCMC) settings included two independent MCMC chains, each with a burnin of 50,000 iterations followed by 100,000 iterations. The optimal valued of K was examined by plotting the average natural log of the probability of the data for K = 1-6 (Ln Pr(X|K); Pritchard and Wen 2003). We also investigated the optimal value of K using the ΔK method as implemented in Structure Harvester; however, due to small standard deviation across replicates for K = 3, the ΔK method could not be performed (Evanno 2005;Earl 2012). Permuted membership assignment coefficients from the ten independent runs were estimated using Clumpp v. 1.1.2 (Jakobsson and Rosenberg, 2007). Structure bar plots were generated using STRU CTU RE PLOT (Ramasamy et al. 2014).
Discriminant analysis of principal components (DAPC) was performed in R using the package 'adegenet' on the complete SNP dataset (p2r60) that included all sampled individuals and also on the dataset that included only F. shoupi individuals from MC and PW (Jombart et al. 2010). This method determines the number of populations or groups without prior knowledge of sample origins. Discriminant analysis of principal components differs from Bayesian clustering in that it does not assume that loci are unlinked and that populations are randomly mating. The optimal number of clusters was determined using BIC as performed by the find.clusters function.

Divergence time estimation
Divergence times estimates for the split between the F. shoupi MC sites and PW were estimated from GBS sequence data using the Bayesian Multispecies Coalescent model with and without introgression (MSC and MSci, respectively) as implemented in the A00 model in BPP v.4.2 (Yang 2015;Flouri et al. 2020). The input tree for the BPP analysis was based on the topology of the mtDNA gene tree. Two different models were used for divergence time estimates. The first model assumed no migration after initial population divergence between F. shoupi populations at MC and PW. The second model allowed bidirectional migration between the two drainages. Each analysis of 50,000 MCMC generations was ran twice from different starting seeds with a burn-in period of 5000 using algorithm 0 (default fine-tuning parameter, e = 2) and an estimated heredity (a = 4, b = 4); this gave consistent parameter estimates between replicate runs and ESS values > 1000 for all parameters. We consider speciation probability values > 0.95 as strong support for all speciation events. The theta prior (3, 0.04), tau prior (3,0.2) phi prior (1,1) finetune automatic update. The posterior probabilities of the model with migration (M1) and without migration (M2) were compared using the stepping-stones method as implemented in the BPPR package (https:// github. com/ dosre islab/ bppr, see also, and Rannala and Yang 2017). The stepping-stone analysis used 16 independent MCMC chains for each model (50,000 generations per chain). Mean divergence time parameters were converted to time in years in the BPPR package using a genome-wide mutation rate estimated from decapod shrimp in the genus Alpheus (μ = 2.64 × 10 -9 ; Silliman et al. 2021).

Mitochondrial haplotypes
Sequencing of the mitochondrial COI gene identified 28 unique haplotypes in a 669 bp alignment across all sampled populations and outgroups (Table 2). Faxonius shoupi sampled from MC-M and PW had the greatest number of haplotypes and haplotypic diversity. The COI sequence alignment used in MSN analysis was trimmed to the shortest sequence in the alignment, for an equal length of 632 bp. Trimming of sequences eliminated several variable sites so that the number of unique haplotypes was reduced from 28 to 21. The resulting MSN network showed four common mitochondrial haplotypes shared at high frequency across all sampled populations of F. shoupi 1 3 (Fig. 3). All haplotypes within F. shoupi were connected by at most four mutational steps. Haplotypes from F. placidus LC, F. placidus BF, and F. durelli were separated from F. shoupi haplotypes by a minimum of 33 steps. Phylogenetic reconstructions based on Bayesian and ML criteria resulted in nearly identical topologies. Both criteria recovered a well-supported monophyletic clade that included all haplotypes from the four sampled F. shoupi populations (Fig. 4). Haplotype sequences were not structured by populations; the 18 haplotypes belonging to F. shoupi formed an unresolved polytomy in both Bayesian and ML reconstructions.
The results of the AMOVA for COI haplotypes (Table 3) showed that, while the majority of genetic variation was found within sampled populations (98.4%), some structuring of haplotype diversity could also be attributed to differences between the populations in MC and PW (3.16%). None of the variance in mitochondrial haplotype diversity could be attributed to differences between sampled sites within the MC watershed.

SNP genotyping
We genotyped a total of 128 individuals from four populations of F. shoupi and two populations of F. placidus, and one population of F. durelli. (Tables 1,2). Three individuals with low sequence coverage (less than 500,000 sequenced tags) were removed from further analysis. The average number of reads across retained individuals was 3,931,286 sequenced tags per individual. The denovo map STACKS pipeline generated 2,031,321 loci, with a mean coverage of 14.3 across individuals and the mean number of sites per locus was 82. The influence of filtering parameters, including the minimum number of populations (p) and the proportion of sequenced individuals (r), greatly influenced the number of genotyped SNPs. A total of 701 loci (433 SNPs, Supplement 1) were retained in the dataset that included all F. shoupi samples from both MC and PW. Although, this is a modest number of SNPs in comparison to SNP panels obtained from GBS datasets where whole genome sequence is available for alignment, multiple simulation-based studies indicate that the resulting SNP panel is sufficient for detecting fine-scale genetic structure of biological significance. With a samples size of only ten individuals, a dataset consisting of 40 independent SNP loci can detect fine-scale populations structure (F ST < 0.01) with a power of 0.93 (Willing et al. 2012;Morin et al. 2009). The mean number of sampled individuals per locus, number of retained loci, variant sites and polymorphic sites for different filtering parameters are listed in Supplement 1.

Within population genetic variation
Estimates of within population genetic variation including Het Obs , Het Exp , and π N at SNP loci showed that the F. shoupi at PW had the highest levels of genetic variation (Table 2). Within the MC drainage, the MC-U population was more variable than downstream sites MC-M and MC-L. Heterozygosity was similar for F. placidus LC and for populations of F. shoupi. Faxonius placidus BF had the lowest estimates of genetic variation overall. Rarefaction estimates of the average A RN and P RN for PW and MC-U populations were similar. Downstream sites within MC had lower estimates of A RN and P RN when accounting for differences in sample size (Supplement 2).

Effective population size and bottlenecks
Estimates of effective population sizes (N e ) based on the LD method ranged from 19 to 41; F. shoupi MC-M had the highest estimate of N e and MC-U had the lowest estimates. Estimates of N e for F. shoupi PW and F. placidus LC were unreliable (N e = ∞), likely due to a lack of signal in the dataset (Table 4). Tests for historical bottlenecks indicated that  Table 1 several populations had experienced recent declines in numbers. The mode-shift test showed a shift in the spectrum of allele frequencies in four of six populations. All MC populations of F. shoupi and the F. placidus LC population showed evidence of a recent bottleneck. Allele frequency spectrums for F. shoupi PW and F. placidus BF were normal. The onetailed Wilcoxon test for excess heterozygosity was not significant for any of the examined populations (p = 1.00 for all populations).

Population differentiation SNPs
Estimates of pairwise population differentiation showed little evidence of genetic structuring among populations of F. shoupi (Table 5). Pairwise F ST ' and Φ ST estimates between populations of F. shoupi were less than or equal to 0.001 and 0.01, respectively. Populations of F. placidus showed moderate levels of differentiation from F. shoupi. The average F ST ' and Φ ST estimates for F. placidus LC/F. shoupi populations were 0.057 and 0.078, respectively; F. placidus BF/F. shoupi averaged 0.018 and 0.033, respectively.
Results from the Bayesian assignment analyses suggested that K = 3 provided the best fit for the data as determined  Table 1) are indicated at tips with haplotype copy number shown in parentheses. Numbers at nodes indicate posterior probability support (below line) and ML bootstrap support (above line) greater than 50%   Fig. 5 Results of Bayesian cluster analysis as performed by STRU CTU RE. Bar plot indicates assignment probabilities for 131 individuals based on the SNP loci generated from GBS. Bar plots shown are for K = 2 through K = 5. A value of K = 3 was found to be the best-fit model for the data as determined by the mean Ln P [D]. Population IDs are defined in Table 1 performed a separate analysis which only included the four sampled F. shoupi populations as suggested by Vaha et al. (2007); however, no further subdivision was found (results not shown). The BIC curve for the DAPC clustering analysis that included F. shoupi and F. placidus indicated that two (K = 2) was the optimal number of explanatory clusters (Fig. 6). Cluster 1 included all populations of F. shoupi as well as the F. placidus BF population. Cluster 2 only included individuals form the F. placidus LC population. Results from the BIC analysis that only included F. shoupi individuals (MC and PW sites) did not show evidence of fine-scale population subdivision. The divergence time estimates for the split between F. shoupi at MC and PW were similar for models with migration (M2) and without migration (M1). The estimate for model M1 was 69,638 YBP (Fig. 7) and M2 was 65,971 YBP (95% C.I. 32.828-100.400). The marginal log likelihood values were greater for model M1 than for model M2 (Supplement 3) and the relative posterior probability for M1 was 1, indicating a greater likelihood for the model with no post-divergence gene flow.

Discussion
The Nashville crayfish has demonstrated remarkable genetic resiliency despite a history of habitat loss, impoundment, and contamination in the MC watershed. Our results suggest that MC populations have recently experienced demographic declines; however, substantial genetic variation is indicated by both mitochondrial haplotypes and genomewide SNP datasets. Genetic connectivity among populations has been maintained and may reflect this species' capacity for dispersal and rapid colonization (Carpenter 2002). Our molecular investigation of the newly discovered PW population shows that this population shares a recent connection with F. shoupi in the MC watershed. Despite this similarity, the PW population also possesses unique genetic variation in both the mitochondrial and nuclear genome, improving the long-term adaptive outlook for this species. Results from this study provide the first glimpse of population-level genetic variation and structure in F. shoupi and give insights into the genetic consequences of the management history of this species.

Genetic variation and demographic history
Allele frequency distributions for genome-wide SNPs indicate that MC populations have experienced recent bottlenecks, likely resulting from habitat disturbance related to urbanization. The mode-shift test detected an upward shift in the distribution of allele frequencies, indicative of a recent population bottleneck, at all three MC sites. Despite evidence of recent population declines, estimates of genetic variation at mitochondrial haplotypes across MC populations were still high compared to estimates obtained from other at-risk North American crayfish species. Comparable surveys of mitochondrial gene diversity (also referred to as haplotype diversity) in the genus Cambarus averaged 0.360 and 0.297, for sampled populations of C. jezerinaci and C. parvoculus, respectively (Thoma and Fetzner 2008). Gene diversity estimates for our MC samples averaged 0.826; estimates obtained here are similar to estimates from native populations of the invasive Procambarus clarkii (H GD averaged 0.844; Yi et al. 2018). Rapid recovery following bottleneck events can limit the loss of genetic variation (Nei et al. 1975). Experimental evidence has shown that F. shoupi can recover quickly following population depletion; their ability to rapidly colonize and recover from loss of numbers may have mitigated the loss of genetic variation that can accompany demographic bottlenecks (Carpenter 2002). This study is among the first to use GBS derived SNPs to assess patterns of population-level genetic variation in crayfish, so it is not possible to place our estimates of heterozygosity in the context of similar studies. Furthermore, bioinformatic parameters can influence estimates of genetic variation, confounding cross-study comparisons. We did see consistent patterns of genetic variation across sampling sites within MC. Individuals sampled from the MC-U site demonstrated higher levels of genetic variation than downstream sites (MC-M and MC-L). This same pattern was observed after correcting for unequal sample sizes (Supp. 2), although MC-U had the lowest estimates of private allelic richness of the three sampled sites. Genetic variation measures for both mitochondrial and nuclear markers were slightly higher for the F. shoupi population at PW than for any of the sampled F. shoupi populations in MC.
The F. shoupi PW population was more variable and harbored unique genetic variation not present in F. shoupi populations from MC sites. The F. shoupi PW population possessed the greatest number of unique mitochondrial haplotypes (5 unique haplotypes; Fig. 4) and high estimates of genetic variation for genome-wide SNPs. It is possible that the Pickwick population has retained genetic variation because it has not experienced the level of urbanization, contamination and siltation that has been documented along MC. Despite impoundments along the stretch of the Tennessee river, water quality and flow regimes near the Pickwick Tailwater have been maintained. Greater demographic stability is also supported by the mode-shift test that did not detect a shift in the allele frequency spectrum for PW.
Estimates of effective population sizes across F. shoupi populations at MC were markedly low; the average N e across F. shoupi populations in MC was 23.7. It is difficult to establish a baseline for comparison of N e estimates in crayfish as few empirical studies have investigated contemporary N e from multi-locus molecular data. A recent study using microsatellite genotypes to estimate N e across populations of Cambarus pristinus recovered much larger N e values than our estimates for F. shoupi. Point estimates of N e for populations of C. pristinus ranged from 137 to 1348 using the same LD method applied here (Grubb 2020). Molecularbased estimates of N e are often much lower than contemporary census numbers; demographic bottlenecks will reduce the N e /N as will variation in reproductive success. Little is known about the reproductive life-history of F. shoupi. Generally, crayfish demonstrate r-selected life histories with females carrying between 200 and 400 eggs that hatch every spring (Bouchard 1976). Organisms with high fecundity typically experience high mortality in early life-stages so that very few parents successfully pass down genetic material to the next generation. This high variance in reproductive success reduces the N e /N (termed the Hedgecock effect; Hedrick 2005). The relationship of the Ne/N ratio as a function of the variance in the number of progeny per parent (V k ) is estimated as follows (Wright 1938): In animals that exhibit this "sweepstake reproductive success", Ne/N can be reduced to a small fraction (<< 0.01; Hedgecock and Pudovkin 2011).

Population structure
Molecular results indicate that F. shoupi populations throughout MC have maintained geneflow despite documented habitat disturbances and impoundments. Historically, drainage from MC into the Cumberland River was impeded by a series of three dams. As of 2018, all dams were removed, restoring connectivity for the entire length of MC downstream to the Cumberland River (Cumberland River Compact 2018). Both mitochondrial haplotypes and genome-wide SNPs suggest that genetic connectivity has been maintained along MC. Analysis of molecular variance indicated that variation in mitochondrial haplotypes could not be explained by differences in sampling sites within the MC watershed. Pairwise measures of population differentiation (F ST and ϕ ST ) were less than 0.01 for all comparisons; for reference, values of F ST less than 0.05 are generally interpreted as little structuring of genetic differentiation (Hartl and Clark 1997a, b). The F. shoupi PW population shared genetic variation with populations in MC for both mitochondrial haplotypes MC; however, we did detect some molecular signatures of recent isolation. The PW population possessed genetic variation that was not present in any of the sampled MC sites. Our Bayesian MSci analysis suggested that connectivity between these sites was severed less than 100,000 YBP.
One explanation for the disjunct distribution of F. shoupi is that the MC and PW are relictual populations of what historically was a more widespread distribution throughout the Cumberland and Tennessee watersheds. Reports of F. shoupi in the Elk and Harpeth River Systems lend support to the hypothesis of a now extant, larger distribution of F. shoupi in Middle Tennessee (Barrociere 1986). Populations at Big Creek and the South Harpeth were presumably extirpated. The population at Richland Creek may have been replaced by the closely related F. placidus, which is the predominant, stream dwelling crayfish in this basin (Bouchard 1984). The disjunct distribution of F. shoupi is shared by other freshwater taxa; specifically, there are about 36 mussel species demonstrating disjunct distributions that are endemic to the Tennessee and Cumberland drainages (Haag and Cicerello 2016). This shared assemblage, long recognized as unique and first termed by Ortmann (1918) as the Cumberlandian Fauna, may be the result of the two drainages previously being isolated from the Mississippi and Ohio Rivers but connected to each other (Galloway et al. 2011). The presence of F. shoupi in the lower Tennessee River may represent a remnant population of this former connection.
Alternatively, the newly found PW population could be the result of a recent, human-mediated introduction as a baitbucket release, however our molecular results do not support this. Introduced populations are typically initiated by a small number of founders and, as a result, demonstrate patterns of genetic variation indicative of a population bottleneck (i.e. an excess of heterozygosity and shifts in the distribution of allele frequencies; Fitzpatrick et al. 2012). This pattern was recently found in the signal crayfish (Pacifastacus leniusculus), where surveys of mitochondrial haplotypes revealed a loss of gene diversity and a reduction in the number of haplotypes in introduced populations compared to populations in their native range (Larson et al. 2012). The PW populations was the only F. shoupi population that did not have a mode-shift in its allele frequency distribution signifying a recent population bottleneck. Colonizers of introduced populations will pass along a small subset of genetic variation from their source, so that genetic variation in the introduced site will be a subset of what is observed in the native population range. The PW population did not reveal any evidence of reduced genetic variation for either mitochondrial haplotype or genome-wide SNP surveys. Furthermore, the PW population possessed unique genetic variation that could not be traced back to any of our MC samples.
The inclusion of samples from multiple F. placidus populations highlights the need for a thorough molecularbased evaluation of this group. Faxonius placidus, as currently described, is a widespread species associated with the Cumberland, Barren, Tennessee, and Lower Ohio rivers (Poly and Wetzel, 2003). Our phylogenetic reconstructions based on mitochondrial haplotypes and pairwise measures of differentiation based on SNPs reveal that F. placidus is paraphyletic with respect to F. shoupi; the F. placidus BF population is more closely related to a monophyletic F. shoupi than it is to the F. placidus population at the type locality (population LC). Integrative species delimitation methods are currently being used to investigate the relationship between populations of F. placidus throughout Tennessee and their relationship with F. shoupi (Hildreth et al. in prep).

Implications for management
Faxonius shoupi has served as the flagship species for the protection and restoration of MC. Prior to the discovery of the PW population, F. shoupi was widely considered to be the only known species whose distribution was restricted entirely to the Nashville basin (Barrociere 1986). Habitat in MC has been compromised by loss of riparian buffer zones, contamination from sewage or chemical spills, and siltation (Miller and Hartfield 1985). Federal listing of F. shoupi has resulted in improvements to this watershed that have not only benefited F. shoupi, but have also restored habitat to a number of co-distributed freshwater fauna. Despite evidence of demographic bottlenecks, F. shoupi have maintained moderate levels of genetic variation and connectivity throughout MC; the ability of this species to rapidly colonize and expand in numbers may have mitigated some of the loss of genetic variation that can accompany sustained demographic declines.
The discovery of a second F. shoupi population in the Tennessee River improves the long-term outlook for this species. While it is not possible to conclusively determine the origins of the PW population, our results demonstrate that the PW population is genetically variable and should be prioritized for conservation planning along with F. shoupi in the MC watershed. Although this species is considered to be a recovery success, the distribution of F. shoupi is still restricted in its range and vulnerable to declines should quality habitat be lost. The disjunct population at PW should be included in management planning for F. shoupi; this site is at risk of degradation associated with urbanization and agriculture, the same factors that threaten MC populations (Bizwell and Mattingly 2010). Potential threats unique to the PW include the introduction of invasive species and pathogens from the Tombigbee waterway, pollution from upstream sources, and changes to flow regimes from impoundments. Further sampling along the mainstem of the Middle Tennessee River is needed to determine the full extent of the range of F. shoupi; the discovery of F. shoupi at PW indicates that this species may inhabit deeper riverine habitat than previously thought. Clarification of the full distribution of the PW population will be useful when reviewing the conservation status of this species.