Ancient DNA from shells reveals delayed genomic erosion and rapid immune adaptation in the critically endangered black abalone

Authors: T. Brock Wooldridge, Joshua D. Kapp, Sarah M. Ford, William E. Seligmann, Holland C. Conwell, Talia Tzadikario, Jonas Oppenheimer, Zachary G. Anderson, Alan Le Moan, Alicia Abadía-Cardoso, Peter Raimondi, Beth Shapiro

Abstract

Predicting the genetic consequences of population decline is a major problem in conservation genomics. The black abalone was a culturally and economically important North American mollusk that declined by ~99% in the 1980s due to a rapid disease outbreak. To investigate the consequences of this event, we sequenced ancient genomes from shells spanning the past 1,500 y, producing the first population genomic dataset derived from this data type. This time-series data show that genomic impacts from this near-extinction event are extremely limited and that natural selection has acted recently on genes related to immunity. Together, these data are encouraging for the genomic future of this species and provide practical insight for how to manage recovering populations.

Rapidly declining populations are susceptible to genomic erosion, or the combined effects of reduced diversity, inbreeding, and mounting genetic load that can increase extinction risk (1). However, the timing and magnitude of these effects is far from predictable, as illustrated by the recent history of the black abalone Haliotis cracherodii. This intertidal mollusk was once abundant along the west coast of North America, serving as a common food item and cultural keystone species (2). In 1985, a bacterial disease known as “Withering Syndrome” (WS) first appeared and drove an estimated ~99% decline within just a few years (3, 4). Despite this well-documented collapse, the genomic consequences of this rapid decline are not clear (5–7). Black abalone across California have a high effective population size [Ne ~ 300,000 (8)], show extraordinarily high contemporary genetic diversity, and exhibit no evidence for genetic isolation between sites (7). All of these features appear inconsistent with expectations for a population that recently underwent near-extinction event (9–11). Without a historical baseline it is difficult to ascertain if any current genomic patterns have been affected by the recent bottleneck, which can impact how management strategies proceed.

It is common to see a mismatch between census population declines and population genomic data. Both empirical and theoretical work shows that many of the signatures of genomic erosion may be slow to reflect population decline (12, 13). This delay is often referred to as a “time lag” (14) or “extinction debt” (15). Compounding this, reconstructing recent changes based solely on modern genomes is challenging (16). Temporal genomics, or the time series analysis of genomes (17, 18), offers a powerful way to directly measure genomic erosion over time (19, 20). Temporal genomic studies of bottlenecks have revealed both expected patterns (21–24) and unexpected patterns, including limited genomic erosion (13, 17, 25, 26) or increases in diversity (27). Emerging data and theory suggest that time lags are linked to life history and long-term demographic processes (28). The life history traits of black abalone make a time lag probable—they can live up to 30 y, reach sexual maturity at 4 y, and have highly overlapping generations (28–30). However, a predictive framework that could estimate the time lag effect remains elusive. Without this framework, prebottleneck genomes provide the best way to compare current populations to a precollapse baseline and forecast future change.

Temporal genomics also provides the opportunity to track allele frequency changes through time, which may indicate natural selection. Temporal genomics studies have revealed genes underpinning adaptation to rapid selection events in wild populations, including cold tolerance in Anolis lizards over a single generation (31), and viral resistance in rabbits over tens of generations (32). In other cases, signals of adaptation have been difficult to detect, which might be attributed to subtle signals of polygenic adaptation or the lack of heritable adaptive variation (17). In black abalone there is some evidence for heritable differences in Withering Syndrome susceptibility (33), which is further supported by complementary experiments in other abalone species (34, 35). If WS resistance has indeed evolved, then it is possible that pre- and postbottleneck genomic comparisons could reveal the gene(s) underlying WS adaptation. Previous work did not identify any obvious targets of selection aside from a chromosomal inversion associated with latitude (7). This inversion remains a compelling candidate as latitude is associated with both temperature and WS spread (3, 36), and chromosomal inversions are frequently implicated in local adaptation (37). The identification of any loci showing local adaptation over space or time will be key in guiding the translocation plans that are only just beginning (38).

Genomes from prebottleneck black abalone are needed to measure genomic erosion and selection, but generating genomic data from ancient and historic mollusk specimens is challenging. DNA damage can accumulate quickly through environmental and chemical mechanisms if preservation (e.g., freezing) is insufficient to slow down this process (39). Because shells and dry preparations account for 91% of all abalone in Malacology (mollusk) collections (40), specimens available for genome sequencing are limited to sample types with inherently damaged DNA. Specialized sample processing and analysis methods are needed to handle the small quantity of degraded DNA present in these specimens. Thus far, less than 100 mollusk shells across four studies (41–44) have been processed for whole genome DNA analysis (i.e., shotgun sequencing), none of which have generated multifold nuclear genome coverage. In contrast, for vertebrate species it is common to generate whole genomes for hundreds of ancient individuals in a single study (45, 46). In mollusks, DNA is entrapped during the shell formation process and higher organic content in the shell structure underlies greater DNA preservation (47, 48). Consistent with this, endogenous DNA content has been reported to be as high as 30% for some specimens (44). Given this evidence and a growing understanding of how different shell features preserve DNA (44), it seems feasible to conduct population genomic studies based on large sets of shells.

