Genetic structuring in the Pyramid Elimia , Elimia potosiensis ( Gastropoda , Pleuroceridae ) , with implications for pleurocerid conservation

The Interior Highlands, in southern North America, possesses a distinct fauna with numerous endemic species. Many freshwater taxa from this area exhibit genetic structuring consistent with biogeography, but this notion has not been explored in freshwater snails. Using mitochondrial 16S DNA sequences and ISSRs, we aimed to examine genetic structuring in the Pyramid Elimia, Elimia potosiensis, at various geographic scales. On a broad scale, maximum likelihood and network analyses of 16S data revealed a high diversity of mitotypes lacking biogeographic patterns across the range of E. potosiensis. On smaller geographic scales, ISSRs revealed significant population structure, even over the distance of a few hundred meters. Unlike other freshwater mollusks like mussels, E. potosiensis showed no evolutionary patterns relating to biogeography. The species does show population-level genetic structure, which may have implications in conservation efforts.


Introduction
The Interior Highlands, separated by the Arkansas River Valley into the Ozark and Ouachita regions, is a major geographical feature of southern North America and includes the Ozark and Ouachita Mountains (Fenneman 1928, Mayden 1985).As a result of the geological processes that formed it, the Interior Highlands has a distinct fauna with high instances of endemism and drainage patterns that have dictated diversity (Mayden 1985, Mayden 1988, Matthews and Robison 1998, Austin et al. 2004, Bonett and Chippindale 2004).The Interior Highlands also served as a refugium for many species during pe-With nearly 170 recognized species (Johnson et al. 2013), the family Pleuroceridae is the second largest group of freshwater snails in North America behind only the Hydrobiidae.As part of the 700 snail species in North America, pleurocerids are increasingly imperiled by river regulation, habitat loss, poor water quality, reduced water quantity, and invasive species (Johnson et al. 2013).Only 22% of all freshwater snail species are stable in a conservation context, and recent extinctions support these rankings (Sada and Vinyard 2002, Hershler et al. 2007, Johnson et al. 2013).Thirty-three pleurocerid species are extinct, and of all federally endangered and threatened species, 55% of the listed snails are freshwater species including eight pleurocerids (Johnson et al. 2013).Only six recognized species occur west of the Mississippi River (NatureServe 2017), with three of them occurring in the Interior Highlands.The Pyramid Elimia, Elimia potosiensis (Lea, 1841), is a morphologically variable species found in springs and rivers throughout the Ozark and Ouachita highlands of Arkansas and Missouri, its range stretching westward into eastern Kansas and Oklahoma (Goodrich 1939, Tiemann andCummings 2007).Unlike many other pleurocerids, populations of E. potosiensis are often large and connected to one another within stream flows (Gordon 1980, Gordon 1982), and the species is secure across its range (global heritage rank of G5; NatureServe 2017).Given that freshwater taxa, including mollusks, in the Interior highlands exhibit genetic structure, and that many freshwater snails are of conservation concern, we aimed to answer three research questions.First, on a broad geographic scale, does E. potosiensis represent a single clade showing genetic structuring by population and river drainage?Second, on narrower scales, do populations show any degree of genetic structure within single streams of closely connected waterways?Finally, what impact can our findings have on pleurocerid conservation?
We utilized two forms of genetic data for our study, mitochondrial DNA sequences at the species level and inter-simple sequence repeats (ISSRs) at the population level.The bulk of freshwater gastropod genetic studies employ one of two mitochondrial DNA fragments derived from either the 16S or cytochrome oxidase c subunit I gene.While these sequences may show utility in population and species-level studies, results are mixed in freshwater cerithioideans including pleurocerids (e.g.Köhler and Deein 2010, Miura et al. 2013, Köhler 2016).Widely divergent mitochondrial haplotypes can be found in single populations and species that can cloud questions of monophyly and relatedness (Dillon andFrankis 2004, Whelan andStrong 2016), and no consistent explanation has been offered as to why these divergent haplotypes exist or persist (Whelan and Strong 2016).Not all studies, however, report this issue (Lydeard et al. 1998, Minton and Lydeard 2003, Minton and Savarese 2005).ISSRs are fragments of DNA separating neighboring microsatellites used as genetic markers amplified using polymerase chain reaction (PCR) (Zietkiewicz et al. 1994).Primers for ISSR amplification are complementary to microsatellite sites, bind to them using a one to three nucleotide anchor sequence on the 5' or 3' end (Bussell et al. 2005), and amplify the ISSR in between (Culley and Wolfe 2001).This approach does not require prior genome sequence information, and leads to banding patterns containing multiple loci and high levels of polymorphism (Tsumura et al. 1996).ISSRs have provided useful population data in a variety of organisms including plants (Lisek and Rozpara 2010), insects (Vijayan et al. 2006), arachnids (Machkour-M'Rabet et al. 2009), and birds (Haig et al. 2003), but have not been utilized extensively in gastropods (Dong et al. 2011, Snegin 2014).

Broad scale phylogenetic analyses
We collected live E. potosiensis from six river drainages throughout the Ozark and Ouachita highlands (Figure 1).We followed a standard CTAB/chloroform protocol (Saghai-Maroof et al. 1984) for DNA extraction.We amplified a 500 bp portion of the 16S ribosomal subunit (primers 16sar and 16sbr [Palumbi, 1996]) using Qiagen Taq PCR Core Kit and the following thermal cycling conditions: an initial denaturation cycle of 95°C for three minutes; 40 cycles of 95°C for 35s, 44°C for 45s, 72°C for 45s; and a final five-minute extension period at 72°C.We sequenced gel-purified products on a Beckman CEQ8000 automated sequencer using Sanger dideoxy sequencing.We used Popart (Leigh and Bryant 2015) to calculate nucleotide diversity and Tajima's D (Tajima 1989), and to visualize haplotype diversity in E. potosiensis by generating a TCS (Clement et al. 2002) network.We then combined our E. potosiensis sequences with other pleurocerid 16S data taken from NCBI (Table 1) and aligned them using the Muscle algorithm (Edgar 2004) with default settings.We used Gblocks 0.91b (Castresana 2000) with default settings to remove poorly aligned positions from the data, and used the ModelFinder (Kalyaanamoorthy et al. 2017) function of Iq-tree (Nguyen et al. 2015) to select the appropriate substitution model by Bayesian information criterion.We then analyzed the data under maximum likelihood in Iq-tree and estimated branch support using 10,000 ultra-fast bootstrap replicates (Minh et al. 2013).We also tested the hypothesis that E. potosiensis sequences from each separate river drainages represented separate monophyletic groups.We used the approximately unbiased test (Shimodaira 2002) with 1000 RELL bootstrap replicates implemented in Iq-tree to test whether topologies were significantly different (p<0.05).

Narrow scale population genetic analyses
We employed ISSRs as population markers on two different spatial scales.In our first study, herein referred to as the 'small' study, ten individuals from each of 12 E. potosiensis populations were collected (Figure 2, white circles).Four populations were in a spring run flowing into Walnut Creek in the Arkansas portion of the Ouachita National Forest.Populations were sampled at the mouth of the spring run (S4) and in 50 m intervals upstream (S3-S1); at the mouth, the spring run flowed over a 1 m rock ledge into Walnut Creek.Populations from Walnut Creek were sampled as well, using the confluence point with the spring drainage (C3) as a starting point.Populations were sampled in 50 m intervals upstream (C2 and C1) and downstream (B1-B5, below where the creek flows under the road through a culvert).
We followed the 'small' study with a second 'large' study.We collected 50 E. potosiensis from each of four localities (Figure 2, gray circles).Three sites were located in the Walnut Creek drainage in the Arkansas portion of the Ouachita National Forest.The first (ARK1) was a small spring feeding the main channel of Walnut Creek.The second and third were in Walnut Creek, 325 m upstream and 350 m downstream of where the ARK1 spring run entered the creek.The final site was Sallisaw Creek, Sequoyah County, Oklahoma near Marble City.For both studies, samples were stored in 95% ethanol or frozen at -20°C.
We isolated genomic DNA as before, purified it on silica filters (UltraClean PCR Clean-up Kit, MoBio), and diluted it with sterile water to 50 ng/μl concentration.We amplified ISSRs by PCR in 50 μl volumes with the Go-Taq PCR Core System I reagents (Promega) at the manufacturer's recommended concentrations.The cycling profile consisted of an initial denaturation cycle of 95°C for three minutes, 30 cycles of 95°C for 30s, annealing for 30s, and 72°C for 60s, followed by a ten-minute extension period at 72°C.We used one primer in each ISSR reaction (Table 2) and optimized each annealing temperature.For the 'small' study we visualized ISSRs by running 10 μl of PCR product on 2% agarose gels in TBE at 100 V for 90 minutes.For the 'large' study, 20 μl of PCR product was run on 10% polyacrylamide-TBE gels at 110V for two hours.In both, a 0.1-1 kb size ladder was loaded, and we stained gels 1% ethidium bromide then de-stained in deionized water.Gels were photographed (Bio-Rad Chemi XRS) and processed (Bio-Rad Quantity One) before analysis.We coded individual bands if they were clear and reliably scored in each individual.Bands were identified and sized using GelAnalyzer2010 (http:// www.gelanalyzer.com/),then binned into markers with 5 bp tolerance, For both studies bands were scored as present or absent in each individual to produce separate data matrices.
We analyzed our data matrices using GenAlEx 6.052 (Peakall and Smouse 2006) and the 'poppr' 2.2.1 package (Kamvar et al. 2014, Kamvar et al. 2015) in R 3.3.1 (R Core Team 2016).For each, we generated a genotype accumulation curve (Kamvar et al. 2015) to determine the minimum number of loci necessary to discriminate  between individuals in a population.We calculated Simpson diversity and evenness (Grünwald et al. 2003) at each locus, and Nei's gene diversity (Nei 1973) for each population.An analysis of molecular variance (Excoffier et al. 1992) was performed on the data by populations in regions (spring [S] versus creek upstream [C] versus creek downstream [B], or Arkansas versus Oklahoma) to partition the overall genetic variation.The Φ statistic was used to describe ratios of between and among population and region variance to total variance; Φ values are comparable to F statistics calculated for dominant-only genotype markers (Excoffier et al. 1992).Significance of Φ was determined through standard permutation (10,000 replicates) of the entire dataset.We also used principal coordinates analysis (PCoA) based on Nei's genetic distance to visualize population structure.
Our parsimony network (Figure 3) connected the 30 resulting 16S haplotypes with a maximum of 83 mutations.
Our final data matrix processed through Gblocks consisted of 110 sequences and 317 nucleotide sites.Under maximum likelihood we generated a single tree (log likelihood = -2454.1915,se = 186.2554)under the HKY+R3 model (Figure 4).The E. potosiensis sequences did not form a single clade but rather resolved in four well-supported clades.
Constraining the monophyly of E. potosiensis, however, did not produce a significantly different topology (p = 0.63).Individuals from single localities were not necessarily recovered as sister taxa, and constraining them to be sister produced a significantly less likely topology (p < 0.001).Also, no grouping of sequences by river drainage was observed in the overall tree nor in the four E. potosiensis clades, and constraining the monophyly of each river drainage generated significantly less likely topologies (p < 0.001).

Population genetics
In the 'small' study, we were able to bin and identify 22 markers.There were 110 unique genotypes identified among the 120 total individual snails.A genotype accumulation curve suggested that the minimum number of markers needed to discriminate individuals was 19, so our analyses were close to the detection limit for those markers.
Simpson diversity for individual loci ranged from 0.033 to 0.500, while evenness ranged from 0.383 to 1. Nei's gene diversity by population ranged from 0.248 (B3) to 0.394 (C3).AMOVA (Table 3) showed that most of the variation in the data was within populations (95.86%), followed by between the regions (spring and creek sites) (3.91%) and within each region (0.23%).Permutation tests suggested significant (p < 0.01) genetic structuring between populations and between regions (spring or creek), but not between populations within the same region (p = 0.28).Pairwise comparisons of populations suggested low genetic differentiation between populations (Φ ≤ 0.133), and PCoA showed no clustering by proximity or region (Figure 5).In the 'large' study, we were able to bin and identify a total of 100 markers.Each of the 200 snails possessed its own unique genotype.A genotype accumulation curve suggested that a minimum of 31 polymorphic markers was needed to discriminate individuals, so our dataset of 98 was sufficient for further analysis.Simpson diversity for individual loci ranged from 0.303 to 0.500, while evenness ranged from 0.376 to 1. Nei's gene diversity for the populations ranged from 0.254 for ARK1 to 0.328 for ARK3.An AMOVA (Table 4) showed that most of the variation in the dataset was seen between individuals of a population (66.8% of total).The least variation (12.1%) was seen between populations in a region (Arkansas or Oklahoma), and the remainder (21.1%) existed between regions.Permutation tests indicated significant genetic structuring within populations, between populations, and between regions (Φ = 0.133 to 0.367); pairwise comparisons at all levels were significant at p<0.01.PCoA suggested complete separation of OK from the ARK populations, and little overlap between the three ARK populations (Figure 6).

Discussion
Freshwater taxa from the same region tend to share evolutionary history, resulting in replicated patterns of biogeography and speciation (Walker andAvise 1998, Carini andHughes 2006).This is true in the interior highlands, where the Ozark and Ouachita drainages represent separate areas of endemism (reviewed in Hoagstrom et al. 2014).Phylogenies of freshwater mussels reflect these patterns (Serb and Barnhart 2008, Inoue et al. 2013, Chong et al. 2016), but they remained unexplored in freshwater snails.Our nuclear ISSR data suggest that Elimia potosiensis exhibits genetic population differen-tiation, even at small geographic scales, while our mitochondrial data show no evidence of phylogeography across the species range.Our analyses of mitochondrial 16S sequences suggested no genetic structuring.Mitochondrial haplotypes did not group together in our TCS network, and we recovered E. potosiensis in four well-supported clades instead of a single clade.However, we could not reject the monophyly of E. potosiensis based on our data.No statistical support existed for genetic structure within populations, within river drainages, or within populations or drainages in each of the four clades.We observed no evidence of heteroplasmy to suggest nuclear analogs of mitochon-drial sequences ("numts") in our data, nor did we find any biogeographic groupings of sequences that would suggest cryptic taxa.What we did observe, however, was a pattern where E. potosiensis sequences comprised multiple well-supported divergent clades.As mentioned previously, this is frequently observed in freshwater snails including pleurocerids using mitochondrial markers.Whelan and Strong (2016) showed that the pattern exists across multiple valid species, is not limited to specific taxa, and becomes more apparent as additional sequences from each population are included in an analysis.The maximum number of individuals we sequenced for a single population was five, and two of these individuals exhibited the highest pairwise sequence divergence.We predict that additional E. potosiensis samples would only increase the observed sequence diversity while not increasing phylogenetic resolution.Whelan and Strong (2016) proposed that balancing selection might keep divergent haplotypes in pleurocerid populations.Our Tajima's D of 0.894 suggests balancing selection of 16S haplotypes in E. potosiensis, though population reduction and migration can generate similar results (Maruyama andFuerst 1985, Simonsen et al. 1995).
Four nominal morphospecies comprise the modern notion of E. potosiensis (Burch and Tottenham 1980): E. potosiensis (Lea, 1841), a narrow-range species from Missouri; E. crandalli (Pilsbry, 1890), endemic to Mammoth Springs, Arkansas; E. ozarkensis (Call, 1886), found in springs in north Arkansas; and E. plebejus (Gould, 1850), the most widespread and the shell form currently associated with E. potosiensis.While E. potosiensis does possess shell variation that may have an environmental component (Minton et al. 2011), no additional data supports recognition of any separate morphological entities.Jones and Branson performed detailed examinations of internal anatomy and radula structure of all four nominal E. potosiensis taxa, noting that, without their shells, internal structures of the morphotypes "...can not be separated" (Jones and Branson 1964: 60), and that the shell variations seen among the forms could be found throughout the species' range.They thus treated E. potosiensis as single species exhibiting shell plasticity, where morphological "...variation is the rule rather than the exception" (Jones and Branson 1964: 60).This anatomical data, combined with our inability to reject the monophyly of E. potosiensis, lead us to support the current notion (Burch andTottenham 1980, Johnson et al. 2013) that E. potosiensis is a single widespread species.
While we saw no broad scale genetic patterns in E. potosiensis, our two ISSR studies suggested population-level genetic structuring.Our 'small' study of ten individuals from each of 12 populations represented a pilot study as well as our first effort with ISSRs.In it, we were able to characterize 22 bands visualized on agarose gels using six primers.The AMOVA suggested that most of the genetic variation was within each population, but that a small (Hartl and Clark 1997) yet significant amount of difference existed between the spring and two creek regions.
We predict that the spring is isolated from the creek by the one meter of elevation, as the spring run flows over a rocky ledge and into the creek below.This normally serves as an isolating mechanism between populations in the two channels.Snails within each channel can migrate upstream and downstream freely, and during heavy rains and flooding conditions, the creek can crest above the spring.This could potentially carry individuals from the spring into the creek and vice-versa, and may explain the low pairwise Φ values observed between populations and lack of PCoA clustering (Wilmer et al. 2011).Our 'large' study showed complete differentiation of the four populations examined with high levels of genetic differentiation at all population and region levels (Hartl and Clark 1997).We used large numbers of individuals and markers, both of which increase the resolution and accuracy of population genetic studies (Ruzzante 1998, Kalinowski 2005), and utilized polyacrylamide gels to better visualize bands.Our AMOVA and PCoA suggested that each population is genetically distinct from the others.
Two issues we identified in our 'small' study were the low number of bands resolved by our primers and the small sample sizes from each population.Our band total for six primers was comparable to the number generated by a single primer in other studies, and is likely indicative of two factors.First, choice of our primer sequences was probably not ideal, employing dinucleotide repeats and/ or degenerate 5'-bases ahead of the repeated portion.Results from the literature suggest that trinucleotide repeats and 3'-degenerate anchors after the repeats increase resolution and repeatability (Godwin et al. 1997).Second, we visualized our bands on agarose instead of polyacrylamide.Agarose has far less resolving power, and multiple bands in a small range of sizes might not have been identifiable on our gels (Stift et al. 2003).In regards to sample size, many studies recommend thirty individuals per population (Pruett andWinker 2008, Hale et al. 2012).Our AMOVA suggested treating the regional groupings as populations, where we were able to resolve significant structuring on a small spatial scale given the now larger sample sizes.The differences in our 'small' and 'large' studies highlighted the importance of having the appropriate number of markers and individuals.
Our study also highlighted the need for novel and useful genetic markers for pleurocerid conservation.Nearly 80% of all pleurocerids are imperiled (Johnson et al. 2013), and most species designations still rely on shell morphology.The genetic species and population structure of most pleurocerid taxa remains unknown.Our study provides the latest example where gene sequences are not useful in delineating neither population nor species-level boundaries.Results from analyses of mitochondrial gene sequences are mixed, generating far less resolution (e.g.Minton andSavarese 2005, Minton 2013) than uncertainty (reviewed in Whelan and Strong 2016).Single mitochondrial gene sequences, even from the same population, can differ by over 20%, which can lead to hypotheses of species non-monophyly and aberrant biogeography (Dillon and Robinson 2009).Nuclear gene sequences, on the other hand, are too highly conserved in pleurocerids to be useful (Whelan and Strong 2016).Other forms of nuclear data, however, such as allozymes and now ISSRs, have demonstrated utility in population genetic studies.Allozymes have been used to show genetic connectivity in various pleurocerid genera (Chambers 1980, Dillon 1984, Stein and Stansbery 1984, Dillon 2014), and herein our data suggests ISSRs show promise across large and small geographic levels.No published applications of more modern population genetic techniques (microsatellites, SNPs, etc.) are available, but we hope that they too will be useful in understanding pleurocerids.The continued development of consistently useful tools is needed to help properly understand genetic population and species structure in pleurocerids, especially in a conservation contextS before more taxa are lost.

Figure 1 .
Figure 1.Map of E. potosiensis collection localities, color coded by river drainage, used in the phylogenetic analysis.Detailed locality information is in Table1.

Figure 2 .
Figure 2. Map showing collecting localities for both ISSR studies.Inset shows position of Oklahoma (OK) relative to Arkansas (ARK1-3) populations.Shell of E. potosiensis from population OK is shown.

Figure 3 .
Figure 3. TCS parsimony network of E. potosiensis 16S haplotypes color coded by river drainage.Hash marks represent mutational steps between haplotypes.

Figure 4 .
Figure 4. Maximum-likelihood phylogeny of E. potosiensis sequences.Individual snails are color-coded by river drainage and labelled according to Figure 1 andTable 1. Branches with bootstrap support ≥ 80% are labeled.

Figure 5 .
Figure 5. Principal coordinates analysis of ISSR genetic distances from the 'small' study showing overlap of populations.Populations labeled as in Figure 2.

Figure 6 .
Figure 6.Principal coordinates analysis of ISSR genetic distances from the 'large' study showing separation of populations.Populations labeled as in Figure 2.

Table 1 .
Locality and NCBI accession numbers for sequences used in this study.

Table 2 .
Sequences, annealing temperatures (ºC), and number of bands produced for ISSR primers used in each study. '

Table 4 .
Results of the 'large' study AMOVA.Regions were either Arkansas or Oklahoma.

Table 3 .
Results of the 'small' study AMOVA.Regions were spring (S populations), and creek above (C) and below (B) confluence with the spring.Source d.