Unveiling the Ecological Applications of Ancient DNA From Mollusk Shells

Der Sarkissian, Clio; Möller, Per; Hofman, Courtney A.; Ilsøe, Peter; Rick, Torben C.; Schiøtte, Tom; Sørensen, Martin Vinther; Dalén, Love; Orlando, Ludovic

doi:10.3389/fevo.2020.00037

ORIGINAL RESEARCH article

Front. Ecol. Evol., 03 March 2020

Sec. Paleoecology

Volume 8 - 2020 | https://doi.org/10.3389/fevo.2020.00037

This article is part of the Research Topic Applied Uses of Ancient DNA View all 14 articles

Unveiling the Ecological Applications of Ancient DNA From Mollusk Shells

$\r\nClio Der Sarkissian*$ Clio Der Sarkissian^1*

Per Möller²

Courtney A. Hofman^3,4

Peter Ilsøe⁵

Torben C. Rick⁶

Tom Schiøtte⁷

Martin Vinther Sørensen⁷

Love Dalén^8,9

Ludovic Orlando^1,5

¹Laboratoire d’Anthropologie et d’Imagerie de Synthèse, UMR 5288, CNRS, Université Paul Sabatier, Toulouse, France
²Department of Geology, Quaternary Sciences, Lund University, Lund, Sweden
³Laboratories for Molecular Anthropology and Microbiome Research, Stephenson Research and Technology Center, Norman, OK, United States
⁴Department of Anthropology, The University of Oklahoma, Norman, OK, United States
⁵The GLOBE Institute, University of Copenhagen, Copenhagen, Denmark
⁶Department of Anthropology, National Museum of Natural History, Smithsonian Institution, Washington, DC, United States
⁷Natural History Museum of Denmark, University of Copenhagen, Copenhagen, Denmark
⁸Department of Bioinformatics and Genetics, Swedish Museum of Natural History, Stockholm, Sweden
⁹Centre for Palaeogenetics, Stockholm, Sweden

The shells of marine mollusks represent promising metagenomic archives of the past, adding to bones, teeth, hairs, and environmental samples most commonly examined in ancient DNA research. Seminal work has established that DNA recovery from marine mollusk shells depends on their microstructure, preservation and disease state, and that authentic ancient DNA could be retrieved from specimens as old as 7,000 years. Here, we significantly push the temporal limit for shell DNA recovery to ≥100,000 years with the successful genetic characterization of one Portlandia arctica and one Mytilus mussel sample collected within a dated permafrost layer from the Taimyr Peninsula, Russia. We expand the analysis of ancient DNA in carbonate shells to a larger number of genera (Arctica, Cernuella, Crassostrea, Dreissena, Haliotis, Lymnaea, Margaritifera, Pecten, Ruditapes, Venerupis) from marine, freshwater and terrestrial environments. We demonstrate that DNA from ancient shells can provide sufficient resolution for taxonomic, phylogenetic and/or population assignment. Our results confirm mollusk shells as long-term DNA reservoirs, opening new avenues for the investigation of environmental changes, commercial species management, biological invasion, and extinction. This is especially timely in light of modern threats to biodiversity and ecosystems.

Introduction

Applications of ancient DNA (aDNA) to ecological studies are plentiful, especially now that high-throughput DNA sequencing (HTS) technologies make it possible to generate genome-wide data at the scale of populations for the last million years (Orlando et al., 2013; Der Sarkissian et al., 2015; Leonardi et al., 2016; Nielsen et al., 2017; Fages et al., 2019; Narasimhan et al., 2019). By allowing high-resolution genomic history reconstructions, aDNA has the power to reveal how changes in the environment (e.g., climate, pollution levels), or interactions with other species and populations (e.g., hybridization, admixture, invasions), have potentially affected past organisms or populations through adaptation, distribution shifts, or extinction (e.g., Green et al., 2010; Lorenzen et al., 2011; Orlando and Cooper, 2014; Palkopoulou et al., 2015). Although most ancient DNA studies have focused on mammal bones and teeth (Green and Speller, 2017), new knowledge in ecology and evolution could be gained by applying ancient genomics to shell-producing mollusks (Coutellec, 2017). These animals indeed grow their carbonate shells in an incremental fashion, simultaneously recording biological and environmental data (Fortunato, 2015; Steinhardt et al., 2016; Butler et al., 2019), which can be tracked through a multitude of proxies, including morphology, isotopes (e.g., ∂¹⁸O_shell, ∂¹³C_shell), trace elements and metal composition (e.g., Ba, Cd, Cu, Mo, Pb, U, Zn). The rich fossil record of mollusk shells, thus, provides an equally rich environmental archive of seawater paleo-temperature and salinity (Surge et al., 2003; Chauvaud et al., 2005; Hiebenthal et al., 2012; Vokhshoori and McCarthy, 2014; Reynolds et al., 2016; Black et al., 2017), pollution levels (Vander Putten et al., 2000; Liehr et al., 2005; Pérez-Mayol et al., 2014), stress (Hiebenthal et al., 2012; Trivellini et al., 2018), infection history (Paillard et al., 2004; Trinkler et al., 2010), as well as food availability and primary productivity (Lartaud et al., 2010; Sadler et al., 2012).

Despite recent advances in ancient genomics and the ecological relevance of mollusks, no mollusk genomic time transect has ever been produced. The potential of applying HTS to ancient marine mollusk shells has, however, been assessed for the first time in a recent study that revealed their metagenomic content (Der Sarkissian et al., 2017). This study identified species with aragonitic shell microstructures, such as Mytilus mussels, Haliotis abalones, Arctica islandica quahogs, and Venerupis/Ruditapes clams, as good candidates for mollusk DNA content and/or preservation, with the longest DNA survival observed for a Mytilus specimen dated to ∼7,000 years Before Present (year BP; Der Sarkissian et al., 2017). This timeframe allows for the examination of many important questions relative to ecological changes and the genomic responses of mollusks, as it is characterized by several climatic transitions (Marcott et al., 2013), e.g., the Little Ice Age in the North Hemisphere (∼14th century to ∼1850 C.E.) or post-industrial global warming (Mann et al., 1999, 2009). During this period of intensifying human impact, mollusks have been exposed to increasing pollution levels that have resulted from rising human population densities, industrial activities, extensive farming, fishing and transportation (Jackson, 2001). They have also been commercially exploited as resources for food, nacre or pearls (Fortunato, 2015), and have been displaced over long distances intentionally for aquaculture purposes, or unintentionally (e.g., through ship transportation) (Carlton, 1999). Although successful DNA recovery may prove challenging, investigating mollusk shell DNA preservation beyond the 7,000 year BP limit holds the potential to provide valuable information about the impact of earlier climatic shifts in the Quaternary, such as the Pleistocene glacial/interglacial/interstadial oscillations, the transition into the Holocene, and subsequent increased human environmental impacts (Stephens et al., 2019). Other important areas of research include the applicability of ancient shell genomics to freshwater and terrestrial mollusks (Lydeard et al., 2004) and determining if ancient mollusk shell DNA enclose useful phylogenetic and population structure information.

To further unlock the potential of ancient mollusk shell DNA for ecological studies, we apply an aDNA HTS approach to a range of samples that extends previous assessments to significantly older specimens, to a larger number of species, as well as to mollusks of the freshwater and terrestrial niches. We then illustrate the importance of ancient mollusk DNA for understanding the impact of environmental change, extinction and biological invasion. We also explore the limitations of such investigations and evaluate the research potential of each genus examined, including genera of commercial interest.

Materials and Methods

Samples