Here, we successfully apply best practices in ancient DNA sequencing to a time series of shells spanning 1,500 y of black abalone history. Using multifold coverage shell genomes, we are able to directly assess changes in genetic diversity, inbreeding, population structure, and natural selection over time, particularly in relation to the Withering Syndrome bottleneck. We then use these findings to predict scenarios of future genomic erosion and identify variation that could be adaptive. These findings are key to managing the recovery of this critically endangered species.

Results

Temporal Genomics of Black Abalone.

We sequenced whole genomes from 59 black abalone shells spanning the species’ range and the past 1,500 y (Fig. 1). Most shell specimens are dated from 1914–1979, prior to the first mass mortalities due to Withering Syndrome in 1986 (4). We also included one shell from the present day to serve as a control sample. To obtain these data, we took ancient DNA protocols optimized for bone and applied them to abalone shell fragments (49, 50). Specifically, we pulverized a ~50 mg piece of each shell, pretreated the resulting powder with bleach (44), extracted DNA from the powder using a silica spin column for small DNA fragment recovery (50), and generated single-stranded DNA libraries from each extraction (49). All steps were performed in a dedicated ancient DNA lab to minimize contamination. Sequencing libraries resulting from this pipeline ranged in endogenous DNA content from 0.5 to 68.6%, although the median content among 20th century museum shells was 49.5%. DNA damage also scaled with sample age. For the samples dated 1500 BP, cytosine deamination frequency at the 5’ termini and average mapped fragment length were 21% and 60 bp, respectively. These values were 4.2% and 104 bp in the 20th century shells, indicating far lower rates of damage (SI Appendix, Fig. S1). Our final shell dataset consisted of 44 samples sequenced at ~2× coverage and 15 sequenced at 20× or greater coverage. This latter group included a 34× genome from a shell midden dated to 1500 BP.

Finally, we combined these ancient and historic shell genomes with 138 modern black abalone genomes from previous work (8) and 16 modern genomes from the Baja California range, new to this study.

No Change in Genetic Diversity, Inbreeding, or Load Following the Withering Syndrome Bottleneck.

A direct comparison of historic and modern black abalone shows no decrease in genetic diversity over time (Fig. 2). After accounting for DNA damage, we found that individual heterozygosity changed little, remaining high at 1.2% across all time periods (Fig. 2A). A slight postbottleneck increase is noticeable with some particularly high diversity individuals, although this effect is dampened in a complementary analysis of high-coverage, transversion-only polymorphisms (SI Appendix, Fig. S2; see Methods for comment on reference bias). The fraction of the genome in runs of homozygosity (FROH~) remained below 1% across all time periods and was 0% for the majority of individuals, although the “Commercial Fisheries” periods exhibited substantially more variation than all other periods (Fig. 2B). Finally, we also recorded limited change in either masked (Fig. 2C) or realized (Fig. 2D) genetic load across time periods (51). Overall, we see no evidence of lost diversity or increased inbreeding depression following the Withering Syndrome bottleneck.

Subtle Changes in Population Structure Occur Following the Withering Syndrome Bottleneck.

We observe only subtle changes in population structure between pre- and postbottleneck black abalone. An initial genetic PCA showed that all samples fell into one of three discrete clusters with some notable outliers (Fig. 3A). Upon closer inspection, we found that all outliers belonged to one of the following 1) shells from a ~1500 BP site at the northern end of the range (blue circle, Fig. 3A), 2) historic and modern samples from Isla Guadalupe where a putative subspecies has been described (53) (red circle, Fig. 3A), 3) historic and modern samples from Faro San Jose, representing the southern end of the species’ continental range (purple circle, Fig. 3A). Samples from this last group were outliers along PC2, but along PC1 aligned with the three primary clusters.

The three PC1 clusters recapitulate the structure described in Wooldridge et al. (7) that was attributed to a 31 Mb chromosomal inversion on scaffold 4 (Fig. 3B). After removing SNPs within the boundaries of the inversion locus, the three clusters disappear and only a diffuse cloud of samples remains (Fig. 3C). Within this cloud we see separation of historic from modern samples along PC2, and this separation is not driven by sequencing coverage. Consistent with this, median Hudson’s FST between prebottleneck (1914–1979) and postbottleneck samples is only 0.005, and the ratio of median pairwise genetic diversity (π) between the two timepoints is only 1.06 (SI Appendix, Fig. S3; transversion sites and high coverage samples only). In light of our consistent spatial sampling over time, these results indicate very limited genetic divergence between pre- and postbottleneck black abalone.

A Time-Lag Effect on Genomic Erosion.

Simulations show that the stability in diversity, inbreeding, and structure we observe is possible even following a severe bottleneck (Fig. 4). To determine this, we first summarized the average change in heterozygosity, FROH and FST between pre-bottleneck, “Commercial Fisheries” era samples (1914–1979) and post-bottleneck samples. We then simulated populations undergoing a bottleneck followed by migration between newly formed subpopulations (SI Appendix, Fig. S4). The declines and recoveries we simulated were informed by records of Withering Syndrome impacts and recent ecological surveys (SI Appendix, Fig. S5) (4). We see that our time-series measurements are close to simulated ones at 10 generations postbottleneck, which is roughly the number of generations since WS first appeared (Fig. 4). This is true even if we simulate a population reduction of 99.9%, which is more extreme than most decline estimates (4) (Fig. 4B). Our simulations also show limited divergence between the newly formed subpopulations that arise postbottleneck. These results indicate that a time-lag may explain why there is little evidence for genomic erosion in the present day.

