Range reduction of Oblong Rocksnail, Leptoxis compacta, shapes riverscape genetic patterns

Sample collection

Leptoxis compacta was collected during two trips to the Cahaba River in June 2018 and June 2019. We collected L. compacta at four sites, and all sites except Cahaba River above Shades Creek were outside the previously documented contemporary range of L. compacta (Fig. 3; Whelan, Johnson & Harris, 2012). At each location, individuals were collected by hand and identified in the field. Despite being a narrow range endemic that has undergone distributional decline, L. compacta was locally abundant where found. Based on qualitative observations, we sampled less than 1% of the population, making our sampling negligible to species survival. Twenty specimens from each site were transported live to the lab, sacrificed following Fukuda, Haga & Tatara (2008), and placed in 96–100% ethanol until tissue clips could be taken. Specimens were collected under an Alabama Department of Conservation and Natural Resources Educational Scientific Collections Permit (License #2019100990068680) or as an agent of the state (P.D. Johnson). Jecca Shumante and Lauren Allred of Living River: A Retreat on the Cahaba provided access to privately owned sites along the Cahaba River. All shells have been cataloged separately to correspond to associated molecular data and deposited at the Auburn University Museum of Natural History (AUMNH 45652-45690; Table 1, Table S1).

Molecular data generation

Tissue clips from 20 individuals per collection site were taken and subjected to a standard proteinase K digestion. DNA was extracted with the Qiagen DNeasy Plant Mini Kit with a minor modification to allow for tissue digested with proteinase K. We used a plant kit because it works well on freshwater gastropods that produce large amounts of mucus polysaccharides (Whelan et al., 2019). The integrity of whole genomic DNA was checked on a 1% agarose gel and quantified with a Qubit DNA assay. Extracted DNA was standardized to 120 ng/µL for 2bRAD library prep.

A reduced representation genomic library was generated for genotyping using the AlfI enzyme and the 2bRAD library prep protocol of Wang et al. (2012). This RAD-seq approach uses a type IIB restriction enzyme that has two recognition sites. AlfI recognizes two sites separated by six base pairs and makes a cut with a one base pair overhang 12 base pairs from the 5′ and 3′ ends of the restriction sites. Following Whelan et al. (2019), we did a 1/16th genomic reduction by using adaptors in the ligation step that had an “NC” overlap, thus only binding to AlfI RAD-loci that had a G base pair at the first base pair of each restriction cut overhang. For more details, see Wang et al. (2012) and the lab protocol on the FigShare repository for this study (DOI: 10.6084/m9.figshare.12014619).

All samples were dual-indexed for pooling. Sequencing occurred in multiple batches. The first batch had 48 L. compacta samples pooled in equimolar concentrations and sequenced on a single lane. The other individuals were pooled in equimolar concentrations with samples from projects on conservation genomics of other pleurocerid species, and 87 individuals were sequenced per HiSeq 4000 lane. Although batch effects in RADseq data have been recently noted in studies that used different read lengths among sequencing runs (Leigh, Lischer & Keller, 2018) and in species introgression studies (Lambert et al., 2019), such issues were not relevant to our sequencing design or study objectives. Nevertheless, we took steps to limit potential batch effects by implementing strict filtering parameters during dataset assembly (see below). Pooled libraries were sequenced on an Illumina HiSeq 4000 with 1 × 75 bp chemistry at University of Oregon Genomics and Cell Characterization Core Facility.

Raw Illumina reads were demultiplexed with the STACKS 1.48 (Catchen et al., 2013) module process_radtags, allowing for one mismatch per barcode. Demultiplexed reads were quality filtered with the script QualFilterFastq.pl (http://github.com/Eli-Meyer/sequence_processing) for any read that had five or more base pairs with Phred quality scores less than 20. Reads were processed with scripts from SHRiMP 2.23 (Rumble et al., 2009) and subsequently trimmed to AlfI RAD-loci with the script AlfIExtract.pl (http://github.com/Eli-Meyer/2bRAD_utilities). As this step removes any sequence that is not part of the RAD-locus, adaptor sequences and non-target sequences are removed from the sequencing reads. RAD-loci, defined as the stretch of DNA cut by the AlfI enzyme, were assembled with the STACKS 1.48 pipeline denovo_map.pl as no reference genome is available for L. compacta. For denovo_map.pl parameters, we set minimum stack depth to five (-m 5), distance allowed between stacks to three (-M 3), and distance between catalog RAD-loci to two (-n 2). These parameters were determined to be most appropriate for our data following Paris, Stevens & Catchen (2017). All other denovo_map.pl parameters were set to defaults.

After assembly, RAD-loci were filtered for missing data using the STACKS program populations. In order to pass filtering steps, a RAD-locus had to be present in 75% of individuals from any given collection site and also present at three collection sites. RAD-loci that had a minimum minor allele frequency of less than 2.5% or heterozygosity higher than 50% were removed to limit the influence of paralogy and misassembly on final datasets. Sequencing coverage of RAD-loci with SNPs was measured with vcftools (Jombart & Ahmed, 2011). Kinship coefficients among individuals were inferred with KING (Manichaikul et al., 2010). Files output by STACKS were formatted for KING with PLINK 1.9 (Chang et al., 2015), and pairwise kinship coefficients were calculated with the KING flag “–kinship”. Given the absence of closely related individuals, no further dataset filtering was done.

After filtering, a dataset that included all SNPs per RAD-locus and a dataset with only one random SNP per RAD-locus were generated. We assume that RAD-loci are unlinked and that the one SNP per RAD-locus dataset had zero linkage disequilibrium. Analyses employed the one SNP per RAD-locus dataset, unless otherwise noted.

Population genetics analyses

Average observed heterozygosity (H_o), expected heterozygosity (H_e), nucleotide diversity (Π), and F_IS at each collection site were calculated by populations. The number of private alleles at each site was also reported by populations. Average allelic richness (A_r) was calculated with the R (R Core Team, 2020) package diveRsity (Keenan et al., 2013). An analysis of molecular variance (AMOVA; Excoffier, Smouse & Quattro, 1992) was done with the R packages poppr (Kamvar, Tabima & Grünwald, 2014) to test genetic structure among collection sites. AMOVA was implemented with the function “poppr.amova” using the ade4 method (Dray & Dufour, 2007) and 10,000 permutations.

We also tested for a pattern of isolation by distance by measuring the correlation between pairwise F_ST values and geographical distance between collection sites. Pairwise F_ST values were calculated using the Weir & Cockerham (1984) method with the R package hierfstat (Goudet, 2005). Stream distance was measured in Google Earth by tracing paths between collection sites along the Cahaba River (see Table 2). River distance was used, rather than straight line distances, because migration over land is impossible for gill breathing pleurocerids. A Mantel test was performed with the R package ade4 (Dray & Dufour, 2007), and significance was tested with 1,000 permutations. However, Mantel tests have been criticized as a method for testing isolation by distance (Legendre, Fortin & Borcard, 2015; Meirmans, 2015), so we also performed a multiple regression on distance matrices with 1,000 permutations using the R package ecodist and its MRM function (Goslee & Urban, 2007). In addition to a pattern of isolation by distance, past studies have shown that many freshwater organisms, including pleurocerids, display a pattern of increased genetic diversity in more downstream populations (Paz-Vinas et al., 2015; Whelan et al., 2019). Therefore, to better assess riverscape genetic patterns of L. compacta, we performed linear regression of distance from the most downstream site against H_o,H_e,A_r, and Π. Linear regressions were done in R.

We examined clustering of L. compacta genetic data with discriminant analysis of principal components (DAPC). We used the multiple SNPs per RAD-locus dataset and the R package adegenet (Jombart & Ahmed, 2011) to perform DAPC. We first used the adegenet function “find.clusters” testing up to 25 clusters and using Bayesian information criteria (BIC) to identify the best-fit number of clusters for our data. Using the number of clusters with the lowest BIC value, we performed a DAPC with the adegenet function “dapc” and plotted the results in R.

We inferred genomic admixture of L. compacta individuals with ADMIXTURE 1.3 (Shringarpure et al., 2016). ADMIXTURE assumes zero linkage disequilibrium, so we used the one SNP per RAD-locus dataset. ADMIXTURE analyses were run with the AdmixPipe pipeline (Mussman et al., 2020). To determine the best-fit number of clusters (K) for our data, K values from 1 to 5 were assessed with 20% cross-validation. Twenty replicates of ADMIXTURE were run at each K, and the best-fit K was determined as the value that had the lowest average CV score across replicates. ADMIXTURE results were visualized with Clumpak (Kopelman et al., 2015).

Genomic co-ancestry among individuals was also assessed with fineRADstructure (Malinsky et al., 2018). Unlike ADMIXTURE, fineRADstructure can use linked SNPs and provides additional information on individual genomic background. Thus, the multiple SNPs per RAD-locus dataset was used for fineRADstructure analyses. First, a co-ancestry matrix was inferred with the script RADpainter. Subsequently, clustering was done with the Markov chain Monte Carlo method of fineRADstructure, running for 500,000 generations and sampling every 1,000 generations; the first 200,000 generations were discarded as burn-in (non-default parameters: -x 200000 -y 300000 -z 1000). We also inferred a tree for visualization with fineRADstructure using the tree-building algorithm of Lawson et al. (2012) with 10,000 attempts (non-default parameters: -m T -x 10000). fineRADstructure results were plotted with R scripts included in the fineRADstructure package.

Code and data availability

All bash and R scripts used for processing and analyzing data are available at github.com/nathanwhelan. Demultiplexed raw Illumina reads have been uploaded to NCBI under BioProject PRJNA631794. Assembled datasets in various file formats (e.g., vcf, genepop) and the 2bRAD library prep protocol are available on FigShare (DOI: 10.6084/m9.figshare.12014619).

Sample collection

During survey work, we collected L. compacta from Cahaba River at old Marvel slab upstream to Cahaba River at Booth’s Ford (Fig. 3). All sites except Cahaba River at Shades Creek are sites where L. compacta was not found during survey work over the last 30 years. Our collections represents a 1.83 km downstream range extension and a 4.76 km upstream extension compared to the previously documented contemporary range of L. compacta (Whelan, Johnson & Harris, 2012). While this study was ongoing, 3 putative L. compacta individuals were collected at Cahaba River at Belle Ellen Shoals (Fig. 3) during a general mollusk survey (Johnson, 2019). However, species identification was uncertain and L. compacta appeared exceedingly rare. Therefore, individuals from Cahaba River at Belle Ellen Shoals were not included in our analyses, and we consider this record unconfirmed without additional positive survey results.

Molecular data and population genetics

DNA yields for two individuals were too low for library preparation so only 19 individuals were sequenced from Cahaba River at old Marvel Slab and Cahaba River at Booth’s Ford. The number of demultiplexed raw reads per individuals varied from 930,062-10,146,649 (mean = 4,836,812). Much of the variation in raw reads can be attributed to whether the individual was sequenced on a HiSeq 4000 lane with 48 or 87 samples. Aside from raw read number, we saw no evidence of batch affects like individuals from one sequencing run all clustering together in analyses (see below). After initial raw read filtering, the number of reads that passed quality filtering steps ranged from 865,314–9,838,187 (mean = 4,632,510). Assembly with the STACKS denovo_map pipeline resulted in 105,542 RAD-loci. Filtering with populations, including removal of 4,009 invariant RAD-loci that passed all filters, resulted in a dataset with 4,962 RAD-loci with at least one SNP. Per individual average sequencing coverage of filtered RAD-loci with at least one SNP, excluding missing genotypes, ranged from 31.7–343.2. Average sequencing coverage across variable RAD-loci, excluding missing genotypes, was 163.7. Kinship coefficients inferred with KING indicated that no individuals were closely related (i.e., half or full siblings).

The number of private alleles at each site ranged from 28–262 (Table 1). H_o at each collection site ranged from 0.0963–0.1568, and H_e ranged from 0.0980 to 0.1801 (Table 1). At each site, H_o was lower than H_e, except at Cahaba River at canoe launch where H_o was 0.001 greater than H_e (Table 1). The difference between H_o and H_e was largest at Cahaba River above Shade Creek and Cahaba River at Booth’s Ford. A_r and Π ranged from 1.4511–1.8241 and 0.1010–0.1829, respectively (Table 1). F_IS values ranged from 0.0134–0.1934 (Table 1), with the highest values being at Cahaba River above Shades Creek and Cahaba River at Booth’s Ford. Overall, genetic diversity was greatest at the most upstream site, Cahaba River at Booth’s Ford, and lowest at the most downstream site, Cahaba River at old Marvel slab. All linear regressions of diversity statistics vs distance from the most downstream site were non-significant (p ≥ 0.169).

Pairwise F_ST values among sites ranged from 0.0–0.055 (Table 2). We found no evidence of an isolation by distance pattern among sites (Mantel test, p = 0.843; multiple regression, p = 0.428). According to the AMOVA, significant genetic structure was present among collection sites (p = 0.004), but only 4.16% of variation was explained by collection site. In contrast, 81.8% of genetic variation was explained by within individual variation, further indicating high amounts of gene flow among collection sites.

DAPC indicated two genetic clusters were present in our data. Data were explained by a single discriminate function, and results are therefore presented as a frequency histogram (Fig. 4). ADMIXTURE analyses indicated that genetic diversity from two ancestral populations were present in our data (K = 2). Most individuals from across the range of L. compacta had a genomic admixture profile that was dominated by a genomic background from a single ancestral population, likely indicating that overall genomic diversity has been lost across the range of L. compacta. Nevertheless, 14 individuals had varying levels of admixture with a second ancestral population (Fig. 3). fineRADstructure analyses corroborated ADMIXTURE analyses as two semi-distinct groupings were recovered by fineRADstructure (Fig. 5, Fig. S1). fineRADstructure groupings did not correspond to collection site or any other obvious variable, indicating gene flow among collection sites. Notably, six individuals with comparably high co-ancestry proportions (upper right of co-ancestry matrix in Fig. 5) correspond to individuals in ADMIXTURE analyses with a large proportion of genetic background from the less common ancestral population (represented by orange in Fig. 3).

Genetic diversity across a small landscape

The two most distant collection sites in this study were separated by a smaller distance (9.2 km) than all but two sites sampled for L. ampla in a previous study (Whelan et al., 2019). Therefore, it is difficult to make direct comparisons between genetic patterns of L. ampla and L. compacta. However, we can leverage differences in geographical scale between the two studies to make inferences about fine-scale versus long-distance genetic patterns in pleurocerids. F_ST values among L. compacta collection sites (Table 2) were much lower than values determined for populations of L. ampla (F_ST 0.377–0.773; Whelan et al., 2019). Furthermore, even though AMOVA indicated significant genetic structure among L. compacta collection sites, the small amount of genetic variation that is explained by collection site probably limits its biological relevance. Overall, these data indicate that pleurocerid riverscape genetic patterns across small distances will not always follow common patterns such as isolation by distance and increased genetic diversity at more downstream collection sites. This is likely attributable to gene flow and random drift that prevent the establishment of genetic patterns typically seen across more geographically separated collection sites. From a historical standpoint, we hypothesize that L. compacta previously displayed an isolation by distance pattern across its range, similar to the patterns determined for L. ampla (Whelan et al., 2019). We think this scenario is likely given limited dispersal abilities of pleurocerids and patterns established for L. ampla, a species that retains a much larger portion of its historic range in the Cahaba River drainage than L. compacta. Whether or not there was a similar historical pattern of increased genetic diversity in downstream populations of L. compacta is more difficult to infer, as such a pattern may not be influenced solely by dispersal ability.

Given the well-documented decline of L. compacta, a small number of individuals with a less common genomic background suggests that the species has lost genetic diversity through bottleneck and drift. Patterns seen in DAPC, ADMIXTURE, and fineRADstructure were not driven by geography as individuals with the less common genomic background were not found in adjacent sites (orange in DAPC and ADMIXTURE plots and upper right corner of fineRADstructure plot; Figs. 3–5, Fig. S1). Although individuals with some admixture from the uncommon ancestral population may be present in unsampled individuals at Cahaba River at old Marvel slab and Cahaba River at canoe launch, they would be uncommon. Recent migration is an unlikely explanation of observed co-ancestry profiles as it would indicate that a sizeable population of L. compacta would have to exist elsewhere in the Cahaba River. The most likely hypothesis for explaining observed clustering and co-ancestry profiles (Figs. 3–5) is a genetic bottleneck resulting from species decline in the 20th century. In this scenario, L. compacta was genetically diverse across its historical range prior to decline, but range contraction caused a considerable loss of genetic diversity. In turn, genetic drift resulted in the observed coancestry pattern of one ancestral population being more common in extant individuals (Figs. 3 and 5).

Broadly, genetic structure across the current range of L. compacta can be characterized by a single population with some subpopulation structure at Cahaba River above Shades Creek and Cahaba River at Booth’s Ford (Figs. 3–5; Fig. S1). The subpopulation structure appears to be causing a Wahlund effect (Wahlund, 1928). That is, the Wahlund effect predicts the lower H_o values compared to H_e values and the higher F_IS values seen in collection sites with inferred subpopulation structure (Fig. 3; Table 2). An alternative explanation for the observed pattern of F_IS and H_e is null alleles. However, null alleles are unlikely as they would increase pairwise F_ST values (De Meeûs, 2018) that are uniformly low across populations (Table 2). Despite the putative presence of a Wahlund effect, Cahaba River above Shades Creek and Cahaba River at Booth’s Ford have greater genomic diversity than the two other sites (Table 1; Figs. 3–5). These sites may have better habitat suitability than the other two, allowing for L. compacta to persist with greater genetic diversity as the species declined in the 20th century.

Conservation of Leptoxis compacta

Leptoxis compacta suffered a massive decline during the 20th century, a period of intense mining, forestry, and urban development in the Cahaba River drainage (Onorato, Angus & Marion, 2000; Pitt, 2000; Shepard et al., 1994; Tolley-Jordan, Huryn & Bogan, 2015). The decline was so drastic that L. compacta was considered extinct less than a decade ago. Fortunately, L. compacta is not extinct. However, conservation efforts are needed to ensure the long-term survival of L. compacta as the species is at risk from both chronic habitat degradation and one-time catastrophic events. Two potential management strategies for L. compacta are habitat restoration and reintroduction with captively reared individuals.

In this study, we report an 8.26 km known range extension for L. compacta. One site, Cahaba River at old Marvel Slab, was previously the focus of intense habitat restoration through the removal of a low-level dam (Johnson et al., 2013). The site may have also benefited from improved water quality in Shades Creek (ADEM, 2007; ADEM, 2012) as the site is just below its confluence with the Cahaba River. Since removal of the low-level dam, increases in fish abundance and diversity have been reported (Bennett et al., 2015). Considering L. compacta was not found at this site by Whelan, Johnson & Harris (2012), we think habitat either improved from a point where L. compacta could not survive or from a point of considerably lower carrying capacity. As the only undammed, major river in the southeastern United States, the Cahaba River is much less modified than most other systems in the southeast. Our findings suggest that imperiled gastropods will benefit from water quality and habitat improvements even in relatively “pristine” river systems. Improving habitat, or identifying suitable habitat, will be a necessary starting point for L. compacta reintroduction efforts.

In addition to having a small range, L. compacta only exists along a single river path. This means that one catastrophic event such as a massive point source pollution event above Cahaba River at Booth’s Ford could result in extinction of L. compacta. Such an event is not merely a hypothetical. In 2016, a gasoline pipeline spill came perilously close to the Cahaba River (Pillion, 2016). To mitigate the risks of a single catastrophic event, reintroduction efforts should emphasize range expansion outside the mainstem Cahaba River. Of course, reintroduction efforts also must be limited by the historical range of any given species. Thus, lower Buck Creek is potentially an ideal reintroduction site if habitat quality is sufficient for the persistence of L. compacta. Once a suitable reintroduction site is chosen, managers will need to choose a broodstock site. This decision should be informed with genetic data. The absence of an isolation by distance effect across the current range of L. compacta indicates that managers do not need to prioritize potential broodstock sites based on whether they are geographically proximate to reintroduction sites. Rather, sites with high genetic diversity and ease of access should be prioritized for broodstock. Therefore, the Cahaba River above Shades Creek is likely an ideal broodstock location. Moreover, L. compacta is easy to sample and relatively easy to distinguish from other sympatric species at Cahaba River above Shades Creek, making it ideal from both a genetic and sampling standpoint.