For aDNA analyses, we selected three fossils identified as mollusk shells, sampled in a river-cut sediment succession along the lower reaches of the Bol’shaya Balaknya River (site BBR17 in Möller et al., 2019b, with close-by subsites 17A and 17B), situated in the southern portion of the Taimyr Peninsula, Russian Federation (N73° 37,084′; E105° 38,178′; Supplementary Table S1). This location is south of the limits of repeated Kara Sea-based glaciations during the last glacial cycle (the Zyryanka, Marine Isotope Stage (MIS) 5d-2) (Svendsen et al., 2004; Möller et al., 2015, 2019b). The samples were all found in a distinct marine silty clay designated sediment unit A, overlaid by fluvial sand and gravel (sediment unit B) (see sediment log in Möller et al., 2019b). The marine unit A sediment is rich in bivalve and gastropod shells (Supplementary Table S2), including subarctic species such as Buccinum undatum, Mytilus edulis, and Macoma balthica, not present today in the Kara Sea (to the north) and the Laptev Sea (to the east), thus suggesting interglacial conditions at sediment deposition (Möller et al., 2019b). Three shells (Portlandia arctica) from the unit A sediments were dated by Electron-Spin Resonance (ESR) to 101–105 kyr (± 9–12 kyr; averaged here to ≥100 kyr). The above-lying unit B fluvial sediments were dated to their deposition by Optically Stimulated Luminescence to 42–43 kyr (± 3–4 kyr), which is within the Middle Zyryanka (MIS 3), while redeposited mollusks (P. arctica) within the unit B sediments, eroded from underlying marine unit A, gave ESR ages of 122 and 123 kyr (± 15 kyr). Möller et al., 2019b thus concluded that the unit A marine sediment can be placed with high confidence to their deposition within the Karginsky interglacial (MIS 5e), equivalent to the North West European Eemian. Site BBR 17 is located ∼50 km from the ocean today, but paleoenvironmental reconstructions showed marine inundation of the Taimyr lowlands in front of the retreating ice margin at the transition between the Taz glaciation (MIS 6) and into the Karginsky interglacial (MIS 5e), with sea level at its maximum reaching levels in excess of 80 m above present sea level in this area (Möller et al., 2019b). The shells subjected to aDNA analyses were sampled from frozen marine unit A sediment (permafrost) at site BBR 17A and have been kept frozen since the time of collection in 2010 to optimize DNA preservation post-excavation.

An additional 39 mollusk shell specimens from nine genera (Arctica, Cernuella, Crassostrea, Dreissena, Haliotis, Lottia, Lymnaea, Margaritifera, Mytilus) were obtained from museum and laboratory collections. When available, genus- or species-level identification based on morphological examination, as well as collection location and date were retrieved from the information curated together with the specimens (Supplementary Table S1). Samples were selected to represent genera from marine, freshwater and terrestrial niches, with particular relevance to important ecological questions, such as those related to economical exploitation, environmental changes, biological invasions and extinctions, as well as to optimize comparison with available genomic resources.

DNA Extraction, Shotgun DNA Library Construction and High-Throughput Sequencing

We carried out shell DNA extraction, shotgun DNA library construction and amplification setup following (Der Sarkissian et al., 2017) in conditions strictly limiting DNA contamination at the ancient DNA facilities of the Centre for GeoGenetics, University of Copenhagen, Denmark. In order to monitor DNA contamination during pre-PCR laboratory work, every step was simultaneously performed on non-template Extraction Blank Controls (EBCs; Supplementary Table S1).

Fragments of the shells’ ventral margin were first decontaminated for 10 min under agitation in one volume of 1% sodium hypochlorite, before being rinsed three times in three volumes of distilled water, air-dried, and reduced to powder with a mortar and pestle (Der Sarkissian et al., 2017). We then performed DNA extraction based on the method developed in Yang et al., 1998; Gamba et al., 2014, 2016. Shell powder (112–2,990 mg) was incubated overnight at 37°C in 15 mL of digestion buffer (0.45 mM EDTA, 0.5% N-laurylsarcosyl, 0.25 mg/mL proteinase K) under constant mixing. Centrifugation at 3,000 RPM for 2 min allowed separating remaining solids from the liquid phase, which was then concentrated to a volume of ∼200 μL using an Amicon Ultra-15 30 kDa centrifugal filter unit (Merck Millipore) at 3,000 RPM for 50–60 min. The retentate was subjected to the MinElute PCR Purification kit (Qiagen), and purified DNA was obtained in a final elution volume of 60 μL EB buffer in presence of Tween (0.05% final concentration, hereafter, EB-Tween).

Double-indexed blunt-ended DNA libraries were constructed following a protocol modified from Orlando et al., 2013; Seguin-Orlando et al., 2013, where adapters carried an identifying 7 bp-index and were paired in a combination unique to each library (Rohland et al., 2015). We used the NEBNext Quick DNA Library Prep Master Mix Set for 454 (New England Biolabs) with a starting volume of 42.6 μL shell DNA extract in reaction volumes of 50 μL for end-repair (12°C for 20 min then 37°C for 15 min) and ligation (20°C for 20 min; 0.5 μM Illumina adapter final concentration), and 25 μL for fill-in (37°C for 20 min then 80°C for 20 min). After end-repair and ligation, the reaction mixes were purified using the MinElute kit with elution volumes of 30 and 20 μL EB-Tween, respectively.

For each library, we used real-time PCR (qPCR) to estimate the optimal number of amplification cycles to obtain DNA amounts compatible with Illumina sequencing, while minimizing library clonality. We performed two independent qPCRs per library in 20 μL of the following reaction mix: 1 μL 1/20 diluted DNA library (in EB buffer), 2.5 units AccuPrime DNA polymerase, 1× AccuPrime mix (Thermo Fisher Scientific), 1 mg/mL BSA, 0.2 μM primer inPE1.0 (5′-AATGATACGGCGACCACCGAGATCTACACTCTTTCCC TACACGACGCTCTTCCGATCT-3′) and 0.2 μM of an Illumina 6 bp-indexed (“I”) primer (5′-CAAGCAGAAGACGGCA TACGAGATIIIIIIGTGACTGGAGTTCAGACGTGTGCTCTTC CG-3′), and 0.8 μL of 1:4:2000 ROX:SybR:DMSO DNA dye mix (Thermo Fisher Scientific). qPCRs were performed on a LightCycler 480 Real-Time PCR System instrument (Roche Applied Science) with the following conditions: activation at 95°C for 5 min; 40 cycles of: denaturation at 95°C for 15 s, annealing at 60°C for 30 s, elongation at 68°C for 30 s; final elongation at 68°C for 5 min.

We analyzed qPCR data with the LightCycler Software 4.0 and the second derivative maximum method to estimate Ct values used as a relative measure of DNA concentration for each library. We found that libraries built from EBCs differed in Ct values by 4–6 additional cycles from those constructed from ancient shells. As the corresponding 16–64-fold increase in the amount of DNA in ancient samples relative to EBCs was not large enough to confidently rule out a significant impact of laboratory-derived DNA contamination, we subjected EBC libraries to the same PCR amplification and sequencing protocols as ancient shell libraries (see below).

Each double-indexed library was PCR-amplified using one Illumina IS4 primer and one unique 6 bp-indexed primer, thus, resulting in the introduction of a third, external index. The same conditions as for qPCR were applied, apart from non-diluted DNA and total reaction volumes of 5 and 25 μL, respectively, and from the absence of ROX:SybR:DMSO. We derived the optimal number of cycles to apply from the qPCR results on a per-library basis (14–29 for ancient samples and 34–35 for EBCs; Supplementary Table S1). In order to limit clonality, all libraries were first amplified for 12 cycles, before MinElute-purified reaction mixes (25 μL EB elution volume) were re-amplified as above for the required remaining number of cycles and in four independent reaction volumes (25 μL each) per library. Amplified libraries were purified using the Agencourt AMPure XP system (Beckman Coulter) and 25 μL of bead solution. Libraries eluted in 25 μL EB were quantified on a 2200 TapeStation Instrument (High Sensitivity D1000 Screen Tape; Agilent).

All triple-indexed libraries were pooled together in equimolar proportions and sequenced in paired-end mode on an Illumina HiSeq4000 platform at the Danish National High Throughput DNA Sequencing Centre. Raw sequencing data were deposited on the European Nucleotide Archive (ENA) public database (project PRJEB35671). Previously published ancient mollusk DNA data from Der Sarkissian et al. (2017) were retrieved from ENA (project PRJEB20113) and reanalyzed here.