However, as postbottleneck time increases in our simulations, we do begin seeing clear signs of genomic erosion. This is most evident in the more extreme 99.9% bottleneck simulation (Fig. 4B). Only recovery rates of 1% or greater are able to stabilize heterozygosity, FROH, and both measures of FST in the long term (SI Appendix, Fig. S6B). In contrast, the 99% bottleneck does not create significant short term (Fig. 4A) or long term (SI Appendix, Fig. S6A) genomic erosion even if population recovery is 0%. These projections highlight how much the initial severity of the bottleneck determines outcomes at short and long timescales. At 99.9% intensity, delayed genomic erosion is likely unless recovery is significant and sustained. A lesser 99% bottleneck, which reflects the range-wide summary of WS decline (4), is less likely to result in genomic erosion even if recovery is negligible.

Evidence for Recent Selection at Immune Loci.

We identified large genomic regions exhibiting postbottleneck balancing selection (Fig. 5). For our selection analyses, we removed the oldest samples, focusing just on prebottleneck (1914–1979) and postbottleneck individuals, then calculated the difference in π between both time periods (“πpost-pre”) along with Hudson’s FST in 50 kb and 10 kb sliding windows. Windows in the top 1% of FST and πpost-pre correspond to genetic divergence between the time periods accompanied by an increase in genetic diversity toward the present day. We tentatively label these windows as under balancing selection (54, 55), although we acknowledge that other selective processes could be responsible (56). We identify outlier windows with these features clustered in large islands on scaffolds 7, 12, 14, and 19. Overlapping regions between the 50 kb and 10 kb outlier windows contain 42 genes, including a large array of innate immunity and gamete recognition genes (Dataset S1). One of these, the scavenger receptor gene DMBT1, is known to be involved in immune responses to bacterial challenges in abalone (57). A phylogeny of the DMBT1 outlier region confirms that prebottleneck diversity is a small subset of modern diversity at this locus (Fig. 5E)

Windows in the top 1% of FST and the bottom 1% of πpost-pre may indicate a recent selective sweep, as these genomic regions have both diverged and decreased in relative diversity towards the present. Of the only six genes consistently recovered from selective sweep windows, one gene—a homolog to SVEP1 (Sushi, von Willebrand factor type A, EGF and pentraxin domain-containing protein 1)—is a compelling candidate. SVEP1 has been associated with total body weight in a QTL study of South African abalone (58), and also acts as a common shell matrix protein and has immune function in mollusks (59, 60). A phylogeny of the outlier region overlapping SVEP1 shows three distinct clades of short branches, suggesting multiple selected haplotypes and a soft sweep (Fig. 5D). Although not an FST outlier, we also see a dramatic drop in diversity at the 31 Mb inversion on scaffold 4. Both the inverted and reference alleles appear to have lost diversity after the bottleneck (SI Appendix, Fig. S7).

In sum, balancing selection in association with the Withering Syndrome bottleneck seems to have been more common than positive selection, although classic selective sweeps occurring on this timescale are difficult to detect (Discussion and Limitations of this study). Nevertheless, these analyses uncovered several genes that are plausibly involved in black abalone immune adaptation over recent generations. We have included a list of these genes and the relevant literature in Dataset S1.

Parallel Increases of Two Inversion Alleles Following the Withering Syndrome Bottleneck.

Two inversions on separate chromosomes showed similar frequency shifts following the WS bottleneck (Fig. 6). Initial genotyping of the two inversions, one on scaffold 4 and the other scaffold 9 (SI Appendix, Fig. S9), suggested that both were more common at higher latitudes. To formally test this, we fitted stable, linear, and sigmoid clines to inversion presence or absence across latitude. Model comparisons showed that sigmoid clines were the best fit for each inversion in both pre- and postbottleneck time periods, although stable clines also had high support in the case of the scaffold 9 inversion (SI Appendix, Table S1). The inflection point of each cline remained at or just north of Pt. Conception, a common biogeographic barrier in marine systems (61) (Fig. 6 A and B and SI Appendix, Table S1). While a clinal pattern is always present, south of Pt. Conception both inversions doubled in frequency following the bottleneck (Fig. 6 A and B and SI Appendix, Fig. S8). Range-wide postbottleneck increases of the scaffold 4 and scaffold 9 inversions amount to 23.9% and 17.6%; these values fall in the top 1.5% of genome-wide allele frequency shifts. In addition to overall changes in frequency, polymorphism within both the inverted and collinear alleles drops following the bottleneck (SI Appendix, Fig. S7).

The two inversion loci are in linkage disequilibrium with each other south of Pt. Conception. An initial analysis showed a significant correlation between having a copy of the scaffold 4 inversion (genotype = A_) and a copy of the scaffold 9 inversion (genotype = B_) south of Pt. Conception (r^2^ = 0.185; P = 1.66e−05). This correlation was absent at northern sites (r^2^ = 0.015; P = 0.22). Explicit tests of two-locus Hardy–Weinberg Equilibrium (HWE) show that, in postbottleneck samples only, there is a significant difference between observed and expected inversion genotypes south of Pt, Conception (X^2^ test, df = 6, P = 7.79E−04; SI Appendix, Table S2). Surprisingly, this deviation appears to be driven by an excess of individuals without either inversion copy (genotype = aabb). We observe no appreciable deviation from HWE for any other combinations of time point and geography (P > 0.048).

Discussion

