Establishing conservation units to promote recovery of two threatened freshwater mussel species (Bivalvia: Unionida: Potamilus)

Abstract Population genomics has significantly increased our ability to make inferences about microevolutionary processes and demographic histories, which have the potential to improve protection and recovery of imperiled species. Freshwater mussels (Bivalvia: Unionida) represent one of the most imperiled groups of organisms globally. Despite systemic decline of mussel abundance and diversity, studies evaluating spatiotemporal changes in distribution, demographic histories, and ecological factors that threaten long‐term persistence of imperiled species remain lacking. In this study, we use genotype‐by‐sequencing (GBS) and mitochondrial sequence data (mtDNA) to define conservation units (CUs) for two highly imperiled freshwater mussel species, Potamilus amphichaenus and Potamilus streckersoni. We then synthesize our molecular findings with details from field collections spanning from 1901 to 2019 to further elucidate distributional trends, contemporary status, and other factors that may be contributing to population declines for our focal species. We collected GBS and mtDNA data for individuals of P. amphichaenus and P. streckersoni from freshwater mussel collections in the Brazos, Neches, Sabine, and Trinity drainages ranging from 2012 to 2019. Molecular analyses resolved disputing number of genetic clusters within P. amphichaenus and P. streckersoni; however, we find defensible support for four CUs, each corresponding to an independent river basin. Evaluations of historical and recent occurrence data illuminated a generally increasing trend of occurrence in each of the four CUs, which were correlated with recent increases in sampling effort. Taken together, these findings suggest that P. amphichaenus and P. streckersoni are likely rare throughout their respective ranges. Because of this, the establishment of CUs will facilitate evidence‐based recovery planning and ensure potential captive propagation and translocation efforts are beneficial. Our synthesis represents a case study for conservation genomic assessments in freshwater mussels and provides a model for future studies aimed at recovery planning for these highly imperiled organisms.