Post-sequencing DNA Read Processing

Post-sequencing, we de-multiplexed DNA reads and only retained those displaying the expected, unique combination of indexes, i.e., unlikely to represent chimeric sequences. We used PALEOMIX v1.2.13 (Schubert et al., 2014) to trim adapters and collapse overlapping pair-end mates with AdapterRemoval2 v2.3.0 (Schubert et al., 2016) as described in Fages et al. (2019), i.e., applying default parameters and quality/length filters –minlength 25 –trimns –trimqualities –minquality 2.

Mapping to Mitochondrial and Nuclear Reference Sequences

Initial molecular taxonomic assignment was performed by mapping shotgun DNA reads to 145,457 reference sequences of the mitochondrial cytochrome c oxidase subunit I gene 5′-extremity barcode region (COI-5P) compiled for the Mollusca phylum from the Barcode of Life Database (BOLD¹). Building on the results, reads were then mapped to reference sequences for mitochondrial complete genomes, nuclear genome or transcript assemblies, when available for the identified genus/species (Supplementary Table S3). Some of the genera studied show Doubly Uniparental Inheritance, where the mitochondrial genome of the mother (or F-genome) is transmitted to both the somatic and gonadic cells of the female offspring, but only to the somatic cells of the male offspring, the gonadic cells of which get their mitochondrial genome from the father (M-genome). As a consequence, we used both F- and M-genomes as reference sequences for mapping (Supplementary Table S3). The end of the circular reference mitochondrial genome sequences was extended using the first 30 bp of the sequence in order to take the circularity of the mitochondrial genomes into account. We used the BWA v0.5.9 aln alignment command (Li and Durbin, 2009) in PALEOMIX (Schubert et al., 2014) with default parameters and disabled seeding as recommended in Schubert et al., 2012. We kept all reads mapping to the BOLD Mollusca database, whereas only reads showing mapping qualities ≥30 against the mitochondrial and nuclear DNA reference sequences were retained. We removed duplicated non-collapsed reads with MarkDuplicates in Picard Tools version 1.119² and duplicated collapsed reads with the PALEOMIX FilterUniqueBAM Python script (Schubert et al., 2014). Average depth-of-coverage (coverage) was calculated using PALEOMIX and considering only unique high-quality reads (after quality filtering and duplicate removal). When multiple reference sequences were available for a given genus, we retained the one yielding the highest coverage for each sample.

Complete Mitochondrial Genome Assembly and Maximum Likelihood Phylogenetic Reconstructions

For samples showing maximal coverage ≥3× when mapping against complete mitochondrial genome reference sequences, we built complete mitochondrial genome consensus sequences from which we constructed maximum-likelihood phylogenies. This was achieved using the perl script wrapper.pl (for circular genomes) as well as the C + ⁣ + programs bam2prof, endoCaller and log2fasta within the schmutzi pipeline (Renaud et al., 2015), taking into consideration rates of cytosine deamination at the 5′- and 3′- read ends, and only retaining bases with quality ≥30. For each sample, we called a first consensus from the rescaled BAM file, then re-mapped all shotgun reads to the newly built consensus with PALEOMIX as described above, and called a final consensus from the resulting rescaled BAM with schmutzi.

We then constructed alignments from the consensus sequences for ancient samples obtained here and in Der Sarkissian et al. (2017), as well as all sequences previously published for the genus/species of interest. We extracted the sequences of 12 genes (CYTB, COX2, NADH1, NADH4, COX3, NADH2, NADH3, COX1, ATP6, NADH4L, NADH5, NADH6) and aligned them independently with PRANK (Löytynoja and Goldman, 2005). For the analyses of Mytilus sp. sequences, we added the sequences of the 12S ribosomal RNA gene (12S rRNA), and 23 transfer RNAs (tRNA). Maximum-likelihood trees were reconstructed using the program PhyML 3.0 (Guindon et al., 2010) with an approximate likelihood-ratio test (aLRT) and a Shimodaira-Hasegawa-like procedure for branch support estimation. Using the SMS heuristic procedure for model selection implemented in PhyML 3.0 (Lefort et al., 2017), we assessed the best substitution models for the concatenated alignments: the Generalized Time-Reversible model, with invariant sites and 4 gamma site-rate categories for the Mytilus DNA alignment, and the Variable Time substitution model with invariant sites, a gamma site-rate distribution and alignment-based estimation of amino-acid frequency equilibrium for the mollusk protein alignment.

Population Structure Inference Based on Mitochondrial DNA

For all samples positively identified through comparison with the BOLD marker database and showing an average coverage ≥3×, a COI-5P consensus sequence was built as described for complete mitochondrial genomes. The obtained consensus was then compared to all the sequences available in the BOLD database for the corresponding genus/species using Median Joining Network (epsilon = 0) in POPART v1.7 (Leigh and Bryant, 2015) and the vouchers’ geographical origin as curated in the BOLD database metadata. When available, published mitochondrial sequences from modern populations for species of interest were retrieved for comparison to ancient sequences using the same procedure. These comprised Mytilus trossulus (Breton et al., 2006; Smietanka et al., 2010, 2013; Zbawicka et al., 2014b; Śmietanka and Burzyński, 2017), A. islandica (Glöckner et al., 2013), Ruditapes decussatus (Cordero et al., 2014; Sanna et al., 2017), R. philippinarum (Cordero et al., 2017), and Crassostrea angulata (Grade et al., 2016).

Population Structure Inference Based on Nuclear DNA Data

The ≥100 kyr Tx101A Mytilus sample from the Taimyr Peninsula could be compared to modern Mytilus populations for which genomic data are publicly available (Fraïsse et al., 2016, 2017). This comparative dataset was previously obtained through capture-based enrichment sequencing targeting 4.3 Mb DNA sequences assembled from M. edulis BAC clones and M. galloprovincialis cDNAs in 1,269 contigs (Fraïsse et al., 2016, 2017). It is composed of 80 individuals from ten populations: M. trossulus from the Baltic Sea (N = 8) and the North American Atlantic coast (N = 7), M. edulis from the North American Atlantic coast (N = 11), as well as from the Wadden Sea (N = 8) and the Bay of Biscay (N = 8; respectively, outside and inside the Atlantic hybrid zone with M. galloprovincialis), M. galloprovincialis from Portugal (N = 6) and Brittany (N = 8; respectively, outside and inside the Atlantic hybrid zone with M. edulis), as well from the Occidental (N = 8) and Oriental Mediterranean basins (N = 8), and M. platensis from the Kerguelen Islands (N = 8). In order to avoid computational analysis biases, the raw data published in Fraïsse et al. (2016, 2017) were re-mapped to the contig reference sequence at the same time as raw data for the Tx101A ancient sample and EBCs. The same BWA mapping parameters were used as described above for both ancient and modern samples, except that the options for mismatch penalty (-M 2), gap open penalty (-O 3) and minimum seed length (-k 10) were added as in Fraïsse et al. (2016, 2017), and that seeding was disabled for Tx101A and EBCs.

Low depth-of-coverage for Tx101A when mapped against the Mytilus contigs precluding accurate genotype calling, subsequent analyses where carried out using genotype likelihood in ANGSD v.0.929 (Korneliussen et al., 2014). Genotype likelihoods were estimated with the following options and filters: -doMajorMinor 1 -doMaf 1 -remove_bads 1 -uniqueOnly 1 -minMapQ 30 -C 50 -minQ 30 -baq 1 -skipTriallelic 1 -SNP_pval 1e-6 -doHWE 1 -maxHWEpval 0.05. All analyses were performed either considering all substitution types, or excluding transitions (-rmTrans 1). We used PCAngsd v.0.98 (Meisner and Albrechtsen, 2018) and the option -minMaf 0.05 to perform Principal Component Analysis (PCA) and determine the optimal number of principal components to describe the Mytilus population structure using the minimum average partial (MAP) test. Admixture proportions were estimated using NgsAdmix (Skotte et al., 2013) for a number of clusters equal to the previously determined optimal number of principal components +1 (Meisner and Albrechtsen, 2018), i.e., four, here.