Mounting evidence shows that time lags are common and can persist for long periods following a demographic bottleneck (28). For black abalone, we see that the nearly 40 y following the Withering Syndrome bottleneck has been insufficient to change genetic diversity, inbreeding, or structure in a meaningful way. We hypothesize that, regardless of the original bottleneck intensity or any ongoing recovery, the “time-lag” effect means that it is too early to see diversity loss, inbreeding depression, or population fragmentation take place. Fortunately, the perspective afforded by ancient DNA also suggests that black abalone may avoid future genomic erosion and that locally adaptive variation is being maintained.

Shells as a Genomics Lens Into the Past.

Our work represents the first population genomic dataset from mollusk shells and demonstrates that it is possible to consistently and affordably generate whole genome data from this sample type. We attribute this success to the application of best practices in ancient DNA analysis. All work from shell powdering to library preparation was performed in a dedicated clean facility to reduce sources of contamination (62, 63). Additionally, bleach treating the shell powder (44, 64), performing extractions with silica spin columns optimized for small DNA fragment recovery (50), and using a single-stranded DNA (ssDNA) library preparation approach (49) produced complex libraries from low input DNA extracts. The combination of these elements allowed us to consistently generate multifold coverage of historic and ancient abalone shells, including a 34× genome from a 1,500-y-old shell. Other recent analyses of shells, including experimentation with shell morphology and treatment methods (44) and the successful capture of nuclear loci from a 100,000 Ka mussel (41), confirm that shells are underutilized reservoirs of DNA. In future work, it would be useful to explore shell DNA preservation in relation to the specimen’s age, deposition environment, and shell mineral composition. It is also important to recognize that practices less optimized for ancient DNA may already be sufficient to recover mitochondrial genomes or target capture loci from shells, providing the data needed to answer key conservation questions (26).

Demography and Selection in Black Abalone Following Withering Syndrome.

A time lag best explains why genetic diversity, inbreeding depression, and population structure have remained stable in black abalone following the Withering Syndrome bottleneck (Figs. 2 and 3) (14). Even centuries-long delays to changes in heterozygosity, inbreeding, and genetic load are supported by theory (12) and have been documented in other threatened populations (13). Time lags are more likely to occur in black abalone because they mature late, have long life spans, and highly overlapping generations, all of which extend the influence of prebottleneck individuals over time (28). Because we were unable to incorporate these life traits into our simulations, the real time lag is probably underestimated (Fig. 4; see Limitations of this Study). Also contributing to time lag is the large prebottleneck Ne and range of black abalone. Formerly large and diverse populations require more generations for heterozygosity to be impacted (65), and the estimated 99% reduction for black abalone would result in Ne greater than 3,000, well above thresholds suggested to minimize genetic drift (66). Consistent with this, our 99% bottleneck simulations, though idealized, show that even at 0% recovery genomic erosion does not occur in the short term (Fig. 4A). For comparison, the timing and magnitude of myxoma virus spread in rabbits mirrors the WS bottleneck in black abalone, yet has also not significantly altered genetic diversity or population structure after roughly 70 generations (32). From this perspective, it is unsurprising that we see stability across our whole time series (Fig. 2), despite the archaeological evidence for human impacts on black abalone over millennia of harvest (67).

Our predictions of genomic erosion rest on accurate assessment of the WS bottleneck. The difference between a 99% bottleneck and a 99.9% bottleneck, the latter of which may better describe WS’s effect at select sites (68), has an impact on whether genomic erosion appears in the short term (Fig. 5). Under the more extreme scenario, only a 10% simulated recovery stabilizes genomic erosion by 100 generations postbottleneck. A recovery rate of this magnitude has been recorded at some sites (SI Appendix, Fig. S5) but is unlikely to be sustained for long. Related to this, our framework conflates Ne at the level of the species range with local (site-level) Ne (66), which may be lower (Limitations of this Study). Conversely, the initial WS impact may have been overestimated if the disease was more harmful to higher intertidal abalone visible to surveyors (69), or if gene flow from less impacted northern areas diminished the impact at southern sites (3). Increased migration, achieved through the indirect effects of population growth or individual translocation (70), may mitigate more serious impacts (12).

We also see evidence for natural selection on immune genes following the Withering Syndrome bottleneck (Fig. 5). Populations can adapt to intense selection pressures within just a handful of generations and leave genomic evidence of this process (31, 32, 71). Adaptation on these timescales is more likely to draw on standing genetic variation, and a selective sweep may involve multiple haplotypes carrying the adaptive allele (72). If adaptation indeed occurred in response to Withering Syndrome, the high standing genetic diversity of black abalone suggests that selection would have acted on multiple haplotypes. Consistent with this scenario, we observe multiple distinct low-diversity alleles of the SVEP1 outlier region, suggesting that multiple haplotypes carrying an adaptive mutation(s) may have undergone selection (Fig. 5D). In addition to this sweep signal, we observe large regions of high diversity and sequence divergence on scaffolds 7, 12, 14, and 19, which could indicate balancing selection (54). Short-term genomic signals of balancing selection are difficult to detect in contemporary samples (19, 54, 73), but the power of temporal samples to detect it remains unexplored except in experimental systems (74). Whether through balancing selection or other processes, the maintenance of diversity at immune loci is well documented in declining populations (75, 76). The high diversity outliers we detect are enriched for genes involved in immunity across mollusks (77, 78), in particular DMBT1 (Fig. 5E) (57). What remains unclear is how diversity at loci like DMBT1 has increased, as there has not been much time for de novo variation to arise. For DMBT1 specifically, the tight clustering of prebottleneck individuals with those as old as 1500 BP (Fig. 5 C vs. E), suggests that this result is not simply due to low sample size.