| INTRODUC TI ON
Population genomics has significantly increased our ability to make inferences about microevolutionary processes (e.g., gene flow, genetic drift, population structure, selection, and mutation) through the use of thousands of genome-wide markers (Allendorf et al., 2010;Luikart et al., 2003). In comparison with microsatellite or Sanger sequencing methodologies (typically 10-20 markers), genotype-by-sequencing (GBS) methods (e.g., restriction enzymebased sequencing approaches) have substantially more power to resolve population dynamics given the vast number of molecular markers. Advancements in genomic technologies have been proven useful for improving protection and recovery of imperiled species (Allendorf et al., 2010;Coates et al., 2018;Funk et al., 2012Funk et al., , 2019, even in nonmodel organisms without established genomic resources (Ellegren, 2014).
Conservation biologists and natural resource managers use information on resiliency, redundancy, and representation to inform management efforts for imperiled organisms (Smith, Allan, et al., 2018;Wolf et al., 2015). However, determining "the three Rs" relies on the a priori designation of population units across the species' range. The U.S. Endangered Species Act (ESA) classifies population units as distinct population segments (USFWS & NMFS, 1996), which are often delineated using natural or manmade barriers. The rationale for this approach is these structures theoretically promote genetic drift between extant populations. However, molecular studies have shown this may or may not be the case (e.g., Clemento et al., 2009;Grummer & Leaché, 2017;Hoffman et al., 2017), which can have substantial impacts on recovery planning. Thus, it is becoming increasingly evident that the use of molecular data is a more thorough approach to delineate population units for imperiled species. Specifically, the implementation of genomic methodologies to define the boundaries of conservation units (CUs), or population units for conservation (Funk et al., 2012), is a necessity in order to infer the status of populations within imperiled species (e.g., genetic diversity, population demographics, population trends).
Freshwater mussels (Bivalvia: Unionida) are a group of aquatic bivalves comprised of approximately 840 species (Graf & Cummings, 2007), but also represent one of the most imperiled groups globally (Lopes-Lima et al., 2018). Widespread alteration to freshwater ecosystems has led to systemic declines in abundance of both common and rare freshwater mussel species (Haag & Williams, 2014;Vaughn & Taylor, 1999). These declines stem from a suite of traits that make freshwater mussels sensitive to ecosystem state change, including filter feeding, limited locomotive capabilities, and an obligate parasitic life cycle that requires codependence with hosts to complete metamorphosis (Haag, 2012;Randklev et al., 2019;Williams et al., 1993). Widespread decline has major ecological ramifications considering the loss of freshwater mussel biodiversity can negatively impact ecosystem function of freshwater systems (Vaughn, 2018;Vaughn et al., 2008). These factors have led conservation biologists and natural resource managers to prioritize freshwater mussels as a group of greatest conservation concern and a keystone for guiding freshwater ecosystem restoration (Ferreira-Rodríguez et al., 2019;Haag & Williams, 2014;Lopes-Lima et al., 2018;Williams et al., 1993).
In this study, we use molecular data and available survey information to delineate CUs and investigate distributional trends of two highly imperiled freshwater mussel species: Potamilus amphichaenus and Potamilus streckersoni. Potamilus amphichaenus is endemic to the Sabine, Neches, and Trinity River drainages in eastern Texas (Howells et al., 1996), while the distribution for P. streckersoni is restricted to the Brazos River drainage in central Texas (Smith et al., 2019). As members of Potamilus, both species are presumed to be host specialists, with glochidia only transforming on Aplodinotus grunniens (Bosman et al., 2014;, a common molluscivorous fish distributed throughout Gulf of Mexico drainages (Page & Burr, 2011 conservation genetics, DPHS, endangered species, population genomics, recovery range of P. amphichaenus and P. streckersoni, and (4) discuss our findings in terms of their impact on future conservation, management, and recovery practices.

| Sampling design and DNA extraction
We sampled P. amphichaenus and P. streckersoni from focal drainages in Texas (i.e., Brazos, Neches, Sabine, and Trinity). Potamilus amphichaenus was collected from 24 localities (Neches = 8, Sabine = 8, Trinity = 8) and P. streckersoni was collected from 22 localities in the Brazos River drainage (Table 1). Outgroups were not included considering multiple phylogenetic studies have resolved P. amphichaenus and P. streckersoni as sister species with strong support (Smith et al., 2019. Genomic DNA was extracted from fresh mantle clips using the Gentra PureGene extraction kit following manufacturer protocol (Qiagen; Hilden, Germany). High molecular weight DNA was ensured by visualizing isolations on a 1% agarose gel stained with GelRed Nucleic Acid Stain (Biotium), and the purity of each isolation was quantified using a NanoDrop™ (Thermo Fisher Scientific).

| Library preparation, sequencing, and filtering
Sequencing for SNP genotyping was performed using DArTseq™ (DArT Pty Ltd). Briefly, DNA samples were processed in digestion and ligation reactions using the Pstl-Sphl restriction enzyme combination following Kilian et al. (2012) but replacing a single Pstl-compatible adaptor with two different adaptors corresponding to each restriction enzyme. The Pstl-compatible adapter was designed to include the Illumina flow cell attachment sequence and the reverse adapter contained the flow cell attachment region and Sphl-compatible attachment sequence (see Elshire et al., 2011). Fragments were amplified using the following thermal cycling conditions: 1-min initial denaturation at 94°C, then 30 cycles of denaturation (20 s, 94°C), annealing (30 s, 58°C), and extension (45 s, 72°C), with a final extension of 7 min at 72°C. The resulting products were subsequently sequenced using 75 base pair (bp) single-end sequencing on an Illumina Hiseq-2500.
Reads were demultiplexed using the individual-specific barcode sequence ligated to the samples. Reads were then processed using proprietary DArT Pty Ltd analytical pipelines (see Georges et al., 2018;Wenzl et al., 2004). At first pass, poor quality reads were filtered (barcode Phred score < 30, read Phred score < 10). Retained sequences were truncated to 69 bp and aggregated using the DArT fast clustering algorithm with a Hamming distance threshold of 3 bp.
Error correction was performed using a proprietary algorithm which corrects low-quality bases (Phred score < 20) with a corresponding high-quality singleton tag (Phred score > 25). Identical sequences were then collapsed, and SNPs were called using DArTsoft14.
Additional filtering was performed in the dartR package v 1.1.11 ) using R v 3.5.1 (R Core Team, 2018 following similar methodologies as Georges et al. (2018). Loci with less than 100% reproducibility, a statistic that measures the concordance of a genotype between a minimum of 27 (~30% of genotyped individuals) technical replicates (see Wenzl et al., 2004), and greater than 30% missing data were removed. Subsequently, individuals with greater than 30% missing data and minor allele frequencies (<0.05) were removed from the dataset. We then removed secondary SNPs by retaining the SNP with the highest degree of polymorphism at each locus. Lastly, we excluded loci found to be in linkage disequilibrium using the R package SNPRelate v 1.16.0 (Zheng et al., 2012).

| Genotype-by-sequencing analyses
A total of 91 individuals of P. amphichaenus (Sabine = 9, Neches = 33, Trinity = 15) and P. streckersoni (Brazos = 34) were used in all GBS analyses (Table 1). Phylogenomic inference was performed on a concatenated alignment of all full loci using BEAST v 2.6.2 . Before the analysis, we used bModelTest v 1.2.1 (Bouckaert & Drummond, 2017) to estimate the best nucleotide substitution model for the analysis using the "transitionTransversionsplit" set of models. A strict molecular clock was paired with a GTR + I + G model of nucleotide substitution. The analysis consisted of 10 9 MCMC generations logging every 1,000 trees with an initial 50% burn-in. Tracer v 1.7 (Rambaut et al., 2018) was used to evaluate the trace log to ensure proper burn-in and convergence of all parameters (ESS > 200). A maximum clade credibility tree was estimated in TREEANNOTATOR v 2.6.0 .
For all downstream GBS analyses, individuals were binned into four groups based on drainage of capture: (a) P. amphichaenus from the Sabine, P. amphichaenus from the Neches, P. amphichaenus from the Trinity, and P. streckersoni from the Brazos. Estimates of genetic diversity and population substructuring were conducted using a custom R script utilizing numerous packages (available at https:// github.com/chase smith 15/Potam ilus_EE). We calculated allelic richness (AR), observed heterozygosity (H o ), expected heterozygosity (H e ), and inbreeding coefficient (F IS ) for each group using diveRsity v 1.9.90 (Keenan et al., 2013). For AR and F IS , 999 bootstrap replicates were used with a critical value of 0.05. Fixed alleles and private alleles for each group were determined in dartR. To visualize genetic structure data relative to geographical distribution, we performed a principal coordinate analysis (PCoA) in dartR. We also performed a discriminant analysis of principal components (DAPC) in adegenet v 2.1.0 (Jombart, 2008;Jombart & Ahmed, 2011) on the first two PCs and DA eigenvalues. Additionally, DAPC predicts group membership probability for each sample, and we compared clustering results to membership designations as stated above. The best-fit number of clusters (K) was determined using k-means.
We used traditional and model-based methodologies to identify patterns of population structure in our genomic data. First, we calculated pairwise F ST values using the R package StAMPP v 1.5.1   TA B L E 1 (Continued) delimitation (BFD) was used to assess support for each of the clustering scenarios with 2lnBF > 10 representing significant support (Kass & Raftery, 1995).

TA B L E 1 Molecular material examined in this study
We used Neestimator v 2.1 (Do et al., 2014) to estimate effective population size (N e ) and the number of effective breeders (N b ) for the four groups based on drainage of capture. The molecular coancestry (Nomura, 2008) and linkage disequilibrium (Waples, 2006) methods were used to estimate N b and N e , respectively. Per developers' recommendations, a critical value of 0.05, singleton exclusion, and jackknifed confidence intervals were used to reduce inflated estimates and address potential linkage.

| Mitochondrial data generation and analyses
We amplified and sequenced the mitochondrial (mtDNA) gene  (Table 1). Phylogenetic inference was performed on the ND1 alignment using BEAST. Before the analysis, we used bModelTest to estimate nucleotide substitution models for each codon position using the "transitionTransversionsplit" set of models. A strict molecular clock was specified for each codon position for both bModelTest and the standard BEAST analysis, and the analyses consisted of 5 * 10 7 and 10 8 MCMC generations, respectively. Tracer was used to ensure convergence of all parameters for both analyses, and a maximum clade credibility tree was estimated in TREEANNOTATOR.
For all subsequent mtDNA analyses, individuals were binned into four groups based on drainage of capture: (a) P. amphichaenus from the Sabine, P. amphichaenus from the Neches, P. amphichaenus from the Trinity, and P. streckersoni from the Brazos. We used DnaSP v 6.12 (Rozas et al., 2017) (Leigh & Bryant, 2015) with the default epsilon value set at 0 (Bandelt et al., 1999). Complete deletion was used for missing data at any given nucleotide position.

| Distribution and abundance estimates
To estimate relative abundance and distributional trends of P. amphichaenus and P. streckersoni throughout their known ranges, we compiled available or generated data from freshwater mussel surveys that detected live specimens in the Brazos, Neches, Sabine, and Trinity River drainages (Ford et al., ,2009(Ford et al., , , 2014(Ford et al., , , 2016Randklev et al., 2011Randklev et al., , 2017Randklev et al., , 2020Smith et al., 2019). For select sites where survey effort was reported (e.g., survey time, number of surveyors), we estimated abundance using catch per unit effort (CPUE), which is calculated by dividing the total number of live individuals by the total person-hours. Resulting estimates of CPUE were plotted to visualize distribution and relative abundance of both species.

| Genotype-by-sequencing analyses
A total of 91 individuals of P. amphichaenus (57) and P. streckersoni  however, F IS was also lower in P. streckersoni than populations of P. amphichaenus. Fixed and private alleles were much higher in comparisons of P. streckersoni and populations of P. amphichaenus, and fixation was low across populations of P. amphichaenus (Table 2).
PCoA showed 4 distinct groupings: P. amphichaenus from the Sabine, P. amphichaenus from the Neches, P. amphichaenus from the Trinity, and P. streckersoni (Figure 3). The first axis defined 33.5% of the variance (P. amphichaenus from P. streckersoni), and the second axis described 3.6% of the variance (populations of P. amphichaenus). kmeans supported the best number of clusters as K = 4, aligning with groupings depicted by PCoA. All clusters were strongly supported by DAPC, with all individuals having 100% membership probability when compared to designations by drainage.
Pairwise F ST was much larger between P. streckersoni and all populations of P. amphichaenus than comparisons within P. amphichaenus populations (Table 2). Pairwise F ST was also found to be larger between P. amphichaenus from the Trinity and P. amphichaenus from the Neches or Sabine ( Table 2). The fastSTRUCTURE analysis using all individuals resolved K = 3 as the best value to explain structure and maximize likelihood in our GBS dataset: P. amphichaenus from the Sabine + Neches, P. amphichaenus from the Trinity, and P. streckersoni ( Figure 4b). For the P. amphichaenus dataset, K = 2 (Sabine + Neches, Trinity) was resolved as the best value to maximize likelihood, while K = 3 (Sabine, Neches, Trinity) was supported as the best value to explain structure. For P. streckersoni, K = 1 was supported as the best value to explain structure and maximize likelihood. The TESS3 analysis using all data supported K = 2 selected by cross-entropy plot, representing P. amphichaenus and P. streckersoni (Figure 4a). TESS3 analyses using species-specific datasets supported K = 2 for P. amphichaenus (Sabine + Neches, Trinity) and K = 1 for P. streckersoni.
The geographic distribution of K = 4 (Figure 4c), which aligned with the best K to explain structure in P. amphichaenus using model-and nonmodel-based approaches, depicted genetic clusters spanning the drainages of capture: Sabine, Neches, Trinity, and Brazos ( Figure 5).
Convergence of path sampling analyses in SNAPP was supported by all steps having ESS values greater than 200. Bayes factor delimitation marginally supported K = 3 as the most likely clustering scenario and rejected K = 2, but BFD could not reject K = 4 ( Table 3).
Estimates of N e and N b for each group are reported in Table 4.
Contrary to genetic diversity statistics, N b was higher in P. streckersoni (9.2) when compared to populations of P. amphichaenus (5.2-7.0).
Estimates of N e were higher in P. streckersoni (5,243.7) when compared to P. amphichaenus from the Neches (538.3) but depicted wide confidence intervals ( Table 4). Estimates of N e for P. amphichaenus from the Sabine and Trinity were infinite, presumably due to low sample size.

| mtDNA analyses
A total of 116 individuals were sequenced for ND1 (900 bp): P. amphichaenus from the Sabine (9), P. amphichaenus from the Neches (38), P. amphichaenus from the Trinity (23), and P. streckersoni (46). All ND1 sequences are deposited on GenBank, and accession numbers can be found in Table 1. The following substitution models were selected by bModelTest: ND1 1st codon position-TPM1 + I + G, ND1 2nd codon position-TrN, and ND1 3rd codon position-HKY + G. The ND1 reconstruction resolved three strongly supported clades representing P. amphichaenus from the Neches + Sabine, P. amphichaenus from the Trinity, and P. streckersoni ( Figure 6). Convergence of the analysis was supported by ESS values for each parameter greater than 200, and a 10% burn-in was deemed appropriate by Tracer.
Estimates for h, Hd, k, and π for each group are reported in P. amphichaenus from the Sabine + Neches, P. amphichaenus from the Trinity, and P. streckersoni ( Figure 6). There was limited divergence between each population of P. amphichaenus accompanied by haplotype sharing between the Sabine and Neches populations.  Figure 7). Potamilus streckersoni appears to be extirpated from the Brazos River upstream of Lake Waco, but downstream of this  Table 6.

| Distribution and abundance estimates
Regressions depicted a general increasing trend in DPHS for each species with respect to time ( Figure 8a); however, survey effort was shown to significantly influence DPHS (r² = 0.509; p = 0.001914) indicating that general increasing trends may be due to increased sampling effort rather than population expansion (Figure 8b,c).

| D ISCUSS I ON
Genetic data are effective at delineating clusters that can be used to inform management and restoration efforts, but our molecular analyses supported disputing number of genetic clusters within P. amphichaenus and P. streckersoni ranging from K = 2-4 ( Figure 4). This is not unexpected because it is an unrealistic expectation that samples will conform to assumptions of each model and the notion of a "true" K may not be applicable to natural populations (e.g., Janes et al., 2017;Jombart et al., 2010;Pritchard et al., 2000;Raj et al., 2014). However, we do not see incongruence as an issue because ecological characteristics and life-history information combined with molecular data can be used to identity appropriate spatial units for conservation (e.g., Funk et al., 2012;Palsbøll et al., 2007). Below, we outline the clustering scenarios supported by our molecular data, possible biogeographic, ecological, and life history-driven explanations for differing clustering scenarios, a methodological approach to choose the "best" K, and how we integrate molecular findings and survey information to guide recovery planning for P. amphichaenus and P. streckersoni.

| Assessing genetic structure in Potamilus amphichaenus and Potamilus streckersoni
The resolution of two genetic clusters (i.e., K = 2; Figure 4a) appears to be driven by species boundaries, with the two clusters TA B L E 2 Summary of genetic diversity statistics for Potamilus amphichaenus and Potamilus streckersoni and pairwise F ST between lineages derived from GBS data representing P. amphichaenus and P. streckersoni. The bias toward clustering algorithms resolving K = 2 has been well discussed (e.g., Janes et al., 2017), and in our case, could be due to the fact that the majority of the variation in our molecular data could be explained by separating the two species (33.5%; Figure 3). Species boundaries have been shown to mask population structure within species when using model-based approaches (e.g., Warner et al., 2015), and we addressed this issue by creating species-specific datasets for P. amphichaenus and P. streckersoni. These analyses, along with traditional nonmodel-based approaches, showed clear evidence for genetic structuring within P. amphichaenus that align with biogeographic hypotheses (Haag, 2010;Figures 2-6; Tables 2 and 4).
The resolution of three genetic clusters from GBS data (i.e., K = 3: P. amphichaenus from the Neches + Sabine, P. amphichaenus from the Trinity, and P. streckersoni) aligned well with results from mtDNA (Figures 4b and 6). Limited divergence between the Neches and Sabine is expected since the two drainages share an embayment, plus the two systems were likely interconnected during the last glacial maxima (Blum & Hattier-Womack, 2009;Blum et al., 2013) leading to an increased possibility of recent gene flow. However, GBS data did show limited, but diagnosable, molecular differentiation between the Neches and Sabine systems, which supported the possibility of K = 4 (Figures 2-3 and 4c; Table 2). Limited genetic divergence between the Neches and Sabine drainages could be reminiscent of repeated isolation and secondary contact, or also driven by human-mediated extirpation events. Given similar phylogeographic patterns seen in congeners ) and life-history characteristics, it is defensible that recent fluctuations in sea level likely led to limited differences at GBS markers rather than by human-mediated extirpation events in the lower Sabine.
Several genetic clustering scenarios were resolved by our molecular methods, as outlined above. To address this issue, we used SNAPP to estimate the marginal likelihood of clustering scenarios resolved by molecular data. Although GBS data were able to diagnose the Neches and Sabine populations, BFD marginally supported three genetic clusters (Table 3)

Potamilus streckersoni Potamilus amphichaenus
Brazos Trinity Sabine N eches Neches + Sabine, P. amphichaenus from the Trinity, and P. streckersoni. We could not reject the model that separated the Neches and Sabine as independent clusters (K = 4; Table 3). Though our GBS dataset clearly contained more variability than our mtDNA dataset, BFD aligned with results supported by mtDNA and emphasizes the utility of mtDNA markers in investigating intraspecific relationships in freshwater mussels. Congruent patterns in mtDNA and nuclear markers may not hold across all freshwater mussel species (e.g., Chong et al., 2016), but our findings support that previous population genetic research using mtDNA should not be discredited.
Although GBS data were able to diagnose all independent river drainages, our data did not depict evidence of intradrainage population structure. Unlike most freshwater mussels, Potamilus species typically reach sexual maturity within 1 year and have a comparatively short lifespan (5-9 years) (Haag, 2012;Haag & Rypel, 2011), which theoretically would lead to an increased possibility of genetic drift. Despite these characteristics, our results are similar to other population genetic studies in freshwater mussels, where hypothetical genetic drift caused by dams constructed within the past 100 years was not detected (e.g., Elderkin et al., 2007;Hoffman et al., 2017). and Neches drainages (Table 2), which have plausibly been separated since the last glacial maxima (Blum & Hattier-Womack, 2009;Blum et al., 2013).
Freshwater mussels are nearly completely reliant on their host for dispersal (Haag, 2012;Vaughn, 2012), and patterns of genetic structuring can be heavily influenced by host use (e.g., Karlsson et al., 2014;Wacker et al., 2019). In the case of P. amphichaenus and P. streckersoni, both species are presumed to be host specialists, with glochidia only transforming on A. grunniens (Bosman et al., 2014;, a wide ranging and mobile species that has been documented to travel over eighty kilometers in several days (Hansen et al., 2020) and up to several hundred kilometers during migration events (Funk, 1957). Dispersal capabilities of A. grunniens could therefore explain the lack of divergence within drainages, considering panmixia was likely historically present before impoundment of river systems. In addition to host use, other life-history

Potamilus amphichaenus (Sabine)
Potamilus amphichaenus (Trinity) characteristics such as fecundity may also explain the lack of genetic structuring within drainages. Fecundity is often positively correlated with the retention of genetic diversity in natural populations (Ellegren & Galtier, 2016;Romiguier et al., 2014), and both P. amphichaenus and P. streckersoni have high annual fecundity (>1,000,000 ;. In addition to a mobile host, high fecundity could be contributing to the lack of structuring within drainages, given the increased possibility of high recruitment and retention of neutral genetic variation.

| Establishing conservation units for Potamilus amphichaenus and Potamilus streckersoni
The demand for water and hydropower, coupled with ongoing changing climate, has led to obvious reductions in streamflow, alteration of sediment transport, increased salinity levels, and exacerbated highly variable seasonal hydrological fluctuations (i.e., flood and drought) throughout the ranges of P. amphichaenus and P. streckersoni (Cañedo-Argüelles et al., 2013;Philips et al., 2004;Randklev et al., 2011Randklev et al., , 2013Randklev et al., , 2017Randklev et al., , 2019Wellmeyer et al., 2005). These impacts combined with changes in land use have likely contributed to the widespread decline of both species throughout significant portions of their historical range (Ford et al., 2014(Ford et al., , 2016Howells et al., 1996;Randklev et al., 2011Randklev et al., , 2017Smith et al., 2019). Our estimates of abundance and contemporary distribution for P. amphichaenus and P. streckersoni support this observation (Figure 7), but DPHS estimates were less clear and depicted an increase in distribution relative to time for both species (Figure 8; Table 6). In recent years, however, there has been a resurgence of sampling effort that has resulted in rediscovery of presumed extirpated populations of numerous mussel species (e.g., Holcomb et al., 2015;Johnson et al., 2016;Randklev et al., 2010Randklev et al., , 2012. We show a similar scenario within our focal species as survey effort significantly influenced DPHS (p = 0.001914), which leads us to hypothesize that the observed increasing trend in distribution for P. amphichaenus and P. streckersoni is a result of increased survey effort rather than range expansion.
The establishment of CUs within P. amphichaenus and P. streckersoni is the logical next step in recovery planning, and we found defensible evidence for the delineation of four genetically diagnosable CUs: P. amphichaenus from the Neches, P. amphichaenus from the Sabine, P. amphichaenus from the Trinity, and P. streckersoni. We do recognize that divergence between P. amphichaenus from the Neches and Sabine is limited (Table 2) and some analyses supported the two drainages as a single CU (Figures 4 and 6; Table 3); however, the genetic distinctiveness of the two drainages (Figures 2 and 3) warrants the recognition of independent CUs. Although the lack of mtDNA diagnosability has been used to justify mixing of CUs (Moritz, 1994), we caution the use of P. amphichaenus from the Neches in recovery efforts in the Sabine. The adaptive significance of the observed genetic distinctiveness between the Neches and Sabine is uncertain, but recovery efforts in these drainages should rely on stock sources within each CU to avoid genetic consequence. Our findings provide valuable information for natural resource managers, especially considering brood stock selection for recovery planning is likely less stringent than previously conceived.  Kilometers Currently, populations of freshwater mussels for conservation and management practices are typically defined as geographic management units (GMUs), or a unit that is geographically or otherwise identifiable by man-made and natural barriers (e.g., USFWS, 2018USFWS, , 2020. In the case of P. amphichaenus and P. streckersoni, there is no molecular support for the subdivision of drainages into GMUs; however, we were unable to include material from multiple stream stretches where the species are presumed extirpated. The lack of molecular diagnosability does not discredit the use of GMUs for management and recovery practices because it is an unrealistic expectation that stressors and habitat suitability will be uniform throughout the range of freshwater mussel species (e.g., Randklev et al., 2019;Strayer et al., 2004;Vaughn & Taylor, 1999).
Contemporary estimates of abundance for both species show stark contrasts with respect to geography (Figure 7), and rather than manage species at a drainage level, we encourage natural resource managers to integrate our CUs within existing frameworks.
For example, population densities in the Neches drainage are variable (Figure 7), and translocation of P. amphichaenus from the lower Neches River (i.e., below B.A. Steinhagen Lake) to augment populations in the upper Neches River (i.e., above B.A. Steinhagen Lake) could be an effective recovery option to improve resiliency of the CU with limited genetic consequence. While CUs offer protection at the drainage level, delineation of GMUs within each CU will allow for more robust investigations of population characteristics such as relative abundance, age-class structure, and threats to long-term sustainability, all of which should be considered in recovery planning.

| CON CLUS ION
Our study represents a model for population genomic assessments of freshwater mussels and provides information for natural resource managers in the development of conservation and recovery strategies for P. amphichaenus and P. streckersoni. Given the wide diversity of host use in freshwater mussels, it is an unrealistic expectation that other imperiled species that co-occur with P. amphichaenus and P. streckersoni will depict similar patterns of genetic structuring.
Thus, there remains a critical need for robust molecular investigations to support recovery planning for many imperiled species. As genomic resources are developed, the identification of potentially adaptive loci through RNA sequencing and whole genome resequencing may improve brood stock selection and species recovery.

ACK N OWLED G M ENTS
The authors wish to thank Kaitlin Havlik for assistance in the labora-

DATA AVA I L A B I L I T Y S TAT E M E N T
All data and materials in this study are freely available in data repositories. Raw reads are deposited in the SRA (BioProject ID: PRJNA663379), and sample accession numbers can be found in