Post-mortem DNA Damage Characterization

We used mapDamage v2.0.1 (Jónsson et al., 2013) to characterize ancient shell DNA in terms of fragmentation patterns, nucleotide mis-incorporation and fragment size distributions. On the basis of 100,000 Markov Chains Monte Carlo iterations, posterior distributions were obtained for the following damage parameters: rates of deamination in double strands (δ_D) and single strands (δ_S), probability of reads not terminating in overhangs (λ, transformed into 1/λ–1, a proxy for the length of overhanging regions). Quality scores of bases likely to represent damaged read positions were then rescaled using default parameters.

Assessing the Potential of Ancient Mollusk DNA for Ecological Studies

We provided a summarizing visualization of the potential for ecological studies of each mollusk genus investigated. In genus-specific radar graphs, we represented: their molecular preservation potential (oldest authenticated mollusk DNA data, DNA recovery rate), the availability or unavailability (encoded as 1 and 0, respectively) of comparative datasets (BOLD barcodes, complete mitochondrial genomes, nuclear transcripts and genome assemblies), and their relevance to ecological questions with regards to the research themes “Environment,” “Commercial” exploitation, and biological “Invasion” (assessed on the basis of bibliography and encoded as 1, relevant and 0, not relevant), and “Extinction.” An “extinction score” was given for each vulnerability category defined by The IUCN Red List of Threatened Species (2019): i.e., 1 = “least concerned,” 2 = “near threatened,” 3 = “vulnerable,” 4 = “endangered,” 5 = “critically endangered,” 6 = “extinct in the wild,” 7 = “extinct,” and a genus extinction score was calculated as the average score across species listed for each genus.

Results

Mitochondrial DNA-Based Taxonomic Identification of ≥100 kyr Mollusk Shells

We tested for long-term preservation of aDNA in marine mollusk shells by constructing shotgun DNA libraries using extracts from three mollusk shell samples (Tx100, Tx101A, Tx103) dated to ≥100 kyr. For the three libraries, qPCR showed copy-numbers compatible with HTS, which yielded 653,177–19,314,607 collapsed and non-collapsed pair-end mate reads after trimming and quality filtering (Table 1 and Supplementary Table S1).

TABLE 1

Table 1. Sample information and sequencing statistics.

Taxonomic identification was achieved by mapping against the BOLD Mollusca mitochondrial COI-5P marker reference sequences. No hit was obtained for Tx100, precluding any classification for this specimen. On the basis of mitochondrial DNA, Tx101A could be identified as M. trossulus since 98.6% of the mapped nucleotides aligned to markers of this species, plus 0.9 and 0.5% aligning to markers of the Mytilus genus and the Mytilida order (Supplementary Table S4). Tx101A had been previously identified as M. edulis on the basis of shell morphology in Möller et al., (2019a,b), but M. trossulus and M. edulis can be difficult to distinguish due to their similarity in shape and size (Innes and Bates, 1999; Riginos and Cunningham, 2005). The morphological assignment of the Tx103 sample to P. arctica could be confirmed as 27.3% of the mapped nucleotides align to markers of this species, and 72.7% to the Nuculanida order (Supplementary Table S4).

We attempted to gain more confidence in our taxonomic assignment, and retrieve more specific information about genomic affinities of Tx101A and Tx103 with modern individuals (e.g., at the population level), by comparing their COI-5P sequences with M. trossulus and P. arctica sequences available in the BOLD database. Median Joining Network analyses, however, proved uninformative due to the limited phylogenetic resolution of the COI-5P marker, the scarcity of geographic information provided for BOLD vouchers (Tx101A), and the lack of comparative P. arctica entries (Tx103) (Supplementary Figure S1).

Mitochondrial and Nuclear Genomic Affinities of a ≥100 kyr Year-Old Mussel

In order to achieve a more robust identification and a characterization of the Tx101A specimen with more resolution, we performed phylogenetic and population structure analyses based on HTS read alignments to all complete mitochondrial genomes published to date, as well as to genome-wide contigs for Mytilus. In absence of such public reference resources for P. arctica, or any close relative taxon, further genomic characterization of the Tx103 specimen was not achievable. The analyses provided further evidence for the mitochondrial identification of Tx101A as M. trossulus. The maximal coverage of the mitochondrial genome was, indeed, obtained when aligning HTS reads to the sequences of M. trossulus F-genomes or recently masculinized F-genomes (24.7–25.7× maximum for GenBank accession number GU936625; Supplementary Table S5). Other Mytilus mitochondrial genomes led to 10–9,267-fold lower coverage (Supplementary Table S5). This pattern was robust to relaxing or strengthening mapping stringency (-n, BWA edit distance option) (Supplementary Figure S2).

A Maximum-likelihood tree reconstructed from concatenated mitochondrial genome sequence partitions confirmed the phylogenetic placement of Tx101A within the clade of M. trossulus (recently masculinized) F-genomes, where it occupies a basal position (Figure 1A). Patterns of coverage are consistent with the phylogenetic position of Tx101A, i.e., the more distant to M. trossulus the mapping reference genome, the lower the coverage (Figure 1B and Supplementary Table S5). When only a 1,171 bp concatenated sequence partition of mitochondrial genes and transfer RNAs from Tx101A was compared to data from 187 present-day M. trossulus individuals (Breton et al., 2006; Smietanka et al., 2010, 2013; Zbawicka et al., 2014b; Śmietanka and Burzyński, 2017), the outlying position of Tx101A was confirmed. Closest relatives equidistant to the Tx101A haplotype were found in the broadly distributed present-day populations of the Baltic, the eastern and western coasts of the North Atlantic and the North Pacific (Figure 1C).

FIGURE 1

Figure 1. Phylogenetic and population affinities of Tx101A within present-day Mytilus. (A) Maximum likelihood phylogenetic tree built from a concatenated alignment of Mytilus mitochondrial sequences for 12 protein-coding genes, the 12S rRNA gene and 23 tRNAs (total length 14,277 bp). The sex and species of modern individuals were determined based on morphology and nuclear DNA markers. Branch support ≥0.70 is shown as Shimodaira-Hasegawa approximate Likelihood Ratio Test. *Recently masculinized mitochondrial genomes. (B) Average depth-of-coverage when mapping against each Mytilus mitochondrial reference genome considering unique high-quality reads (≥30). (C) Median-Joining network built from a 1,171 bp concatenated alignment of mitochondrial genes (COX3, ND2) and transfer RNA (Ser, Met) from 187 present-day M. trossulus individuals. (D) Principal Component Analyses of the top two principal components in a dataset comprising 11,882 variable sites (considering transversions only) for 80 individuals. The proportion of the variance explained by each component is indicated for each axis. (E) Ancestry proportions estimated considering four ancestral populations (transversions only). Individual color code is the same as in (D).

As hybridization can occur amongst members of the M. edulis complex (M. trossulus, M. edulis, M. galloprovincialis), we examined the nuclear genomic background of Tx101A in order to validate its taxonomic identification and/or estimate potential admixture proportions. PCA of 11,882 variable sites (considering transversions only; 29,440 variable sites considering all substitutions in Supplementary Figure S3) in 80 individuals and along the optimal number of three components (as determined by the MAP test in PCAngsd) confirmed that Tx101A belongs to the M. trossulus population cluster, clearly segregating from the other clusters of the M. edulis complex and of the M. platensis outgroup (Figure 1D and Supplementary Figure S3). This result was confirmed by clustering analyses considering four ancestries, in which a genome-wide affinity of Tx101A for M. trossulus modern populations and no admixture with another M. edulis complex or outgroup population were detected (considering transitions alone; Figure 1E and Supplementary Figure S4).

Authentication of the DNA Data Obtained From the ≥100 kyr Old Mollusk Shells