Other evidence for selection comes from the clinal shifts of the scaffold 4 and scaffold 9 chromosomal inversions (Fig. 6). Prior to this work, we had characterized the 31 Mb scaffold 4 inversion and its significant association with latitude, hypothesizing that it would be involved in local adaptation (7). With prebottleneck genomes, we now see that the scaffold 4 inversion has recently doubled in frequency at sites south of Pt. Conception—a common biogeographic barrier (61)—while remaining stable at northern sites (Fig. 6 and SI Appendix, Fig. S8). Surprisingly, a second inversion locus on scaffold 9 shows the same pattern, although a clinal fit is only slightly more supported than a stable fit (SI Appendix, Table S1). These parallel shifts could result from interactions between inversion alleles, as supported by our observation of significant present-day LD at southern sites (SI Appendix, Table S3). Inversions are common in wild populations (79, 80) but instances of epistasis between inversions are few (81, 82). Here, it is tempting to speculate that epistatic interactions are shaping WS adaptation at the southern, more WS-afflicted sites (3). Alternatively, rather than epistasis, both inversions may be responding independently to selection pressures that vary with latitude and have changed on the same timescale (80, 83). Regardless, the stark differences in inversion distributions across Pt. Conception and through time suggest that these loci may be important for local adaptation (84).

Limitations of this Study.

Our current simulation approach confirms many of the predictions of population genetic theory (85). A more thorough approach would incorporate selected alleles and non-Wright–Fisher population dynamics [i.e., (86)], would provide more confidence in our predictions for genomic erosion. Currently, incorporating these factors requires a forward simulation framework that records mutations and individuals over the course of simulation (87). This becomes prohibitively memory and time intensive when simulating the Ne of prebottleneck populations (3.0e5), and rescaling population parameters (e.g., Ne~, μ) for computational speed can produce bias (88). Moreso, many aspects of black abalone biology remain uncertain because attempts to rear multiple generations in captivity have been unsuccessful. Larval duration and its effect on dispersal is uncertain (6), and dispersal ability will impact time lag (28), local Ne (66), and recovery. Finally, the functions of the putative targets of selection also remain uncertain. Functional genomic data that might validate candidate genes is limited for black abalone and its closest congeners, all of which are endangered. Sequence homology to experimental mollusk systems remains the best way to speculate on these selection targets.

Implications for the Conservation of Black Abalone.

Future genomic erosion in black abalone is not guaranteed. If assessments of decline and recovery rates are accurate, even minimal growth may be sufficient to avoid the loss of genetic diversity, inbreeding depression, and loss of connectivity between sites. Knowing that genomic erosion is not a foregone conclusion can increase enthusiasm for conservation efforts (86) and better direct limited resources. Our genomic baselines show that black abalone lack genetic structure across much of their range and that the diversity that existed prior to Withering Syndrome is still represented today. Therefore, translocations of reproductive adults from growing sites (38) could be an effective way of fostering recovery that does not disrupt long term baselines. While we do not know the fitness consequences of the inversions that segregate across Pt. Conception, a conservative strategy might restrict translocations to sites on the same side of this boundary. Whether considering inversions or other loci like SVEP1, these results demonstrate that data connecting genotype to phenotype to fitness are sorely needed for this species. Experiments to collect such data (e.g., transcriptome responses to heat stress) should be prioritized alongside other conservation efforts. This will require increased support and flexibility from state and federal regulatory agencies. Finally, while we find these results encouraging for the genomic future of black abalone, we would like to acknowledge this is only one component of conservation. Other components, for example the preservation of high-quality habitat, are essential to ensure that healthy populations will be around for future generations (68).

Methods

Sample Collection.

For prebottleneck samples, 44 black abalone shells were loaned from the Malacology Collection at the Natural History Museum of LA County, ranging in age from 1914 to 1979. seven shells ranging from ~800 BC to ~1880 BC were loaned by Todd Braje at the University of Oregon Museum of Natural and Cultural History, and seven ranging from ~570 BC to 1770 BC loaned from the Amah Mutsun Tribal Band. All shells were identified as black abalone based on morphology. For postbottleneck samples, we combined 138 black abalone published in Wooldridge et al. (7) with 16 samples from Baja California (89).

Sequence Data Generation.

Shells.

All shell sequencing procedures were performed in a dedicated ancient DNA facility at UC Santa Cruz, following standard clean room criteria (Poinar and Cooper 2000). We first used a dremel to obtain a shell fragment from the anterior end of each specimen. For smaller or more fragmented shells, we used whatever material was available. We then sampled approximately 50 mg of shell powder after pulverizing each shell using a Mixer Mill MM 400 (Retsch). We incubated shell powder at room temperature in a 0.5% bleach solution followed by three washes in 1 mL of molecular grade water to remove contaminants [Korlević et al. 2015; Boessenkool et al. (64)]. Next, each powder sample was incubated overnight at 37 °C in digest buffer (1 mL: 0.45 M EDTA, 0.25 mg/mL Proteinase K). Finally, DNA was isolated using the silica column-based method described in Rohland et al. 2018 using Binding Buffer D with a final elution of 35μL buffer EBT. All extracts were quantified using a Qubit 4 (Invitrogen) and the Qubit 1× dsDNA HS assay kit (Invitrogen).

