A Comparative Phylogeographic Approach to Facilitate Recovery of an Imperiled Freshwater Mussel (Bivalvia: Unionida: Potamilus inﬂatus )

: North American freshwaters are among the world’s most threatened ecosystems, and freshwater mussels are among the most imperiled inhabiting these systems. A critical aspect of conservation biology is delineating patterns of genetic diversity, which can be di ﬃ cult when a taxon has been extirpated from a signiﬁcant portion of its historical range. In such cases, evaluating conservation and recovery options may beneﬁt by using surrogate species as proxies when assessing overall patterns of genetic diversity. Here, we integrate the premise of surrogate species into a comparative phylogeographic framework to hypothesize genetic relationships between extant and extirpated populations of Potamilus inﬂatus by characterizing genetic structure in co-distributed congeners with similar life histories and dispersal capabilities. Our mitochondrial and nuclear sequence data exhibited variable patterns of genetic divergence between Potamilus spp. native to the Mobile and Pascagoula + Pearl + Pontchartrain (PPP) provinces. However, hierarchical Approximate Bayesian Computation indicated that the diversiﬁcation between Mobile and PPP clades was synchronous and represents a genetic signature of a common history of vicariance. Recent ﬂuctuations in sea-level appear to have caused Potamilus spp. in the PPP to form a single genetic cluster, providing justiﬁcation for using individuals from the Amite River as a source of brood stock to re-establish extirpated populations of P. inﬂatus . Future studies utilizing eDNA and genome-wide molecular data are essential to better understand the distribution of P. inﬂatus and establish robust recovery plans. Given the imperilment status of freshwater mussels globally, our study represents a useful methodology for predicting relationships among extant and extirpated populations of imperiled species. phylogeographic techniques to evaluate range-wide genetic diversity within P. inﬂatus and sympatric congeners P. fragilis and P. purpuratus using multi-locus sequence data. Next, we integrate the premise of surrogate species into a comparative phylogeographic of exchange


Introduction
Due to anthropogenic alterations to the environment, the world is losing species at rates comparable to mass extinctions during major transitions of geological time periods [1,2]. North American freshwaters are among the world's most threatened ecosystems [3], and freshwater mussels (Bivalvia: Unionida) are among the most imperiled groups of organisms inhabiting these systems with 65% of all recognized species considered to be of conservation concern [4][5][6]. Several inherent biological characters (e.g., limited locomotive capabilities in many species, extreme sensitivity to pollutants, obligate parasitism, and filter feeding) have disproportionately impacted mussels in anthropogenically framework to hypothesize genetic relationships between extant and extirpated populations of P. inflatus by characterizing genetic structure in co-distributed congeners with similar life histories to better inform conservation and recovery planning.