Considering that the age of the Tx101A and Tx103 samples dramatically pushes back the limit of marine mollusk shell DNA preservation from ∼7,000 years to ≥100 kyr, particular attention was paid to the line of evidence supporting data authenticity: (1) Tx101A and Tx103 were excavated at the same site, but yielded DNA reads that were identified as representative of different Bivalvia sub-classes; (2) our molecular classification confirms previous morphological identification of Tx101A as Mytilus and Tx103 as P. arctica (Möller et al., 2019b, a); (3) the detection of M. trossulus and P. arctica at the BBR 17A site was ecologically coherent with these being considered as (sub-) arctic species of the littoral or deeper marine environments (Möller et al., 2019b); (4) none of the HTS datasets obtained for the extraction blank controls (EBCs) led to the identification of M. trossulus or P. arctica DNA sequences. Compared to those observed for Tx101A and Tx103, significantly lower coverage were attained when mapping EBC reads against reference sequences for: M. trossulus/P. arctica BOLD COI-5P (no hit), Mytilus complete mitochondrial genomes (≥2,335-fold coverage decrease), the M. galloprovincialis genome assembly (≥16,563-fold coverage decrease), Mytilus sp. genome-wide contigs (≥5,704-fold coverage decrease); (5) for both Tx101A and Tx103, high-quality HTS reads aligned to the various reference sequences displayed molecular signatures characteristic of post-mortem DNA fragmentation and cytosine deamination-induced mis-incorporation. When considering the COI-5P alignment, Tx103 showed reads of ∼66.6 bp average length and a 9.4% C-to-T mis-incorporation rate at the 5′-read extremities. As for Tx101A, when considering mitochondrial genome/nuclear contig alignments, HTS reads exhibited ∼82.8/78.5 bp reads, maximum 3.9/6.4% C-to-T mis-incorporation rates, post-mortem deamination rates in overhangs 4.9/5.2-larger than in double stranded DNA, and overhanging ends of 3.0/2.0 bp (1/λ – 1), as expected. In addition, a 10–bp periodicity in the size distribution was observed in the nuclear DNA of Tx101A, but not in its mitochondrial DNA, a specific pattern that was previously proposed to reflect nucleosome protection in authentic and degraded ancient nuclear DNA (Pedersen et al., 2014; Figure 2 and Supplementary Table S6); and (6) all phylogenetic and population structure analyses of mitochondrial and nuclear DNA were consistent with Tx101A being part of M. trossulus. Combined, these results provide compelling evidence that the DNA obtained from ≥100 kyr mollusk shells is authentic.

FIGURE 2

Figure 2. Post-mortem DNA damage patterns. (A) For Tx101A, rates of C-to-T nucleotide mis-incorporation at 5′-read ends, (B) rates of G-to-A at 3′ read-ends. (C) Size distribution of reads mapping to mitochondrial, and (D) nuclear reference sequences. (E) For Tx103, rates of C-to-T nucleotide mis-incorporation at 5′-read ends, (F) Rates of G-to-A at 3′ read-ends. (G) Size distribution of reads mapping to mitochondrial reference sequences. Only unique high-quality (MQ≥30) reads were considered and mapDamage v.2 was run on the full alignment (100,000 iterations).

Expanding the Exploration of Ancient DNA Recovery From Mollusk Shells

We next extended the assessment of DNA recovery from ancient mollusk shells to 39 new specimens, 14 of which led to shotgun HTS datasets containing 2,477,750–30,000,220 quality-filtered collapsed and non-collapsed pair-end mate DNA reads (Table 1 and Supplementary Table S1). In subsequent analyses, these were added to 32 previously published ancient shell data (Der Sarkissian et al., 2017) and the two ≥100 kyr Taimyr samples described above to build a dataset comprising 48 specimens from 13 species and 9 genera, with 73% of the samples dated to the last 300 years of the post-industrial era (Figure 3 and Supplementary Table S1). Success rates were variable across genera, with the highest rates observed for Ruditapes/Venerupis sp. (100%, N = 20), Haliotis sp. (70%, N = 7) and A. islandica (57%, N = 12). No successful DNA recovery was achieved for the Lottia gigantea (N = 1), Margaritifera margaritifera (N = 3), and the terrestrial Cernuella virgata specimens (N = 2), thus limiting investigations to marine (N = 13) and freshwater individuals (N = 3) in this study. We found no indication that the overall success rate depends on the amount of shell powder analyzed, as specimens yielding successful and unsuccessful DNA recovery did not show significant differences in input shell median weight (Wilcoxon signed-rank test, p-value = 0.3).

FIGURE 3

Figure 3. Geographical location and age of the shell samples screened for DNA. (A) Map showing the number and average age of specimens by geographical location. (B) Chronology by genus (Supplementary Table S1). “kyr,” kilo years.

Mitochondrial DNA-Based Taxonomic Identification

Mapping to the Mollusca database of BOLD COI-5P markers provided sufficient read alignments for taxonomic identification, as all 48 samples, but the LEFB Pecten sp. shell, could be classified at least to the genus level (Figure 4 and Supplementary Table S4). A. islandica, C. angulata, Haliotis tuberculata, R. decussatus, R. philippinarum, and Dreissena polymorpha specimens were unambiguously identified at the species level, and even subspecies level for H. tuberculata coccinea and H. tuberculata tuberculata. We noted, however, 6–30% false positive rates for H. rufescens, V. corrugata, and Lymnaea stagnalis at the genus level. Through resequencing, we also assembled 29 additional complete mitochondrial genomes with coverage ≥2.8-fold that were used to reconstruct a maximum-likelihood tree of mitochondrial genes. The L. stagnalis and V. corrugata specimens could not be included in this analysis as no reference sequences are available for these species. The obtained tree retraces known phylogenetic relationships within gastropods and bivalves, but also confirms the taxonomic identification based on COI-5P sequences (Figure 5). Interestingly, with five new R. decussatus mitochondrial genomes assembled here from ancient specimens, we could confirm a previous hypothesis based on one modern R. decussatus mitochondrial genome that R. decussatus is more closely related to the Paphia sp. than to R. philippinarum despite being named as members of the same genus (Ghiselli et al., 2017).

FIGURE 4

Figure 4. COI-5P barcode-based identification of Mollusca taxa from shell DNA. Proportion of total nucleotides mapped to each reference marker combining all datasets for samples of the same genus and considering unique high quality (≥30) reads and reference markers covered by ≥150 bp (Supplementary Table S4).

FIGURE 5

Figure 5. Maximum likelihood phylogenetic tree of mollusk mitochondrial genes. Tree built from a 60-individual alignment of 12 concatenated gene sequences. Sequences assembled from shell DNA in this study are highlighted in red and their coverage is indicated in the tip label, GenBank accession numbers are provided for previously published sequences. Branch support ≥0.70 is shown as Shimodaira-Hasegawa approximate Likelihood Ratio Test. The subtrees for the R. decussatus and R. philippinarum clades were enlarged and branches collapsed in order make branching order more apparent.

Population Affinities Between Ancient and Present-Day Mollusks