Single-stranded library preparation was performed for all extracts following the protocol outlined in Kapp, Green, and Shapiro (2021), with modifications as described in Nguyen et al. (2023). Single-stranded libraries were indexed and amplified in 50 μL reactions containing 20 μL preamplified library, 25 μL AmpliTaq Gold 360 Master Mix, 2.5 µL of 20 μM i7 indexing primer, and 2.5 µL of 20 μM i5 indexing primer. Libraries were amplified in a Bio-Rad C1000 thermocycler using the following 95 °C for 10 m, followed by 11 to 19 cycles of 95 °C for 30 s, 60 °C for 30 s, and 72 °C for 60 s, followed by 72 °C for 7 m. Postamplified libraries were purified using a 1.2× SPRI clean. Finally, libraries were visualized on an Agilent Fragment Analyzer. Libraries were pooled and sequenced on an Illumina NextSeq2000 at UC Santa Cruz (2 × 61 bp) to assess endogenous DNA content. Libraries with sufficient endogenous DNA content and complexity were sent for deeper sequencing at UC San Francisco Center for Advanced Technology on an Illumina NovaSeq 25B (2 × 100 bp).

Live samples.

All modern samples in this study derive from collection efforts published in Wooldridge et al. (7) and Delgadillo-Anguiano et al. (89). DNA extracts from the Delgadillo-Anguiano study were provided by Alicia Abadia-Cardoso, and libraries were generated following the NEBNext Ultra II FS DNA Library Prep Kit for Illumina (NEB) using the recommended protocol with Y-Adapters in place of the NEBNext Adapters. We incubated the samples for 5 to 6 min during enzymatic fragmentation, performed a single-sided 0.8× SPRI bead mixture prepared according to Rohland and Reich (2012). Libraries were amplified for seven cycles using dual unique indexes. Libraries were eluted in 21 μL 0.1× TE and quantified using the Qubit dsDNA HS Assay (Invitrogen) and an Agilent Fragment Analyzer. All libraries were sequenced on an Illumina NextSeq 2000 before being sent for greater sequencing at the UC San Francisco Center for Advanced Technology on an Illumina NovaSeq 25B (2 × 150 bp).

Alignment and Variant Calling.

We adapter-trimmed and merged overlapping reads from all shell libraries using fastp with default parameters (90). Merging was necessary given the smaller average fragment size of DNA in the shell libraries. Modern abalone libraries were also processed by fastp with the read merging step omitted. Next, we aligned all shell and modern samples to the primary haplotype of the black abalone reference genome [GCF_022045235.1; (91) with bwa mem using default parameters (92)]. We elected to use the bwa mem approach for all sample types to reduce batch effects, as analyses with mapDamage2 (93) showed minimal damage and fragmentation in the vast majority of shells (SI Appendix, Fig. S1). We marked and removed duplicates using sentieon driver --algo LocusCollector --fun score_info followed by sentieon driver --algo Dedup --rmdup. Finally, for all bam files we filtered for only primary alignments using sambamba view -F “not (unmapped or secondary_alignment or supplementary).”

We next generated per-sample gVCF files using sentieon driver –algo Haplotyper --emit_mode gvcf. We then performed joint genotyping on this set of gvcfs using sentieon driver –algo GVCFtyper --emit_mode ALL, which produced invariant + variant sites across the genome. Finally, we filtered these variant sites using GATK VariantFiltration (Van der Auwera et al. 2013). We performed initial filtering on SNPs and INDELs independently, excluding SNPs with QUAL < 30.0, QD < 2.0, FS > 60.0, MQ < 40.0, MQRankSum < −12.5, ReadPosRankSum < −8.0 or SOR > 3.0 and excluding INDELs with QUAL < 30.0, QD < 2.0, FS > 200.0, ReadPosRankSum < −20.0, SOR > 10.0 (94). For invariant sites, we filtered based on site quality (“QUAL > 30”). Finally, from this set of filtered variants we selected only biallelic SNPs at transversion sites. This transversion SNP dataset was used for all subsequent analyses unless otherwise specified.

Genomic Masks.

We generated three sets of genomic masks to remove regions that could bias downstream analyses.

First, we created a strict mappability mask with GenMap (95, 96). We ran GenMap with -K 60 -E 2 to score regions based on mapping of 60mers with up to two mismatches. From this, we created a mask to exclude regions with scores less than 1. These excluded regions represent areas of the genome where only the most damaged and fragmented shell libraries would be susceptible to mismapping. Second, we created a mask to exclude two putative chromosomal inversions identified in Wooldridge et al. (7). Third, we created masks based on depth. For this, we examined the distribution of site-level read depth in our variant + invariant site vcfs and excluded sites with less than 48 total reads (bottom 2.5% of sites) or more than 2,884 reads (top 2.5% of sites).

Pseudohaploid Genomes.

