INTRODUCTION
Speleothems, secondary cave carbonates such as stalagmites and flowstones, are among the most complete terrestrial paleoenvironmental archives that can be accurately dated back to ca. 600 ka (Cheng et al., Reference Cheng, Edwards, Shen, Polyak, Asmerom, Woodhead and Hellstrom2013). While speleothem-based inferences commonly include the reconstruction of hydrological, paleoclimatic, and biogeochemical cycles from seasonal to millennial time scales (Fairchild and Baker, Reference Fairchild and Baker2012), the use of speleothems as bio-archives remains largely unexplored (Blyth and Frisia, Reference Blyth and Frisia2008; Lechleitner et al., Reference Lechleitner, Dittmar, Baldini, Prufer and Eglinton2017; Luetscher et al., Reference Luetscher, Moseley, Festi, Hof, Edwards and Spötl2021).
In ventilated caves, the transport of aerosols favors incorporation of pollen grains and other particles into speleothems (Bastin, Reference Bastin1979; Dredge et al., Reference Dredge, Fairchild, Harrison, Fernandez-Cortes, Sanchez-Moral, Jurado and Gunn2013), and the remarkably resistant pollen walls consistently slow down degradation processes (Li et al., Reference Li, Phyo, Jacobowitz, Hong and Weng2019). Although the analysis of pollen, spores, and other palynomorphs (organic-walled microfossils) encapsulated in sedimentary deposits remains the principal technique used to reconstruct long-term vegetation changes (Faegri et al., Reference Faegri, Kaland and Krzywinski1989), species-level identification remains difficult (Richardson et al., Reference Richardson, Lin, Sponsler, Quijia, Goodell and Johnson2015) because several morphological features are shared within genera and sometimes even families. Moreover, pollen concentrations in these sedimentary archives are often low (Matley et al., Reference Matley, Sniderman, Drinnan and Hellstrom2020, Luetscher et al., Reference Luetscher, Moseley, Festi, Hof, Edwards and Spötl2021).
Thus, alternative methods based on DNA analysis, allowing the investigation of small amounts of biological particles and providing good to outstanding taxonomic resolution, hold great promise. Environmental DNA (eDNA) metabarcoding (Taberlet et al., Reference Taberlet, Bonin, Coissac and Zinger2018) has been successfully applied for pollen identification in different environments (e.g., Bell et al., Reference Bell, de Vere, Keller, Richardson, Gous, Burgess and Brosi2016; Leontidou et al., Reference Leontidou, Vokou, Sandionigi, Bruno, Lazarina, De Groeve and Li2021; Varotto et al., Reference Varotto, Pindo, Bertoni, Casarotto, Camin, Girardi, Maggi and Cristofori2021), but the exploitation of terrestrial archives has been restricted mainly to lake sediments, permafrozen soils, and ice cores (Jørgensen et al., Reference Jørgensen, Haile, Möller, Andreev, Boessenkool, Rasmussen and Kienast2012; Giguet-Covex et al., Reference Giguet-Covex, Pansu, Arnaud, Rey, Griggo, Gielly and Domaizon2014; ter Schure et al., Reference ter Schure, Bajard, Loftsgarden, Høeg, Ballo, Bakke and Støren2021). These sedimentary archives have given access to relevant paleovegetation records, but dating beyond the radiocarbon age-span remains problematic (Willerslev et al., Reference Willerslev, Hansen and Poinar2004). Occasionally, speleothems have also been used for molecular-based biological investigations, but limited to very specific applications (e.g., the study of microbial communities associated with biogeochemical processes; Dhami et al., Reference Dhami, Mukherjee and Watkin2018) or resulting in a very low number of detected organisms, with taxonomic complications due to the high degradation of DNA (Zepeda Mendoza et al., Reference Zepeda Mendoza, Lundberg, Ivarsson, Campos, Nylander, Sallstedt and Dalen2016; Stahlschmidt et al., Reference Stahlschmidt, Collin, Fernandes, Bar-Oz, Belfer-Cohen, Gao and Jakeli2019; see Discussion).
Here, we applied sedimentary ancient DNA (sedaDNA) metabarcoding (Pedersen et al., Reference Pedersen, Overballe-Petersen, Ermini, Sarkissian, Haile, Hellstrom and Spens2015) for identifying plant taxa in three 7600-year-old stalagmite sections from the Milchbach cave system (Switzerland, 1840 m asl; Luetscher et al., Reference Luetscher, Hoffmann, Frisia and Spötl2011). The periglacial Alpine environment and the reasonably fast calcite deposition rate (ca. 40 μm/yr) strictly limit DNA leaching, while a strong ventilation favors the incorporation of biological particles. Indeed, a remarkably high pollen concentration was found in Milchbach speleothem sections in a previous palynological analysis (564 palynomorphs/g; Festi et al., Reference Festi, Hoffmann and Luetscher2016). In this study, we applied a rigorous approach following a set of recently defined recommendations (Capo et al., Reference Capo, Giguet-Covex, Rouillard, Nota, Heintzman, Vuillemin and Ariztegui2021) to test if results from PCR-based sedaDNA metabarcoding are reliable and comparable to those of pollen spectra derived from the same speleothem sections by means of classical morphological analyses (Festi et al., Reference Festi, Hoffmann and Luetscher2016). We highlight several intrinsic limitations of the PCR-based amplicon sequencing technology, which can make mean its application to speleothem samples (and other sedaDNA matrices) incurs considerable uncertainties, and we suggest potentially more effective alternative strategies.
MATERIALS AND METHODS
Samples
We investigated for the presence of plant DNA in three coeval stalagmites (MB3, MB4, MB5-D) from Milchbach Cave, a periglacial cave system that opens along the western flank of the Upper Grindelwald Glacier valley, Switzerland (1840 m asl; see Fig. 1). U-series dating previously confirmed deposition 7000–8000 yr BP (Luetscher et al., Reference Luetscher, Hoffmann, Frisia and Spötl2011), and palynological analysis highlighted the presence of pollen entrapped in the calcite fabrics (Festi et al., Reference Festi, Hoffmann and Luetscher2016). Carbonate sub-samples of ~2–4 cm3 in size were cut along the stalagmite growth axis using a diamond wire saw. The surface of each section was cleaned in an ultrasonic bath using deionized water and stored in a sealed sterile bag until further processing. The taphonomy of Alpine speleothems, with peculiar chemism and steady deposition of calcite, together with the temperature conditions of Milchbach Cave (0–2°C; Luetscher et al., Reference Luetscher, Hoffmann, Frisia and Spötl2011) are potentially favorable to the long-term preservation of DNA (see Appendix S1 for more details on the cave system and speleothem samples).
Figure 1. Upper left: location of the sampling site. Upper right: schematic representation of the Milchbach cave system. Bottom: precise location of each speleothem sample (red) from where DNA has been extracted.
DNA extraction, PCR amplification and sequencing
All laboratory work was carried out at an ancient DNA-dedicated facility at Fondazione Edmund Mach, following the strict procedures required for ancient DNA studies to avoid contamination with modern samples (see Appendix S1: Sample preparation and pre-processing).
Speleothem specimens were made into a fine powder with a cryogenic grinding mill (using liquid nitrogen; 6770 Freezer/Mill, SPEX SamplePrep, Metuchen, New Jersey). During the procedure, each sample was placed in a closed grinding vial. DNA was extracted from the speleothem powder using the DNeasy PowerSoil DNA extraction kit (Qiagen, Hilden, Germany). In this protocol, cell lysis occurs by mechanical (beads-beating) and chemical processes, allowing extraction of DNA even from the toughest and most difficult samples. The final DNA extract was eluted in 100 μl. We processed one specimen for each speleothem section, except for MB5-D, for which two distinct specimens were available and were therefore processed separately. For each specimen, the obtained powder was split into two samples, from which DNA was extracted in two temporally separated extractions (except for MB3: the low quantity of starting material resulted in only one sample for extraction). Moreover, for each powder sample (i.e., for each single extraction), two replicates were performed using a different amount of starting material (speleothem powder): 250 mg and 500 mg. A scheme of processed samples, DNA extraction, and extraction replicates is reported in Table 1. Two negative controls (blanks) were included in each extraction in order to check for the presence of contaminants in reagents and/or in the laboratory environment.
Table 1. Replicates scheme: number and type of extraction and PCR replicates (repl.) per speleothem section and specimen. Suffixes E1, E2, and E3 identify different (temporally separated) DNA extractions, while suffixes a and b identify extraction replicates performed using 250 mg and 500 mg of speleothem powder, respectively.
For identifying plant DNA, a short fragment of the chloroplast gene trnL was amplified using the universal primers c-A49325 and h-B49466 [5′-CGAAATCGGTAGACGCTACG-3′ and 5′-CCATTGAGTCTCTGCACCTATC-3′] (Taberlet et al., Reference Taberlet, Coissac, Pompanon, Gielly, Miquel, Valentini, Vermat, Corthier, Brochmann and Willerslev2007) produced by Sigma-Aldrich (Milan, Italy). This trnL fragment has been previously tested in a methodological paper on DNA metabarcoding of airborne pollen (Leontidou et al., Reference Leontidou, Vernesi, De Groeve, Cristofolini, Vokou and Cristofori2018) and subsequently used in more recent applications of eDNA analysis of pollen in different Italian Alpine habitats (Leontidou et al., Reference Leontidou, Vokou, Sandionigi, Bruno, Lazarina, De Groeve and Li2021) and Alpine ice cores (Varotto et al., Reference Varotto, Pindo, Bertoni, Casarotto, Camin, Girardi, Maggi and Cristofori2021). PCR amplification mixture and thermal conditions are reported in Appendix S2. Band sizes and concentrations of PCR products were checked on a gel electrophoresis using QIAxcel with DNA High Sensitive cartridge (Qiagen, GmbH, Hilden, Germany) and analyzed by QIAxcel ScreenGel Software, in a separated post-PCR room. Multiple PCR replicates were performed for each specimen/extraction, according to the availability of extracted DNA after the PCR set-up (see Table 1). A maximum number of 16 samples was amplified for each separate PCR reaction, including extraction negative controls (blanks) and at least one PCR blank every eight samples, to track potential contaminations. PCR reaction products were gel purified using the MinElute Gel Purification kit (Qiagen, GmbH, Hilden, Germany) and eluted in 50 μl. Each individual PCR replicate was treated separately and processed as an independent sample for sequencing and bioinformatics analysis. Approximately 100 ng of the purified PCR product were used for library preparation using TruSeq DNA sample preparation kit V2 (Illumina Inc., San Diego, CA, USA), pooled in equimolar ratio, and sequenced using MiSeq Reagent Kit v3 in Illumina MiSeq platform.
Bioinformatics analysis and post-processing of metabarcoding data
Raw Illumina data were processed using the MICCA v. 1.7 (Albanese et al., Reference Albanese, Fontana, De Filippo, Cavalieri and Donati2015) pipeline and the derived operational taxonomic units (OTUs) were classified using the Sintax algorithm (Edgar, Reference Edgar2016) as implemented in VSEARCH (Rognes et al., Reference Rognes, Flouri, Nichols, Quince and Mahé2016; for the parameters used, see Appendix S3). Taxonomic assignment of OTUs was performed by comparison with a custom-made reference database including trnL DNA sequences for most of the Alpine seed plants (Spermatophyta). This local database was built starting from a previously developed trnL database (Leontidou et al., Reference Leontidou, Vernesi, De Groeve, Cristofolini, Vokou and Cristofori2018), which was focused on the flora of a southeastern Italian Alpine region, subsequently complemented with public sequences (NCBI GenBank database; National Center for Biotechnology Information, www.ncbi.nlm.nih.gov/genbank) for the other most frequent species of the Alpine flora (Aeschimann et al., Reference Aeschimann, Lauber, Moser and Theurillat2004).
After the bioinformatic analysis, a strict, multistep post-processing procedure was implemented on the dataset in order to remove negligible OTUs (choosing a stringent threshold: number of reads <100) and identify major and minor (potential) contaminants. Details on the post-sequencing removal of potential contaminations are described in Appendix S3 (Supplementary Information).
RESULTS
Eighteen out of 88 PCR replicates (20.45%) gave a positive outcome based on gel electrophoresis (at least one visible band) and were therefore sequenced, along with 42 negative controls. We obtained at least one PCR replicate with positive outcome for each of the three speleothem sections (MB3, MB4, MB5-D) and for all four processed specimens (two specimens were available for the MB5-D section). Thus, 4/7 of the extraction replicates (powder samples) resulted in the detection of at least one plant taxon (see Table 1 for details on replicates scheme). No systematic differences in PCR outcome were observed among extraction replicates performed from different initial amounts of calcite powder (250 mg versus 500 mg). After sequencing and bioinformatics filtering, a total of 1,353,919 merged reads were obtained and 15 operational taxonomic units (OTUs) were identified as present with a non-negligible (>100) number of reads (for details, see Appendix S3 in Supplementary Information).
Negative controls
We obtained 243,340 sequence reads in negative controls (average: 5794 ± 19,066), with three OTUs being recorded with non-negligible quantities, corresponding to the following taxa: Mediterranean cypress (Cupressus sempervirens L.), European yew (Taxus baccata L.), and alder (Alnus sp. Mill.) (Figure 2). Specifically, Mediterranean cypress was recorded in three PCR replicates of extraction blanks (with 43,175; 19,166; and 37,019 number of reads, respectively), all corresponding to the same extraction batch (E3). Alder was recorded in two PCR replicates of one extraction blank (E3; 51,680 and 33,352 number of reads). These two taxa were absent from all PCR negative controls. European yew was recorded only in one PCR blank (55,264 number of reads) and was not recorded in extraction blanks, despite the latter being amplified in multiple PCR replicates.
Figure 2. Results of sedaDNA metabarcoding. (a) Plant taxa recorded in the three speleothem sections from Milchbach Cave (Switzerland); only OTUs recorded with No. of reads >100 were considered (see Methods); (b) number of replicates and (c) number of negative controls analyzed (EXT and PCR refer to extraction and PCR blanks, respectively).
Overall, we recorded contamination in 1/24 PCR negative controls and in 3/18 PCR replicates of negative control of extractions, all of them deriving from the same extraction batch. Notably, the three taxa recorded from negative controls were absent from filtered speleothem samples (with a 100-reads threshold for negligible presence; see Appendix S3). Excluding these three OTUs from the dataset, the average number of reads per blank decreased to 49, showing no other relevant contaminants.
Samples
We obtained 1,099,305 reads from the speleothem samples (average: 61,072 ± 17,484), corresponding to 12 OTUs detected with non-negligible number of reads (>100). The 12 OTUs were assigned to 11 plant taxa (each one reported at the lowest reliable taxonomic level): (1) Pinaceae, Picea sp. (spruce); (2) Pinaceae, Pinus sp. (pine); (3) Abies sp. (fir); (4) Cupressaceae, Juniperus communis (common juniper); (5) Urticaceae, Parietaria sp. (pellitory; two OTUs); (6) Poaceae (grasses); (7) Fagaceae, Fagus sylvatica (European beech); (8) Plantaginaceae (the plantain family); (9) Oleaceae (the olive and ash family); (10) Betulaceae, subfamily Coryloideae, which includes the following genera: Corylus (hazel), Carpinus (hornbeam), and Ostrya (hop-hornbeam); and (11) Caprifoliaceae, subfamily Valerianoideae, which in the European Alps includes genus Valeriana, Valerianella, and Centranthus. All these taxa were considered absent from negative controls according to the adopted cut-off threshold (<100 reads), except for Fagus sylvatica, which was recorded with up to 150 reads in four blank PCR replicates. According to our criteria, the recorded taxa should therefore be considered reliable endogenous presences (i.e., not deriving from external contaminations), with the above-mentioned exception of Fagus sylvatica (see further discussion below). For all of the 18 sequenced PCR replicates, at least one plant taxon was detected. A synthetic representation of the recorded taxa composition for each stalagmite section is reported in Figure 2, together with number of extraction and PCR replicates (for details on number of reads/taxa for each PCR replicate, see Table S1 in Supplementary Information).
Four taxa were recorded in the MB3 section (in four positive PCR replicates, all corresponding to the same extraction): (1) Betulaceae, subfamily Coryloideae; (2) Caprifoliaceae, subfamily Valerianoideae; (3) Picea sp.; and (4) Fagus sylvatica. Each taxon was recorded in only one PCR replicate, except for Caprifoliaceae, subfamily Valerianoideae (two PCR replicates) and always with high number of reads (>40,000), except for Fagus sylvatica, which was recorded with 125 reads. Fagus sylvatica was recorded with comparable abundance in four negative controls and therefore, based on DNA data, cannot be considered a reliable presence in the MB3 section.
Two taxa were recorded in the MB4 section (one positive PCR replicate): (1) Betulaceae, subfamily Coryloidea (~69,000 reads) and (2) Abies sp. (733 reads). Abies sp., even if recorded with relatively low abundance, was however almost completely absent from all negative controls (max number of reads = 3).
The MB5-D section resulted in the highest rate of success both in terms of positive PCR replicates (13; representing 2/4 different extractions; i.e., powder samples) and plant diversity, with the following taxa: (1) Pinus sp.; (2) Picea sp.; (3) Juniperus communis; (4) Parietaria sp.; (5) Betulaceae, subfamily Coryloideae; (6) Poaceae; (7) Plantaginaceae; (8) Oleaceae; and (9) Fagus sylvatica. However, only 1–4 taxa were recorded in each of the PCR replicates. Pellitory (Parietaria sp.) was found in six PCR replicates, always with high number of reads (>10,000), except for one case (number of reads = 547). Common juniper (Juniperus communis) and pine (Pinus sp.) were recorded in five PCR replicates, always with number of reads >10,000, except for pine in two cases (but always >3000 reads). Spruce (Picea sp.) was recorded in two PCR replicates, always with number of reads >20,000. Each of the remaining taxa was recorded in only one PCR replicate, in all cases with a high number of reads (>40,000), except for Fagus sylvatica (European beech), which was recorded with 1825 reads. The abundance of Fagus sylvatica in this speleothem section is more than an order of magnitude higher than in negative controls, but its presence in blanks still prevents us from considering it an endogenous presence with high confidence.
DISCUSSION
Indirect authentication of sedaDNA sequences using palynomorphological data
To test the reliability of PCR-based metabarcoding of sedaDNA from speleothems, we focused on speleothems already known to contain high concentrations in biological particles. For two of the three 7600 yr-old Milchbach stalagmite (Luetscher et al., Reference Luetscher, Hoffmann, Frisia and Spötl2011) sections analyzed in this study, namely MB3 and MB5-D, palynomorphological-based data that were available from Festi et al. (Reference Festi, Hoffmann and Luetscher2016) allowed us to make a direct comparison between the two methodologies (see Table S2 in Supplementary Information for pollen counts of the recorded plant taxa). This comparison has the potential to reveal whether or not there is a similarity between the speleothems’ biological content determined from sedaDNA and that from pollen data. Considering the two sections together, all the taxa identified by target DNA barcode region (a small fragment of the chloroplast gene trnL, ~100–185 bp long) have also been detected from the palynomorphological analysis, with the only exception of Caprifoliaceae, subfamily Valerianoideae (found only in two PCR replicates of the MB3 sample).
Focusing on MB5-D, the most pollen-rich section (3767 pollen counts versus 256 for MB3), four of the nine most abundant plant taxa highlighted by the palynomorphological analysis (number of pollen counts >100) corresponded to those retrieved by sedaDNA (Figure 3): (1) Pinus spp.; (2) Corylus avellana (based on the taxonomic level reached by the target DNA barcode region we identified the presence of Betulaceae, subfamily Corlyloidea, which includes Corylus); (3) Poaceae; and (4) Fraxinus excelsior (DNA-based taxonomic level: family Oleaceae, represented in the native Alpine flora only by the genera Fraxinus and Ligustrum).
Figure 3. Speleothem section MB5-D: palynomorphological versus sedaDNA analysis. Plant taxa are ordered by pollen abundance (No. of counts); taxa that also were detected by sedaDNA are highlighted in red. Pollen counts are derived from Festi et al. (Reference Festi, Hoffmann and Luetscher2016). Note that, due to the different taxonomic resolutions, we adopted the following approximated correspondences: Corylus (pollen data) = Betulaceae, subfamily Coryloidea (sedaDNA); Fraxinus (pollen) = Oleaceae (sedaDNA). See text for further explanations.
In general, the palynological analysis resulted in detection of a higher number of taxa, compared with sedaDNA: 46 versus 9 in the MB5-D section. Many diverse and not mutually exclusive reasons for explaining this stark difference can be invoked. In morphological analyses, even a single pollen grain can be spotted and attributed to a specific taxon. In contrast, eDNA metabarcoding identifies taxa only if a discrete number of sequences belongs to a specific taxon. Furthermore, molecular investigations are limited by DNA degradation, while pollen walls are known to be one of the most resistant cellular portions (Brooks and Shaw, Reference Brooks and Shaw1978). Pollen grains may thus be suitable for morphological analyses, even when DNA is not preserved. In addition, standard DNA extraction protocols for sedimentary samples are designed for processing a specific (rather low) amount of material per sample (in our case, 500 mg), which limits the possibility of using multiple primers and performing multiple PCR replicates, which in turn reduces the probability of detecting rare taxa. Ficetola et al. (Reference Ficetola, Pansu, Bonin, Coissac, Giguet-Covex, De Barba and Gielly2015) suggested that at least eight PCR replicates should be performed if detection probability is not high, but their conclusions were based on empirical studies targeting specific groups of organisms and using specific primers.
The problem of rare species detection in eDNA metabarcoding is complex and still debated. For example, Shirazi et al. (Reference Shirazi, Meyer and Shapiro2021) found that, for soil samples and using generic primers (ITS1 and ITS2), even if performing a very high number of PCR replicates (24), low-abundance taxa were often unique to one or a few replicates. Indeed, one major flaw when using generic primers (such as for past vegetation reconstruction) is that the results might be skewed towards preferential amplification of certain taxa, while others, particularly rare species (Gonzalez et al., Reference Gonzalez, Portillo, Belda-Ferre and Mira2012) remain undetected (see Pedersen et al., Reference Pedersen, Overballe-Petersen, Ermini, Sarkissian, Haile, Hellstrom and Spens2015, and reference therein). In sedaDNA samples, a further complication may derive from the presence of enzymatic inhibitors in sediment extracts, as well as from their removal using reagents included in commercial kits (which can also remove short, damaged DNA fragments). The combination of these two drawbacks can severely affect the retrieval of low-abundance ancient DNA molecules (see Murchie et al., Reference Murchie, Kuch, Duggan, Ledger, Roche, Klunk and Karpinski2021, and reference therein).
Further insights may derive from a quantitative point of view. The palynomorphological analysis highlighted the MB5-D section as much richer than MB3 in terms of total pollen grains (number of pollen counts: 3767 and 256, respectively; MB5-D/MB3 ratio: 14.71). This difference in pollen content is reflected by the DNA-based results: considering the overall results of the different PCR replicates, MB5-D yielded many more sequence reads identified as plant taxa than MB3 (825,593 versus 199,961; MB5-D/MB3 ratio: 4.13). A correspondent quantitative difference emerged in the total number of detected plant taxa: the palynomorphological analysis found 46 taxa in MB5-D and 28 in MB3 (MB5-D/MB3 ratio: 1.64); similarly, our DNA-based results found nine taxa in MB5 and four in MB3 (MB5-D/MB3 ratio: 2.25). Here it is worth noting that the opportunistic experimental design of our study makes a fair comparison difficult for the DNA-based results of these two sections (i.e., different PCR replicates due to different availability of starting material; see Table 1).
In summary, the plant sedaDNA sequences quantitatively and qualitatively seem to reflect endogenous presences rather than potential external origins. None of the plant taxa we found in the speleothem samples was detected in the negative controls with a significant number of reads, despite the large numbers of extraction and PCR blanks included in our experiment. This can be viewed as rule-of-thumb evidence for excluding contamination from modern DNA as the most likely explanation of sedaDNA results.
Issues associated with contaminations from modern DNA and low PCR success rates
We adopted a stringent approach for contamination detection and removal. The rationale is straightforward: a totally contamination-free experimental setting is almost impossible to achieve, even under the most strict conditions (Weyrich et al., Reference Weyrich, Farrer, Eisenhofer, Arriola, Young, Selway, Handsley-Davis, Adler, Breen and Cooper2019), and therefore, if contamination is unavoidable, it has to be carefully monitored. In our experiments, we sequenced a considerably high number of negative controls (42), only four of which were found to be contaminated with significant quantities of plant DNA. The contaminating taxa (cypress, yew, and alder) were not detected in any of the PCR replicates from the speleothems above the adopted 100 reads cut-off threshold. The results for Fagus sylvatica lie in a “gray zone”: adopting a cautionary approach, its detection in four blanks, even in low quantities (<150 reads), prevents us assigning 100% credibility to its presence in the ancient samples. Beyond this single exception, results from speleothem samples and negative controls are not overlapping, indicating that the plant DNA found in the examined speleothem sections is very unlikely derived from contaminations.
Nevertheless, a more critical evaluation of the plant species recorded in speleothems by sedaDNA highlights two points: (1) all of them can be found commonly in modern vegetational assemblages of Alpine areas, including the surroundings of the analytical laboratory; and (2) several of them are characterized by high pollen production and/or persistence (e.g. Pinus, Corylus, Fraxinus, Poaceae, etc.).
We recorded PCR success rates in sample replicates (18/88; 20.5%; see Table 1) that are not much higher than those in the negative controls (extraction and PCR blanks: 4/44; 9.1%). These two critical aspects leave our results open to question. It should be noted, however, that low rates of success and high heterogeneity of results among PCR replicates of the same samples (from zero to high number of reads, for a considered taxon) are not unusual in many sedaDNA studies (Capo et al., Reference Capo, Giguet-Covex, Rouillard, Nota, Heintzman, Vuillemin and Ariztegui2021).
Lastly, when comparing sedaDNA-based data with those from morphological identification, our results again are not unexpected. In this study, the taxa identified by sedaDNA (11) are fewer than those (46) identified by morphological pollen analyses. This is consistent with the outcome of several other sedaDNA studies that were performed using the same barcode gene, trnL (Jørgensen et al., Reference Jørgensen, Haile, Möller, Andreev, Boessenkool, Rasmussen and Kienast2012; Parducci et al., Reference Parducci, Matetovici, Fontana, Bennett, Suyama, Haile, Kjaer, Larsen, Drouzas and Willerslev2013, 2015; Pedersen et al., Reference Pedersen, Ginolhac, Orlando, Olsen, Andersen, Holm, Funder, Willerslev and Kjær2013). For example, Parducci et al. (Reference Parducci, Väliranta, Salonen, Ronkainen, Matetovici, Fontana, Eskola, Sarala and Suyama2015) identified a total of 32 pollen taxa in peat samples dating back at ca. 8500–5900 yr BP, and thus spanning the age of our speleothems, whereas sedaDNA identified only nine taxa. This ratio of 28% is comparable to the 24% ratio measured in our study. It should be noted, however, that the performance of sedaDNA compared with pollen data may vary according to the nature and conditions of the sediment/site and the taxonomic group: opposite examples also exist, with more taxa detected by molecular than morphological analysis, particularly for lake sediments and cold environments (e.g., Clarke et al., Reference Clarke, Alsos, Edwards, Paus, Gielly, Haflidason and Mangerud2020; see Capo et al., Reference Capo, Giguet-Covex, Rouillard, Nota, Heintzman, Vuillemin and Ariztegui2021, for a more detailed discussion). In summary, in the context of sedaDNA and with specific regards to “difficult” mineral deposits, low PCR success rates and limited overlap with results deriving from morphological identification cannot be considered a unique peculiarity of our study.
Limitations of PCR-based amplicon sequencing for speleothem sedaDNA and future perspectives
To our knowledge, only three studies have been published so far on aDNA extracted from speleothems. Dhami et al. (Reference Dhami, Mukherjee and Watkin2018) used an amplicon sequencing approach, but the authors only focused on the identification of microbial communities associated with biomineralization processes. Zepeda Mendoza et al. (Reference Zepeda Mendoza, Lundberg, Ivarsson, Campos, Nylander, Sallstedt and Dalen2016) and Stahlschmidt et al. (Reference Stahlschmidt, Collin, Fernandes, Bar-Oz, Belfer-Cohen, Gao and Jakeli2019) presented the first attempts to exploit speleothems as past biodiversity archives. Both studies adopted a shotgun metagenomics approach. With only six identified genera (including animals and plants), Stahlschmidt et al. (Reference Stahlschmidt, Collin, Fernandes, Bar-Oz, Belfer-Cohen, Gao and Jakeli2019) admitted that their samples were likely characterized by low DNA preservation. Zepeda Mendoza et al. (Reference Zepeda Mendoza, Lundberg, Ivarsson, Campos, Nylander, Sallstedt and Dalen2016) found a large spectrum of taxa, but concluded that the analyzed speleothems could not be used confidently as paleoarchives due to the high degree of DNA damage (perhaps induced by fungal activity) that negatively affected sequence quality and therefore the taxonomic resolution (for plants, results are limited to the class level). The authors therefore suggested the use of target amplicon metabarcoding to provide more reliable and accurate results. In the present study, we followed this direction, but, despite strictly adopting the best practices recommended for sedaDNA (Capo et al., Reference Capo, Giguet-Covex, Rouillard, Nota, Heintzman, Vuillemin and Ariztegui2021), some uncertainties remain.
Some drawbacks are typical of the PCR-based metabarcoding methodology, such as (1) the unavoidable trade-off between “universality” (i.e., extensive taxonomic coverage) and taxonomic resolution in the choice of barcode markers (Hollingsworth et al., Reference Hollingsworth, Graham and Little2011); (2) preferential amplification of certain taxa over other ones, which leads to taxonomic bias (Elbrecht and Leese, Reference Elbrecht and Leese2015; Zinger et al., Reference Zinger, Bonin, Alsos, Bálint, Bik, Boyer and Chariton2019); and (3) the extremely high sensitivity (due to PCR amplification) to contamination by DNA from external sources. As demonstrated in previous work (e.g., Weyrich et al., Reference Weyrich, Farrer, Eisenhofer, Arriola, Young, Selway, Handsley-Davis, Adler, Breen and Cooper2019), contaminating DNA molecules are ubiquitous, even in ultra-clean labs and commercial reagents. The risk of detecting contamination from modern DNA sources is particularly severe when processing ancient eDNA samples (Ficetola et al., Reference Ficetola, Pansu, Bonin, Coissac, Giguet-Covex, De Barba and Gielly2015, Reference Ficetola, Taberlet and Coissac2016), which generally are characterized by low concentrations of damaged DNA and require extreme PCR amplification protocols, such as a high number of cycles: under these conditions, even a few modern intact DNA molecules can be amplified preferentially over a larger amount of ancient damaged molecules (Zinger et al., Reference Zinger, Bonin, Alsos, Bálint, Bik, Boyer and Chariton2019). To partially overcome these drawbacks, we (1) selected a short region of a barcode marker (trnL), which is known for its “universal” detection power (Taberlet et al., Reference Taberlet, Coissac, Pompanon, Gielly, Miquel, Valentini, Vermat, Corthier, Brochmann and Willerslev2007); and (2) extensively monitored contamination, from DNA extraction to PCR and post-metabarcoding analysis. Nevertheless, the outcome was still not optimal.
We think that a significant improvement can be offered by alternative approaches that are not directly based on PCR such as hybridization capture, also called target capture (Shokralla et al., Reference Shokralla, Gibson, King, Baird, Janzen, Hallwachs and Hajibabaei2016; Wilcox et al., Reference Wilcox, Zarn, Piggott, Young, McKelvey and Schwartz2018). In hybridization capture, a library derived from eDNA is hybridized with specific probes. After discarding the libraries with no sign of hybridization, the resulting libraries are massively sequenced. Pioneer attempts to use targeted capture for ancient environmental DNA, such as in Murchie et al. (Reference Murchie, Kuch, Duggan, Ledger, Roche, Klunk and Karpinski2021), which aimed at reconstructing plant and animal communities from permafrost samples dating to the Pleistocene-Holocene transition, yielded promising results with a generally better performance compared with the standard metabarcoding approach. We strongly encourage this direction, also for speleothem and other ancient sedimentary samples. Target capture is particularly well suited for sedaDNA for the following main reasons: (1) it is PCR-free, thus avoiding limitations such as high sensitivity to contaminations from modern DNA, preferential amplification, and heterogeneity among PCR replicates; (2) it offers higher taxonomic resolution due to the possibility to design probes with high sequence homology to the target taxa; and (3) it preserves the typical damage patterns found at or near the ends of aDNA molecules (such as cytosine deamination; Armbrecht et al., Reference Armbrecht, Hallegraeff, Bolch, Woodward and Cooper2021). The latter point is of crucial importance, providing the possibility to authenticate ancient DNA by means of specific bioinformatic tools (see for example MapDamage; Jónsson et al., Reference Jónsson, Ginolhac, Schubert, Johnson and Orlando2013) and to unambiguously identify contaminations from modern DNA. This is not possible for amplicon sequencing given that, during PCR amplification, primers either bind to or exclude the ends of the template molecules where damage is present, leading to the loss of the damage signal.
The aforementioned drawbacks of PCR-based metabarcoding prompt us to suggest that this approach should be avoided when analyzing DNA retrieved from ancient sediments such as speleothems. It can be argued that our case study, with the retrieval from speleothems of a limited number of plant taxa that commonly also are found in modern environments, may represent a particularly unfortunate example. If some particular plant species, currently absent in Alpine latitudes and altitudes, but which might have occurred in (warmer) past paleoenvironments, had been found in the samples, one probably would have been tempted to consider these results reliable. However, if a method can be considered valid only ex post, that is after checking the results, this method is severely biased. We conclude that speleothems could hold promise as biological archives to be studied with molecular methods, provided different approaches than those based on direct PCR amplification of the extracted DNA are used.
Acknowledgments
The authors thank the anonymous reviewers for their valuable and constructive comments on an earlier version of the manuscript.
Supplementary Material
The Supplementary Material for this article can be found at https://doi.org/10.1017/qua.2022.46.
Open access