Little information about affinities between ancient and modern populations could be retrieved from the phylogenetic tree owing to the limited number of published complete mitochondrial genomes in mollusks. Although more mitochondrial COI-5P sequences are available, analyses proved poorly informative for most genera due to the paucity of sequenced present-day individuals (Pecten maximus) and corresponding geographical metadata (P. maximus, L. stagnalis). Additionally, many ancient haplotypes were found similar to the most abundant and widely distributed haplotypes in modern populations partly due to the limited resolution of the COI-5P marker (C. angulata, D. polymorpha, R. decussatus; Supplementary Figure S1). The mitochondrial phylogenetic tree could, however, show a segregation within R. decussatus between the ∼30 year BP Adriatic Sea MURp individual from the 400–5,500 year BP Atlantic Ocean auzay1B, auzay3B, lmc1B, and lmc3B individuals (Figure 5). Within the other clam species R. philippinarum, both the mitochondrial gene tree and the COI-5P network showed a differentiation between the 1988 C.E. POSp French specimen and the 1983–2012 C.E. LAN1p, LAN2p, LAN3p, NEG, MATp, AKKp and KORp samples from France, Japan and Korea (Figure 5 and Supplementary Figure S1). For A. islandica quahogs, we compared our 73–297 year-old dataset to 24 COI-5P sequences from the BOLD database (Supplementary Figure S1) as well as to 3,029-bp long mitochondrial sequences for 20 modern individuals (Glöckner et al., 2013). Modern quahogs belong to two main clades: one comprising individuals from Iceland, the Baltic and the North Seas, and one composed of individuals from Iceland and the Baltic Sea, with the Icelandic and the Baltic Sea populations displaying a larger diversity than the North Sea population (Glöckner et al., 2013). We found the same pattern in the past as our sequences from ancient quahogs fell within the two clades described previously (Supplementary Figure S5). More present-day data in the form of complete mitochondrial genome sequences are needed to rule out possible sampling biases in both the number and geographical origin of the sequenced individuals.

Authentication of Ancient Mollusk DNA Data

Authenticity of the aDNA retrieved from the mollusk shells was supported by: (1) congruent molecular and morphological identification, at least at the genus level; (2) increased coverage when mapping against the complete mitochondrial genome reference sequence for the identified species (Supplementary Figure S6); (3) phylogenetic placement of the COI-5P barcodes and mitochondrial genes (Figures 4, 5 and Supplementary Figure S1); and (4) the presence of typical molecular signatures of mitochondrial DNA post-mortem damage in the form of short fragment sizes (45.1–95.2 bp), increased rates of C-to-T mis-incorporation rates at read 5′-ends (0.0–32.2%), increased mis-incorporation rates in single- versus double-strand contexts ∂_S/∂_D (2.3–70.3), and decreased values for the overhang-length proxy 1/λ – 1 (0.1–4.7). Similar patterns were observed in reads mapping against nuclear reference sequences: DNA fragmentation (41.9–107.8 bp), C-to-T mis-incorporation rates at read 5′-ends (0.0–17.4%), ∂_S/∂_D (0.1–184.6), and 1/λ – 1 (0.2–16.3) (Supplementary Table S6).

Discussion

Preservation of ≥100 kyr DNA From Marine Mollusk Shells in Siberian Permafrost Sediments

In this work, we broaden the assessment of mollusk DNA recovery from ancient shells in terms of time scale, species, and applications for ecological studies. The successful retrieval of aDNA from ≥100 kyr shells significantly extends the temporal scale of mollusk aDNA analyses, as the oldest specimens that had previously yielded aDNA were ∼7,000 year-old Mytilus shells collected from shell middens in Denmark (Der Sarkissian et al., 2017). The constantly frozen Siberian permafrost marine sediments from which the ≥100 kyr shells were retrieved (Möller et al., 2019b, a) may have provided exceptionally favorable DNA preservation conditions minimizing damage induced by both water and microbial activity (as reviewed in Pedersen et al., 2015). As examples of long-term DNA preservation in similar conditions, the oldest remains for which animal aDNA could be recovered in Siberia were >63.5 kyr and ∼50 kyr mammoth hair, and tooth/femur (Gilbert et al., 2008; Pečnerová et al., 2017), and, outside Siberia, a >560 kyr horse bone from the Yukon Territories, Canada (Orlando et al., 2013). Although earlier studies demonstrated permafrost as one of the best environments for long-term preservation of animal DNA in keratin- and calcium phosphate-based remains, we show here for the first time that permafrost also preserves calcium carbonate-based remains over significant timescales. Past marine diversity was formerly recovered from <43 kyr sediments that yielded DNA from protists, including haptophytes and foraminifera, which are other organisms producing calcium carbonate exoskeletons surviving or dissolving within sediments (Boere et al., 2009; Coolen et al., 2013; Lejzerowicz et al., 2013; Pawłowska et al., 2014; Morard et al., 2017; More et al., 2018; Armbrecht et al., 2019).

Recovery of Degraded Ancient Mollusk DNA Molecules From 13 Marine and Freshwater Species

Compared to the previous assessment of eight marine species (Der Sarkissian et al., 2017), we demonstrate here that aDNA can be recovered from an extra five species (P. arctica, H. rufescens, H. cracherodii, D. polymorpha, and L. stagnalis), including two freshwater species. Considering the relevance of mollusks to ecological questions about the terrestrial niche (Lydeard et al., 2004), aDNA content of terrestrial mollusk shells is worth further investigation. In the summarizing visualization of the potential of each mollusk genus investigated for ecological studies, radar-plot maximal surface areas indicate Mytilus, Haliotis, and Arctica as the three genera with the most potential for ecological studies based on aDNA (Figure 6). Overall, we observed a DNA recovery success rate of 61.25% for the 80 shells examined (Figure 6), and endogenous contents of 0.09–33.26% for those genera with available reference nuclear genome assemblies (Crassostrea, Dreissena, Haliotis, Lymnaea, Mytilus), and we estimated that no less than 71,000 sequencing reads were required here to obtain 1× coverage of 16,358–18,653 bp mitochondrial genomes. In line with these results suggesting low numbers of shotgun-sequenceable DNA templates in ancient mollusk shell extracts, we were not successful at recovering DNA from three 42–191 year BP M. margaritifera shells, although positive nuclear Small Tandem Repeat (STR) PCR amplification could be achieved from fresh dry shells of the same species by Geist et al. (2008). Our negative results could be explained by the particular storage conditions of the specimens, extensive post-mortem DNA degradation and/or concentrations of recoverable DNA below the sensitivity threshold of our methodology. Target-enrichment protocols may help circumventing these limitations in the future. We, however, caution that such methods will result in the loss of the vast majority of off-target DNA molecules. Their benefits in terms of information gain should be weighed against the risk of heritage loss for each (rare) ancient specimen considered. In addition, the type of capture should carefully be selected and designed to maximize the overlap with datasets generated from comparative modern or ancient populations. Importantly, neither capture nor shotgun HTS aDNA data can, at this stage, provide much information about population dynamics in cases where modern population data available for comparison are nuclear STR genotypes. These indeed typically encompass >100 bp regions (e.g., Knott et al., 2003; Cruz et al., 2005; Gardeström et al., 2008; Feldheim et al., 2011; Van Wormhoudt et al., 2011; Beldade et al., 2012; Morvezen et al., 2013; Borrell et al., 2014; d’Auriac et al., 2017; Jiang et al., 2018) that cannot reliably be typed from aDNA extracts due to their characteristic high fragmentation (here, 41.9–107.8 bp).

FIGURE 6

Figure 6. Radar plot summarizing application potential of ancient shell DNA for mollusk genera. Marine genera are represented in blue shades, freshwater genera in green shades and the terrestrial genus Cernuella in purple. “Max. Age” refers the age of the oldest sample for which DNA could successfully be recovered; “kyr,” kilo years; “year BP,” years Before Present.

Limits Associated With the Availability of Modern Mollusk DNA Data for Comparative Analyses

We observed a substantial heterogeneity among the examined genera in the availability of comparative mitochondrial sequences, complete mitochondrial genomes, nuclear transcripts and genomes (Figure 6). As expected for mostly non-model species, the present study highlights a general paucity of published mitochondrial and nuclear genome assemblies, population-scale data, as well as geographical metadata submitted to barcode databases (BOLD, GenBank), with some mollusk populations better described at the population genomic level, i.e., Crassostrea sp., Mytilus sp., Haliotis sp., and Pecten sp. (Zbawicka et al., 2014a, 2018; Fraïsse et al., 2016; Mathiesen et al., 2016; Wenne et al., 2016; Gutierrez et al., 2017; Harney et al., 2018; Wilson et al., 2018; El Ayari et al., 2019; Masonbrink et al., 2019; Paterno et al., 2019; Vendrami et al., 2019a, b). It is worth noting that the phylogenetically robust data reported here significantly increases the number of released complete mitochondrial genome sequences by 6-fold for R. decussatus, five-fold for R. philippinarum, three-fold for A. islandica, and D. polymorpha, and finally two-fold for H. rufescens, H. tuberculata, and P. maximus. Although mitochondrial DNA is very informative for genus-level classification and investigating the evolution of organellar genomes, it should be stressed that species-level taxonomic identification and population affinity inference should be confirmed by nuclear sequences for those mollusks in which interspecies hybridization is observed, e.g., Crassostrea and Mytilus. Importantly, as genomic data production is expected to increase significantly in the near future, it is crucial that appropriate metadata are made easily accessible to the research community, so that published results can be reproduced and future analyses can incorporate ever-growing datasets. Despite limitations associated with the degraded nature of mollusk shell aDNA and the occasional scarceness of comparable modern population data, phylogenetic and population structure signals could nevertheless be recovered, thus emphasizing the promising potential of mollusk shell aDNA for ecological studies.