As a complementary way to account for differential coverage across samples, we called pseudohaploid and pseudodiploid genotypes for all samples. Genotypes were called by random read sampling at a set of previously ascertained sites. These sites included transversion SNPs in the modern (undamaged) samples which passed quality and mappability filters. Pseudohaploid genotypes were called with SAMtools (97) v1.9 using mpileup -B -q25 -Q30 and pileupCaller from sequenceTools v1.5.2 (https://github.com/stschiff/sequenceTools) with the --randomHaploid and --singleStrandMode options, which allowed for excluding genotypes calls potentially originating from ancient DNA damage. This enabled us to include transitions in downstream analyses. We created separate pseudohaploid genotype calls with major inversion regions included and excluded.

Reference Genome Bias.

In order to assess whether the choice of reference genome could be introducing bias in our analyses of the ancient shell genomes, we estimated D-statistics to test for excess allele sharing with the reference between all pairs of shell samples. Specifically, using the pseudohaploid calls, we calculated D-statistics of the form D(shell 1, shell 2; modern abalone, reference allele) for all combinations of shells, using several high coverage modern individuals from different populations. These statistics test whether either shell sample shares significantly more alleles with the reference, as compared to the modern abalone baseline. Negative skews in the resulting D-statistics may be interpreted as reference genome bias with shell 1, relative to the other shell used in comparison. D-statistics were calculated using ADMIXTOOLS2 (97).

We observed some evidence for reference bias in the oldest samples (SI Appendix, Fig. S10), several of which happened to be population outliers including the 1500 BP “San Vicente” site (Fig. 3). This reference bias was mitigated considerably by increased sequencing coverage (e.g., sample SC23.SF063 at 2× vs. SC23.SF059 at 37×; SI Appendix, Fig. S10). The majority of shells, in particular the younger samples from the “Commercial Fisheries” period (1910–1980), exhibited distributions of D-statistics centered on 0, indicating minimal relative reference bias. Because the only methods using these oldest, lowest-coverage samples—heterozygosity, fROH, and PCA—were a) specifically designed to account for DNA damage and sequencing depth, and b) corroborated by high coverage genomes from the same sites and time periods, we consider it unlikely that reference bias substantially influenced these results.

Population Genetic Analyses.

All following population genetic analyses implement different strategies based on sequencing coverage and the extent of DNA damage. For clarity we have included a table outlining the sample sets and data types used for each analysis (SI Appendix, Table S3).

Genetic Diversity and Population Structure.

Heterozygosity.

To analyze individual genetic diversity, we performed DNA-damage aware inference of heterozygosity in both full and 1× downsampled genomes using ROHAN (52). We first inferred DNA damage profiles for all shell samples using bam2prof -minq 20 -both. After confirming that the DNA damage profiles met those inferred by MapDamage2, we proceeded with heterozygosity inference using the command rohan –rohmu 1e-4 –tstv 1.06 –size 500000 –chains 10000, providing the deamination profiles with –deam5p and –deam3p. We then ran the same command for all modern samples, but without any deamination profiles. The choice of 1.06 for the transition-to-transversion ratio (--tstv) was informed by polymorphism in modern samples. For all commands, we also restricted our analyses to regions with good genome mappability (see Methods: Genome masks) with the –map flag.

PCA and inversion genotyping.

We performed a genetic principal components analysis (PCA) on three data types for 1) genotype likelihoods, 2) pseudohaploid genomes, and 3) filtered variant calls. Given that all individuals clustered by inversion genotype in Wooldridge et al. (7), and preliminary analyses with only high-coverage samples showed the same pattern, we had a strong a priori expectation that the inversion structure would be recapitulated in some form here.

However, an initial analysis of genotype likelihoods with PCAngsd (98) showed samples clustering strictly by whether the DNA came from shells or live individuals, regardless of coverage. The shell cluster also included the modern “control” shell that we sequenced, leading us to believe that this structure might be artificial. To test this, we generated pseudohaploid genomes as is common with low coverage degraded DNA samples (see above), then pruned variants in these genomes to randomly sample 1 transversion SNP every 1 kb. We used these variants as input to PCA with plink (99) and smartPCA (100). Complementary to this, we obtained the variant calls themselves at these same pruned sites to serve as input to plink and smartPCA as well. Both plink and smartPCA returned similar structure on both the pseudohaploid genomes and the variant calls. This structure showed the effect of the inversion genotype in modern samples and grouped shell + modern individuals from outlier populations together (Fig. 3), leading us to believe that this PCA approach was accurately capturing population structure.

We then used our PCA analyses to assign inversion genotypes to each individual. Following Wooldridge et al. (7), we generated PCAs based on variants from the scaffold 4 and scaffold 9 inversion loci. After confirming the three-cluster structure indicative of chromosomal inversions (79), we defined inversion heterozygotes (0/1) as individuals belonging to the central cluster (SI Appendix, Fig. S9). Finally, we defined inversion homozygotes (1/1) as individuals belonging to the less diverse outer cluster, and reference homozygotes (0/0) as those belonging to the more diverse outer cluster.

Pi, FST and DXY.

We generated diversity metrics from variant + invariant site VCFs with pixy (101), which estimates sequence diversity while accounting for the pitfalls in generating such estimates from heterogeneous data with high rates of missingness. We ran pixy --stats pi, dxy, fst on two sliding window 1) 50 kb × 25 kb sliding windows, and 2) 10 kb × 5 kb sliding windows. Both sets were masked based on mappability and overall sequencing depth (see above). For further comparison, we ran pixy for both the full sample set and high coverage subset, as well as on all polymorphisms and transversion polymorphisms only (SI Appendix, Fig. S3). Given the similarity in pi (π) between high coverage pre- and postbottleneck samples, we defined the statistic πpost-pre as the difference between postbottleneck π and prebottleneck π. We used πpost-pre for further inference of differences in selection between the two time periods (see below).