Taxon Sampling
We examined genetic diversity from co-distributed members of Potamilus native to the Mobile, Pascagoula, Pearl, and Pontchartrain drainages (Table 1; Figure 1). A total of 103 individuals were examined in this study (Table 1), and more details on the specimens and collections are available on ScienceBase (https://doi.org/10.5066/P9Q3CFL5) and Johnson and Smith (Under Review). Mantle tissue clips from vouchered individuals in public museums were used to extract genomic DNA using the Qiagen PureGene DNA extraction kit with the standard extraction protocol (Qiagen, Hilden, Germany). We amplified and sequenced two mitochondrial (mtDNA) loci commonly used in freshwater mussel phylogenetic studies: a partial portion of cytochrome c oxidase subunit 1 (CO1) and NADH dehydrogenase subunit 1 (ND1). For all P. inflatus and a subset of P. fragilis and P. purpuratus, we sequenced three nuclear (nDNA) loci: the commonly used internal transcribed spacer 1 (ITS1), and two protein-coding loci fem-1 homolog C (FEM1) and UbiA prenyltransferase domain-containing protein 1 (UbiA). A subset of individuals representing the geographic range of P. fragilis and P. purpuratus were chosen for the additional nDNA loci due to the high prevalence of multiple copies at ITS1 and low genetic diversity at FEM1 and UbiA (Table 1; Figure 1). We utilized two recently developed primer sets from Johnson and Smith (Under Review) to amplify FEM1 and UbiA based on data generated in phylogenetic studies using the recently developed anchored hybrid enrichment probe set Unioverse [28,43,44]. Primers for all loci and thermal cycling conditions are reported in Table 2.  PCRs were conducted using a 25 µL mixture of the following: molecular grade water (9.5 µL), MyTaq TM Red Mix (12.5 µL; Bioline), primers (1.0 µL each) and DNA template (100 ng). Products were sent to Molecular Cloning Laboratories (McLAB, South San Francisco, CA, USA) for bi-directional sequencing on an ABI 3730. Geneious v 10.2.3 was used to assemble contigs and edit chromatograms [45], and sequences were aligned in Mesquite v 3.31 [46] using MAFFT v 7.311 [47]. Loci were aligned independently using the L-INS-i method in MAFFT and translated into amino acids to ensure absence of stop codons and gaps. Incomplete codons at each terminal end were removed. The total number of individuals included for each locus are as follows: CO1-102, ND1-103, FEM1-29, UBiA-29, and ITS1-31. Novel GenBank accessions for this study were as follows: CO1: MT662002-MT662099; FEM1: MT669773-MT669799; ITS1: MT661766-MT661792; ND1: MT669647-MT669745; and UbiA: MT669746-MT669772 (Table 1).
Coalescent-based approaches have been repeatably criticized to delimit populations and not species [52], including in freshwater mussels [53][54][55]. Here, we use the Bayesian coalescent-based model STACEY [56] in BEAST v 2.6.2 [57] to define genetic clusters in our molecular dataset for downstream analysis. STACEY allows for the inclusion of individuals with missing data, so we included all available data for the 5 loci in the analysis. Potamilus spp. were binned by drainage of capture, and we allowed the model to freely assign drainages to appropriate clusters. A substitution model for each locus alignment was determined using ModelFinder, a strict molecular clock was set at 1.0 for CO1, and clock rates for the four additional loci were estimated by STACEY. The Epi Tree prior was used as the species tree prior with a collapse height of 0.0001. Our analyses executed 10 9 generations and logged every 5000 trees with an initial 10% burn-in. Effective sample size (ESS) was ensured using Tracer v 1.7 [57], and the most likely number of clusters was calculated by SpeciesDelimitationAnalyser (SpeciesDA) v 1.8.0 [56] with a collapse height of 0.0001, a 1.0 simcutoff, and an initial 10% burn-in (2000 trees).
To estimate divergence times among well-supported clusters, we used the Bayesian coalescent-based model *BEAST [58] in BEAST. We chose a coalescent approach to account for concatenation methods, which typically overestimate the divergence times across species trees [59,60]. Similar to STACEY, *BEAST allows for the inclusion of individuals with missing data and all available data for the five loci were included in the analysis. For each species, individuals were grouped based on the most likely clusters resolved by STACEY: (1) Mobile; and (2) Pascagoula + Pearl + Pontchartrain (herein referred to as PPP). A strict molecular clock and an HKY model of nucleotide evolution was fit to each locus to better match priors for comparative phylogeographic analyses (see below). The substitution rate for CO1 was set to 2.56 × 10 −9 ± 0.6 × 10 −9 substitutions per site per year [61], and substitution rates were estimated for the four additional loci. Yule process was used as the species tree prior paired with a piecewise linear and constant root population size model. The analysis was run for 1.5*10 9 MCMC generations sampling every 5000 generations and a 10% burn-in. Effective sample size (ESS) was ensured using Tracer v 1.7 [57], and a maximum clade credibility tree was created using TreeAnnotator v 2.5 [57].

Phylogeographic Analyses
To visualize genetic divergence with respect to geographic distribution, we created a median joining haplotype network [62] for each of the three Potamilus spp. independently in PopART 1.7 [63] with the default epsilon value set at 0. Additionally, an analysis of molecular variance (AMOVA) was conducted for each species independently in PopART to further evaluate the geographic distribution of genetic diversity. Each analysis was performed on a concatenated alignment of CO1 and ND1, and missing data in both PopART analyses was handled using complete deletion. To further assess genetic variation within Potamilus spp. with regard to geography, we calculated DNA sequence divergence between clusters of Potamilus spp. using uncorrected pairwise genetic distances in MEGAX [64]. Partial deletion was used to handle missing data in MEGAX calculations. For haplotype networks, species were grouped by drainage and groups for all other analyses were as follows: P. fragilis from the Mobile and Pearl + Pontchartrain, P. inflatus from the Mobile and Pontchartrain, and P. purpuratus from the Mobile and PPP.

Comparative Phylogeography
We tested for simultaneous divergence between clusters of Potamilus spp. defined by STACEY under a hierarchical Approximate Bayesian Computation (hABC) approach as implemented in the PyMsBayes package [65]. Specifically, we tested if divergence between Mobile and PPP clusters of P. fragilis, P. inflatus, and P. purpuratus was synchronous. PyMsBayes implements a modified version of msBayes [66] that specifies a Dirichlet-process prior (dpp) to compare fit of empirical data to simulated data under user-informed priors [14]. We used dpp-msbayes to test for synchronous divergence between Mobile and PPP clusters of Potamilus spp. using alignments from all available loci. We used results from our *BEAST divergence time analysis to guide prior selection for dpp-msbayes as follows: the concentration parameter [1000, 0.00141] in which there was prior probability for one, two, or three divergence events, population size (θ) [1, 0.0005], and divergence times (τ) [1, 0.01].
To allow dpp-msbayes to freely explore different divergence scenarios, we allowed the model to estimate independent parameters for each species (θ parameter = 012) and the number of divergence events (τ classes = 0). Transition-transversion rate of the HKY substitution model was estimated for each alignment independently using IQ-TREE. Our dpp-msbayes run performed a total of 10 7 simulations with 10,000 standardizing samples and reported every 20,000 simulations. We retained the 1000 simulations with the best fit to empirical data to estimate posterior probability (PP) values for each divergence scenario. To measure support for the number of divergence events, Bayes factors were measured using twice the difference of −ln likelihood [67].

Molecular Analyses
Five partitions and substitution models were determined by ModelFinder for phylogenetic reconstruction in IQ-TREE: TN + F + I for mtDNA codon 1 and nDNA codon 3, TN + F + I for mtDNA codon 2 and nDNA codon 2, K3Pu + F + G4 for mtDNA codon 3, F81 + F for nDNA codon 1, and K2P + I for ITS1. All species-level relationships had full support (BS = 100) and the only two major nodes that were not strongly supported (i.e., BS ≥ 95) were the PPP clade of P. fragilis (BS = 94) and the Mobile clade of P. purpuratus (BS = 92; Figure S1). All three taxa were resolved as monophyletic with P. inflatus, the sister to P. fragilis and P. purpuratus, aligning with findings in a previous phylogenetic study [28].
Substitution models determined by ModelFinder for locus alignments in the STACEY analysis were: HKY + I for CO1, HKY + I for ND1, JC for FEM1, F81 + I for UbiA, and K2P + I (=K80 + I) for ITS1. Convergence of the analysis was supported by all parameters having ESS values > 200, and all nodes were strongly supported (PP ≥ 95). SpeciesDA supported six clusters (54%): (1) P. inflatus from the Mobile, (2) P. inflatus from the Pontchartrain, (3) P. fragilis from the Mobile, (4) P. fragilis from the Pearl + Pontchartrain, (5) P. purpuratus from the Mobile, and (6) P. purpuratus from the PPP ( Figure S2). The second most likely clustering scenario supported seven clusters (18.5%), with the Pearl population of P. purpuratus recognized as a distinct cluster.
The topological reconstruction from *BEAST was congruent with IQ-TREE and STACEY topologies, and all nodes were strongly supported (Figure 2A). Mobile and PPP clusters of Potamilus spp. were resolved as monophyletic with full support (PP = 100; Figure 2A  Mean uncorrected p-distances between Mobile and PPP clusters for all species were larger than 1% and are reported in Table 3. Distance values were larger in P. inflatus (2.33%) when compared to P. fragilis (1.11%) and P. purpuratus (1.31%). Haplotype networks were concordant with phylogenetic analyses and showed clear separation between the Mobile and PPP groupings of all three Potamilus spp. (Figure 3). However, within the PPP province there was haplotype sharing between drainages in P. fragilis and P. purpuratus (Figure 3). AMOVAs indicated the majority of molecular variation was distributed between Mobile and PPP groups for all Potamilus spp. (Table 3). Molecular variance was higher within P. fragilis (19.1%) than P. inflatus (1.1%) and P. purpuratus (3.7%).  The dpp-msbayes analysis supported synchronous divergence between clusters of Potamilus spp. (Figure 2B-D). Support for a single divergence event was 55.7 PP, with the next best supported scenario of two divergence events (P. inflatus and P. purpuratus equal, and P. fragilis subsequently diverged independently) at 15.7 PP (Figure 2C,D). Similarly, Bayes factors indicated positive support for one divergence event (2lnBF = 1.7), and negative support for two (2lnBF = −0.74) and three (2lnBF = −2.19) divergence events ( Figure 2B). The overlap of confidence intervals for divergence estimates in the *BEAST analysis and dpp-msbayes further supports evidence of synchronous divergence between Potamilus spp. (Figure 2A).

Discussion
Accurate evaluations of genetic diversity are a critical component in developing effective conservation and recovery strategies. The specific goal of our study was to characterize range-wide genetic variation of P. inflatus. Given the overall rarity of the species and plausible extirpation from multiple river systems, estimating genetic relationships across the historical range of P. inflatus is completely dependent on understanding the genetic composition of closely related and co-distributed species with similar dispersal capabilities and life histories. Our comparative phylogeographic approach integrated the premise of surrogate species to predict relationships among extant and extirpated populations. Although the use of comparative phylogeography has been used to characterize genetic diversity in common and rare species within freshwater mussels [22][23][24][25][26], the use of surrogate species within a comparative phylogeographic framework to hypothesize relationships among extant and extirpated populations of imperiled species is a novel approach. Below, we discuss the evolutionary forces driving congruent patterns of genetic divergence within Potamilus spp., and how our findings may impact future conservation and recovery efforts for P. inflatus.

Patterns of Genetic Variation in Potamilus Species
Large-scale environmental change has substantial effects on communities of species and associated microbiota [14,72,73]. This is certainly the case in mussels and their hosts, as biogeographical processes are a driver of faunal structure and genetic diversity [28,29,53,[74][75][76]. Given biogeography is a critical driver of genetic variation, identifying faunal provinces is the first step toward understanding specific patterns of phylogeography [77]. Multiple attempts have been made to classify North American mussel fauna into biogeographic provinces [78][79][80][81][82], and understanding the processes that have driven faunal shifts across these regions has been integral toward understanding the evolution of the group [28,74,83]. In the case of the Mobile and PPP provinces, the drainages have been linked in hierarchical classifications of mussel diversity based on species composition [79]. Prior to our study, however, these relationships have not been tested in a molecular context. Our molecular analyses align with the hypothetical historical connection between the Mobile and PPP, as our phylogenetic and coalescent-based species delimitation analyses strongly supported Potamilus spp. in these biogeographic provinces as distinct clines. These results align with other mussel species showing genetic distinctiveness across these drainages [53,69,83,84], as well as other aquatic species [85][86][87][88][89][90].
The geological connection between the Mobile and PPP drainages has been hypothesized by numerous authors (reviewed by [91]) and a vicariance event between the two systems has likely driven the observed genetic differentiation in Potamilus spp. If a vicariance event was the causation of molecular diversification for all the species, we would expect to see similar patterns of divergence across Potamilus spp. Results from our molecular analyses, however, deviated from these expected patterns of genetic drift and showed varying levels of sequence divergence (Table 3). Specifically, genetic distance values between populations of P. inflatus were larger than those in P. fragilis and P. purpuratus (Table 3). However, it is an unrealistic expectation to assume that rates of evolution are identical between species, especially across geographically isolated populations [18,92,93]. Variable rates of molecular diversification within Potamilus spp. could be indicative of a variety of confounding variables, such as differing population demographics (e.g., population size, age structure), evolutionary processes (e.g., mutation rate, genetic drift, selection), or species-specific traits (e.g., habitat preferences) rather than multiple hypothetical vicariance events. To address this issue, we used a hABC approach to explicitly test whether divergence between Mobile and PPP populations of Potamilus spp. occurred synchronously. Our results suggest that the divergence between Potamilus spp. in the Mobile and PPP occurred simultaneously and further support previously described biogeographic provinces [79]. The causative event driving genetic differentiation between these groupings is uncertain, but additional molecular investigations in other freshwater mussels, as well as host fishes, may further elucidate the timing and patterns of faunal exchange between these two provinces.
Despite extensive geographic range within the PPP, our molecular data showed no diagnostic divergence between drainages within the province (Figure 3; Table 3). Limited genetic diversity was suspected within P. inflatus given there is only one extant population; however, the more common and wide-ranging species, P. fragilis and P. purpuratus, both showed haplotype sharing between drainages and no evidence of drainage specific structuring within the PPP (Figure 3A,C; Table 3). A signal for incomplete lineage sorting at nDNA loci is expected due to the effective population size being nearly four times that of mtDNA loci [94,95]; however, incomplete lineage sorting of mtDNA loci likely indicates relatively recent gene flow between populations. Approximately 18 Kya during the last glacial lowstand; geological evidence suggests the PPP drainages were connected [91,96], permitting gene flow to occur. Subsequent sea level rise from deglaciation began to form modern fluvial systems in the PPP [96], causing genetic isolation among contemporary populations of Potamilus spp. Given the hypothetical mtDNA mutation rates of freshwater mussels [61,97], it is an unrealistic expectation that mtDNA markers would become fixed across these drainages, and using more rapidly evolving markers (genotype-by-sequencing, whole genome resequencing) would be necessary to delineate fine-scale genomic differentiation among Potamilus spp. inhabiting these drainages or test for ongoing gene flow. However, only one extant population of P. inflatus occurs within the PPP (Amite River-Pontchartrain drainage) and it is a realistic expectation that the presumed extirpated populations of P. inflatus in the Pontchartrain and Pearl drainages would be most closely related to members of the Amite River given the patterns of genetic diversity seen in P. fragilis and P. purpuratus.

Implications on Conservation
Numerous practices have been proposed for delineating population units using genetic data for conservation and management [98][99][100][101][102], and in particular, ESA listed species have been partitioned into distinct population segments (DPSs), evolutionary significant units (ESUs), or management units (MUs) to facilitate recovery practices [101][102][103][104]. However, under the ESA, DPSs only apply to vertebrate species [103] and formal recognition of population units remains rare in freshwater mussels [105]. This is particularly concerning because information regarding population units is often required to facilitate conservation and recovery efforts [106]. In the case of P. inflatus, we observed high levels of molecular divergence at mtDNA loci. Formal recognition of ESUs are diagnosed based on reciprocal monophyly [104] and significant differences in allele frequencies at both mtDNA and nDNA loci [107]. Although individuals from the Amite River and Mobile drainage show evidence of fixation at mtDNA markers, we saw no evidence of fixation at nDNA loci, which would rule out the recognition of the two populations as ESUs. However, it is an unrealistic expectation that the nuclear loci used in this study would diagnose population units within species, and assessments with more rapidly evolving nuclear data (e.g., microsatellites, genome-wide SNPs) would facilitate delineation of ESUs. Nonetheless, molecular data from this study paired with available distributional information [34,35] provide ample evidence for the delineation of the Amite River and Mobile drainage populations of P. inflatus as two distinct MUs [102,104,107]. The designation of these MUs ensures the protection of irreplaceable genetic variation, and in particular, emphasizes conservation needs in the highly susceptible Amite River given its limited geographical distribution [34] and presumed extirpation of populations from adjacent systems that were hypothetically closely related based on our comparative phylogeographic approach. Future long-term monitoring efforts will be useful to identify specific population characteristics such as abundance, age-class structure, dispersal capabilities, and reproductive timing within these MUs and may lead to fine-scale delineation of population units.
Captive propagation of freshwater mussels is a critical component of recovery planning for many species [106,108] and likely the only viable option for re-establishing extirpated populations of P. inflatus [35], especially in the Pearl River drainage. Our assessment provides defensible justification for natural resource managers to use individuals from the Amite River rather than the Mobile drainage as a source of brood stock for recovery efforts for P. inflatus that include translocation or captive propagation in the Pearl and Pontchartrain drainages. Based on the likely scenario that extant populations of P. inflatus are restricted to the Amite River and Mobile drainage, possible re-establishment sites include historically occupied reaches in the Bogue Chitto, Comite, Pearl, and Tangipahoa rivers.

Conclusions
Our study revealed congruent patterns of molecular diversification within a group of freshwater mussels with analogous life history traits and dispersal capabilities. Our findings suggest synchronous diversification between Potamilus spp. in the Mobile and PPP, which advances knowledge regarding the drivers of molecular diversification and biogeography of freshwater mussels in these provinces. Patterns of genetic variation in Potamilus spp. recovered by our comparative phylogeographic approach provided defensible justification for the use of the Amite River brood stock for re-establishing P. inflatus in PPP drainages. This finding provides direction for natural resource managers to develop appropriate recovery practices that may include captive propagation and translocation. Although a useful tool, without proper guidance and planning efforts, introduction of captive raised or translocated individuals has the potential to harm existing populations or nontarget species [109,110]. Recovery planning would greatly benefit from robust distributional information for P. inflatus, and future efforts utilizing both eDNA sampling and traditional surveys would help resolve whether the species is truly extirpated from select drainages. We also encourage further evaluations using fine-scale genomic markers and detailed genetic management planning to characterize genetic diversity of brood stock and captively bred individuals in efforts to maximize genetic diversity in augmented or re-establish populations.
Given the imperilment status of freshwater mussel species globally [111], our study represents a useful methodology for hypothesizing the genetic relationships of extant and extirpated populations of imperiled species to facilitate recovery planning. The use of mtDNA may be limited on a regional scale in most species; however, comparative phylogeographic approaches incorporating more rapidly evolving genome-wide markers introduces a more robust methodology for evaluating population dynamics within drainages and even at a local scale using surrogate species. As the understanding of phylogenetic relationships and life history characteristics continues to improve, utilizing comparative phylogeographic methodologies is a promising tool toward effective species recovery and long-term viability of freshwater mussels.

Acknowledgments:
The authors wish to thank Caitlin Beaver for assistance in the laboratory; John Pfeiffer for supplying unpublished data; Matt Cannister for help with SceinceBase; along with Carla Atkinson, Michael Buntin, Matthew Duplessis, Jeff Garner, Paul Johnson, and Kevin Kocot for assistance with collections. Special thanks to Jeff Powell for help obtaining funding, which was provided to Nathan A. Johnson by the U.S. Fish and Wildlife Service and U.S. Geological Survey. Specimens utilized in this study were either from museum collections or collected under the U.S. Fish and Wildlife Service permit TE 697819-4. This work was authored as part of the Contributor's official duties as an Employee of the United States Government and is therefore a work of the United States Government. In accordance with 17 U.S.C. 105, no copyright protection is available for such works under U.S. Law. This is an Open Access article that has been identified as being free of known restrictions under copyright law, including all related and neighboring rights (https://creativecommons.org/publicdomain/mark/1.0/). You can copy, modify, distribute and perform the work, even for commercial purposes, all without asking permission. Any use of trade, firm, or product names is for descriptive purposes only and does not imply endorsement by the U.S. Government.

Conflicts of Interest:
The authors declare no conflict of interest.