Potential of Ancient Mollusk DNA for Ecological Studies

Our results directly hint at taxonomic identification based on mitochondrial barcoding, which could reveal helpful in cases where morphological examination of fragmentary or undiagnostic shell remains preclude classification. Mollusk aDNA analyses are compatible with shotgun- or mitochondrial metabarcoding-based (Bush et al., 2019) reconstructions of past communities from bulk samples (e.g., Murray et al., 2013; Grealy et al., 2015; Seersholm et al., 2018) or sediment cores (e.g., Giguet-Covex et al., 2014; Pedersen et al., 2016), which could be of relevance to the study of shell middens. Temporal changes in ecosystems detected through such approaches could reveal interesting insights into the timing and impact of environmental changes, biological invasions and extinctions.

Potential of Ancient Mollusk DNA for Studying Environmental Changes

Time series DNA data from mollusk shells hold the potential to illuminate patterns of species or population distribution and diversity. These can be compared to past global and local paleo-environments reconstructed through stable isotope, metal or trace element analyses of the very same shells in order to characterize the temporal and spatial scales of biological responses to climatic events, in particular prior to access to instrumental data or historical records (reviewed in Jones et al., 2009). Marine bivalves such as A. islandica and P. maximus, for which we reported successful aDNA recovery in this study, have previously allowed sclerochronological reconstructions of past seawater temperature and salinity changes at millennial-length annual and seasonal resolutions, respectively (Surge et al., 2003; Chauvaud et al., 2012; Vokhshoori and McCarthy, 2014; Reynolds et al., 2016; Black et al., 2017). At a more local scale, the freshwater L. stagnalis, for which we report shell aDNA recovery here for the first time, is a model commonly used in ecotoxicology (Amorim et al., 2019). Its abundant and widespread distribution in temperate limnic systems and its low dispersal ability can be leveraged to investigate the impact of pollutants in freshwater systems, in particular pesticides, molluscicides, algaecides or industrial (petro-) chemicals pollutants (Amorim et al., 2019). Adding temporally sampled genomic data would allow for controlling the genomic background of L. stagnalis population subjected to toxic stress and provide baselines prior to the use of pesticides and chemicals (Coutellec et al., 2013; Bouétard et al., 2014).

Integrating aDNA and environmental information provides new perspectives to better understand mollusk responses to global and local changes. These are expected to be complex, heterogeneous through time, space and across species. Many gaps still remain in our knowledge of the principles and mechanisms driving these responses; for example, the relative contribution of migration (routes and speed) and in situ tolerance in the form of phenotypic plasticity (morphology, behavior, ecophysiology, reproductive strategies) or evolutionary adaptation (speciation, selection) to the survival of species facing environmental changes (e.g., climate, pollution) of varying magnitude and rate.

Potential of Ancient Mollusk DNA for Studying Species of Commercial Interest

Mytilus, Haliotis, Crassostrea, and Ruditapes species represented 73% (8,816,367 tons) of the world marine mollusk production in 2017 (FAO, 2017). Their production can be managed in large-scale hatcheries and has resulted in massive translocation across large geographical scales. For example, Haliotis were commercially fished in Europe in the 19th century, which led to drastic wild stock depletions and the implementation of management plans and regulations for fishing and cultivation (Huchette and Clavier, 2004). European Crassostrea has had a distinct history, as the Portuguese oyster C. angulata was first introduced unintentionally from Asia, probably in the 16th century, and cultivated until it underwent high mortalities associated with iridovirus outbreaks in the 1960–1970s (Grizel and Héral, 1991). It was then substituted by the more resistant Pacific cupped oyster sister species C. gigas with distributions in Asia overlapping those of native C. angulata (Grizel and Héral, 1991). Haliotis and Crassostrea in North America have similar histories of overfishing and stock depletion across the 19 and 20th centuries (Braje et al., 2009, 2016; Rick et al., 2016). As for the Manila clam R. philippinarum, it has a natural distribution on the West coast of the Pacific from the Philippines to Russia. According to historical records, R. philippinarum was first unintentionally introduced from Japan to North America through oyster “hitchhiking” in 1936 C.E. (Quayle, 1964), and then imported from there into Europe partly to overcome the over-exploitation, irregular yields, recruitment failures, and outbreak of bacterial infections in the endemic R. decussatus in 1972–1974 C.E. (Flassch and Leborgne, 1992; Sanna et al., 2017). Human-driven displacements, intense harvest, over-exploitation, and controlled reproduction in hatchery are expected to have impacted exploited mollusk species, in particular with regards to genomic patterns of diversity (e.g., inbreeding) and/or admixture (Astorga, 2014), as well as allele frequencies for loci associated with economically important phenotypic traits (e.g., survival, growth rate, body weight, and yield) (Gutierrez et al., 2018). These could be tracked through time using the same mollusk shell aDNA approach as applied here, which could also potentially reveal the existence of undocumented transfers.

Additionally, as human-mediated transfers may subject farmed mollusks to new pathogens, mollusk shell aDNA could also help track and monitor the spread of infections. For example, in R. philippinarum, high mortalities were observed in 1987 C.E. in Brittany from the so-called “Brown Ring Disease” (BRD) for which the bacteria Vibrio tapetis was identified as the etiological agent (Paillard et al., 1989, 2008; Borrego et al., 1996; Allam et al., 2000). DNA recovered from the shells of R. philippinarum showing substantial brown conchiolin deposits diagnostic of acute BRD (shells labeled as POSp and KORp in the present study; Supplementary Figure S1) was previously found in Der Sarkissian et al., 2017 to show highest affinity for the virulent V. tapetis strain RP2-3, whereas DNA from shells displaying weak BRD infection (LAN1p, LAN2p, and LAN3p) showed highest affinity for the less virulent V. tapetis strain HH6087 (Der Sarkissian et al., 2017; Dias et al., 2018). Here, phylogenetic analyses of the R. philippinarum COI mitochondrial gene in specimens from Landéda, Brittany, reveal a segregation between the 1988 C.E. acute BRD specimen POSp falling within the mainly Asian diversity (China, Japan), and the 1983–1988 C.E. asymptomatic or weak BRD specimens NEG, LAN1p, LAN2p, and LAN3p closely related to mainly European and American individuals (Supplementary Figure S1). These results are in line with a tentative scenario, in which lower virulence V. tapetis strains may have already been present in Brittany before the 1987 C.E. virulent BRD outbreak, which might have been triggered by a second introduction of R. philippinarum stocks, possibly of more recent Asian origin, containing highly virulent strains. The presence of the 2003 C.E. KORp shell in Korea showing acute BRD and a COI sequence closely related to European and American clams could possibly reflect a later transfer from Europe to Korea, where BRD was first identified by molecular methods in 2006 C.E. (Park et al., 2006), and where stocks have underwent severe reductions due to overexploitation and coastal pollution (Cordero et al., 2017).

Potential of Ancient Mollusk DNA for Studying Biological Invasions