Genetic Load.

We estimated the load of deleterious mutations with snpEff (102) and snpSift (103, 104). For this, we used the NCBI-generated annotation for the black abalone reference genome (GCA_022045235.1). We restricted this analysis to only high coverage (>15×) samples with confident genotype calls at transversion sites, allowing us to capture relative load between historic and modern samples. After scoring these variants with snpEff, we selected for loss-of-function mutations (filter “(exists LOF[].PERC) & (LOF[].PERC > 0.9)”) and synonymous mutations (filter “ANN[0].EFFECT has ‘synonymous_variant’”) using snpSift. Finally, to obtain relative measures of load, we divided the loss-of-function variants by the number of synonymous variants for each sample.

Gene Content and Enrichment Analyses.

Particular candidate genes (e.g., Fig. 5 C–E) and enrichment results were only reported if they were significant for both 50 kb and 10 kb window sizes. We took shared genomic windows from our FST and πpost-pre analyses and intersected them with the H. cracherodii annotation from NCBI. We defined outliers as windows that were in the top 1% of FST and πpost-pre (e.g., balancing selection candidates) or the top 1% of FST and bottom 1% of πpost-pre (e.g., selective sweep candidates). We next added 5000 bp of to the start and end of each gene’s coordinates prior to this intersection in order to capture cis-regulatory regions under selection. We manually inspected each gene associated with an outlier window and also used Orthofinder v.2.5.5 (105) to identify orthologous gene families with the more studied Eastern oyster Crassostrea virginica (GCF_002022765.2). Finally, we constructed gene trees for individual outlier windows by first converting the VCF region to a phylip file (https://github.com/edgardomortiz/vcf2phylip) and performing phylogenetic inference with IQ-TREE 2 (105) using the GTR+I + G substitution model and default parameters.

Population Genetic Simulations.

We used msprime for all coalescent simulations (106). We simulated the long-term demographic history for black abalone inferred by Wooldridge et al. (8) and added a recent bottleneck and recovery (SI Appendix, Fig. S4). We also generated a population split concurrent with the bottleneck to explore the effects of the bottleneck on isolation between subpopulations. Bottlenecks of 99%, and 99.9% intensity were applied, and these bottlenecks were followed by growth occurring at rates of −10%, −1%, −0.1%, −0.01%, 0%, 0.01%, 0.1%, 1%, and 10%. We set these growth rates to return to zero if Ne reached the prebottleneck state (~3.0e5) or dropped to 100. All simulations were of 10 Mb chromosomes, and we generated 100 replicates of each parameter combination. We assumed a mutation rate of 8.60e−9 (8) and recombination rate of 1e−8.

For our short-term simulations (Fig. 4), we emitted diploid samples from this simulation at 10 generations prior to the bottleneck, the start of the bottleneck, and at intervals of 10 generations afterward. From these samples we computed 1) heterozygosity, 2) fraction of the genome in runs of homozygosity (FROH), 3) FST between the prebottleneck baseline and postbottleneck timepoints (Temporal FST), and 4) FST between postbottleneck subpopulations (Subpopulation FST) (SI Appendix, Fig. S4). Finally, we repeated the above procedure but calculated the same statistics at intervals of 100 generations postbottleneck in order to investigate long-term effects (SI Appendix, Fig. S6). All statistics were calculated with tskit and custom functions (107–109). Code for these steps is present in the scripts expo_split_shortterm.py and expo_split_longterm.py.

Inversion Clines and Associations.

Following inversion genotyping through PCA, we aimed to quantify spatial changes in the distributions of the scaffold 4 and scaffold 9 inversions over time. First, we applied a log transformation on the spatial data and a logistic transformation on inversion presence/absence following Westram et al. (110). Then, we fitted stable, linear, and sigmoid clines, based on the formula encoded in the R packages HZAR (111), to these data using a maximum likelihood search with mle2 function from the R package bbmle (112). For the search parameter space, we limited the lower and upper bounds of the inversion frequency to −1e-5 and 1e5, the lower and upper bounds of the cline center to 28- and 37-degrees latitude, and the lower and upper bounds of the cline width at 0.1 and 10 degrees latitude. Following maximum likelihood estimation, we determined the best fitting model of the three using Akaike’s Information Criterion (AIC) via the R stats function AIC (SI Appendix, Table S1).

Complementary to this cline estimation, we tested for two-locus Hardy–Weinberg equilibrium and linkage disequilibrium between each inversion locus using a chi square test. Specifically, we calculated population-level allele frequencies for each inversion by timepoint and position relative to the Pt. Conception barrier (i.e., Prebottleneck and North). From these frequencies, we calculated the expected genotype counts and derived the chi2 statistic using sum((data$obs_counts − data$exp_counts)^2/data$exp_counts). We then evaluated the significance of this statistic using the base R function pchisq(...,df = 6, lower.tail = F). All code for these analyses can be found in cline_fitting.R.

GBIF Analysis.

We aimed to quantify what proportion of all museum abalone specimens consisted of shells. To do so, we downloaded all Specimen records under TaxonKey “Mollusca” from GBIF on 21 October 2025. We removed entries with no listed preparation and filtered for entries with “Haliotis” in the species field. We then summed the “IndividualCount” field for all specimens with a hit for the search term “shell|dry|dried|concha|valve” in the “Preparations” field.

Supplementary Material