In a context of global trade, human activities greatly facilitate rapid species displacements over large geographical distances, whether it be intentionally (e.g., for cultivation in hatcheries or direct commercial purposes) or unintentionally (e.g., through waterways, ship and seaplane hull fouling, ballast waters, fishing) (Carlton, 1999). Thanks to a competitive advantage (e.g., long larval phases, rapid growth, early sexual maturity, high fecundity, broad salinity, and temperature tolerance), newly introduced species can proliferate from a point of introduction into new ecosystems in which they become dominant as invasive species (Valéry et al., 2008). Biological invasions are recognized as the second most important threat to biodiversity after habitat destruction due to their substantial negative ecological (e.g., food web alteration and extirpation of native species) and economic impacts [e.g., biofouling, blocking of water pipes (McCarthy et al., 1997)].

Belonging to the mollusk genera examined here, C. gigas and R. philippinarum threaten the European native species of Ostrea edulis and R. decussatus, respectively, and M. galloprovincialis has been listed as one of the “100 world’s worst invasive alien species” (Lowe et al., 2004). D. polymorpha is another well-identified invasive species for which shell aDNA could be retrieved. Today, D. polymorpha has a worldwide distribution and its invasive status is a recognized threat (Lowe et al., 2004). It initially spread from the Ponto-Caspian Basin to Europe in the mid 1800s (or before), and then reached North America from Europe in the mid 1980s (Kinzelbach, 1986; Carlton, 1999). Although it is one of the most studied invasive species with regards to ecological, toxicological and physiological aspects, population genomic data are currently scarce. Early detection of invasive species is of foremost importance in efforts to eradicate them, or at least slow their spread down (Williams et al., 2017). Designing management plans for the control of invasions can greatly benefit from the reconstruction of their history: geographical origin, number, and timing of the introduction(s), as well as dispersal and connectivity patterns. This can be informed by phylogeography in both the native and invaded areas. The timing of invasions, however, can be difficult to estimate with precision. One reason for this is that some mollusk (cryptic) invasive and native species show similar shell morphology due to their close relationships and phenotypic plasticity (Morais and Reichard, 2018). Some of the local species may also be impacted by hybridization with invasive sister species in semi-reproductive isolation (Harrison and Larson, 2014), which is facilitated by high fecundity and dispersal potential in broadcast-spawning aquatic mollusks Dreissena, Crassostrea, and Mytilus species (Voroshilova et al., 2010; Fraïsse et al., 2016; Gagnaire et al., 2018). Mollusk shell DNA is a useful tool to assess the presence and genomic landscape of invaders at given points in time and space. Considering the dynamic and stochastic nature of biological invasions, shell aDNA could help evaluate their impact on native populations and ecosystems, thus showing potential for future studies in ecology and conservation.

Potential of Ancient Mollusk DNA for Studying Extinctions

A total of 1,526 species of gastropods and bivalves are classified as threatened by The IUCN Red List of Threatened Species (2019). Among them, two species of Haliotis abalones of the North American Pacific coast have been classed as “endangered,” i.e., H. kamtschatkana, and “critically endangered,” i.e., H. cracherodii or black abalone. The latter was shown here to be amenable to aDNA analyses, similarly to another Pacific abalone species, H. rufescens or red abalone. Both species were found in shell middens and in artifacts at prehistoric Native American sites (Erlandson and Rick, 2008; Braje et al., 2009). H. cracherodii experienced critical population declines following intensive fishing in the 1850-1900s and then again after the peak in abalone farming of the 1950–1960s, which severely affected all abalone populations in the Pacific, including H. rufescens (Rogers-Bennett et al., 2002). Abalone population declines were aggravated by the recovery of sea otters (Enhydra lutris nereis and E. l. kenyoni) in the 20th century, abalone predators that had nearly gone extinct following overexploitation by the historical fur trade (Lee et al., 2016). High mortalities due to infections by the bacterium Candidatus xenohaliotis responsible for the Withering syndrome devastated abalone populations (Crosson et al., 2014) and contributed to the low genetic diversities and population differentiations observed in H. cracherodii and H. rufescens today (Gaffney et al., 1996; Burton and Tegner, 2000; Hamm and Burton, 2000). Conservation efforts are focused in restoring and protecting abalones in the Pacific: shore picking of H. rufescens is allowed only under strict rules, and fishing of H. cracherodii (commercial or recreational) have been suspended (Haas et al., 2019). Management and conservation plans would greatly benefit from shell aDNA analyses that could help better understand the dynamics of genomic diversity and distribution for the two species, the effects of exploitation, infection and predator recovery, and complement archeological studies aimed at helping in abalone restoration (Braje et al., 2009, 2016; Hofman et al., 2015).

Conclusion

In this study, we illustrate the potential of ancient mollusk shell mitochondrial and nuclear DNA analyses using geological, archeological and museum specimens from 13 species, spanning ≥100 kyr. Ancient mollusk DNA is of growing relevance for ecological studies as present-day societies are increasingly concerned by on-going global environmental changes (e.g., warming, pollution, sea level rise, ocean acidification, eutrophication) and associated biodiversity loss (Cardinale et al., 2012) continue to increase at an alarming pace. Some environmental stressors are expected to influence directly shell-producing mollusks in their ability to calcify, and hence, their integrity and vulnerability to predators (Kroeker et al., 2013), with some of the affected species playing the role of keystone species with large impacts on ecosystems. Our work highlights the importance of mollusk shell aDNA as a powerful and novel tool for ecological studies.

Data Availability Statement

The datasets generated for this study can be found in the European Nucleotide Archive public database (project PRJEB35671).

Ethics Statement

Ethical review and approval was not required for this study because no living animals were collected or examined.

Author Contributions

CD and LO conceived and coordinated the study. CD performed molecular and computational analyses, with input from LO. PM, CH, PI, TR, MS, TS, LD, and LO provided samples and reagents. CD wrote the article with input from PM, LD, and LO.

Funding

This work was supported by the Danish Council for Independent Research, the Natural Sciences (FNU, 4002-00152B), the Danish National Research Foundation (DNRF94), and the Initiative d’Excellence Chaires d’attractivité, Université de Toulouse (OURASI). LO has received funding from the European Research Council (ERC) under the European Union’s Horizon 2020 Research and Innovation Programe (Grant agreement no. 681605-PEGASUS). PM and LD received funding from the Swedish Natural Science Research Council (VR) (Grant agreements nos. 621-2008-3759 and 621-2008-4896, respectively), and expedition logistics, including for the Taimyr 2010 expedition, from the Swedish Polar Research Secretariat (SPRS).

Conflict of Interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Acknowledgments

We thank Martin Sikora at The GLOBE Institute, University of Copenhagen, for transporting samples. We also thank the PALEOMIX group at the University of Copenhagen and the AGES group at the Laboratory AMIS UMR 5288 (CNRS, University Toulouse 3), as well as the staff of the Danish National High-Throughput DNA Sequencing Centre and of AMIS for support.

Supplementary Material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fevo.2020.00037/full#supplementary-material

Footnotes

References

Allam, B., Paillard, C., Howard, A., and Le Pennec, M. (2000). Isolation of the pathogen Vibrio tapetis and defense parameters in brown ring diseased Manila clams Ruditapes philippinarum cultivated in England. Dis. Aquat. Organ. 41, 105–113. doi: 10.3354/dao041105

PubMed Abstract | CrossRef Full Text | Google Scholar

Amorim, J., Abreu, I., Rodrigues, P., Peixoto, D., Pinheiro, C., Saraiva, A., et al. (2019). Lymnaea stagnalis as a freshwater model invertebrate for ecotoxicological studies. Sci. Total Environ. 669, 11–28. doi: 10.1016/j.scitotenv.2019.03.035

PubMed Abstract | CrossRef Full Text | Google Scholar

Armbrecht, L. H., Coolen, M. J. L., Lejzerowicz, F., George, S. C., Negandhi, K., Suzuki, Y., et al. (2019). Ancient DNA from marine sediments: precautions and considerations for seafloor coring, sample handling and data generation. Earth Sci. Rev. 196:102887. doi: 10.1016/j.earscirev.2019.102887