A pilot study of eDNA metabarcoding to estimate plant biodiversity by an alpine glacier core (Adamello glacier, North Italy)

Current biodiversity loss is a major concern and thus biodiversity assessment of modern ecosystems is compelling and needs to be contextualized on a longer timescale. High Throughput Sequencing (HTS) is progressively becoming a major source of data on biodiversity time series. In this multi proxy study, we tested, for the first time, the potential of HTS to estimate plant biodiversity archived in the surface layers of a temperate alpine glacier, amplifying the trnL barcode for vascular plants from eDNA of firn samples. A 573 cm long core was drilled by the Adamello glacier and cut into sections; produced samples were analyzed for physical properties, stable isotope ratio, and plant biodiversity by eDNA metabarcoding and conventional light microscopy analysis. Results highlighted the presence of pollen and plant remains within the distinct layers of snow, firn and ice. While stable isotope ratio showed a scarcely informative pattern, DNA metabarcoding described distinct plant species composition among the different samples, with a broad taxonomic representation of the biodiversity of the catchment area and a high-ranking resolution. New knowledge on climate and plant biodiversity changes of large catchment areas can be obtained by this novel approach, relevant for future estimates of climate change effects.


Results
Visual stratigraphy of the ice core. The 573 cm of ADA15 core embodied the surface snow accumulation and the upper part of ice/dense firn of the Pian di Neve, the accumulation basin of Adamello Glacier. Stratigraphically, the upper 212 cm were characterized by single snowfalls showing texture changes (particularly in grain size); the snow layers were clearly visible in the thin sections as variations in gray tones (under transmitted light) (Fig. 1a).
At 212 cm-depth, a more complex snow facies started, containing ice lenses with vertical structures, probably related to water percolation. These facies stopped at 245 cm-depth in a 5 cm-thick ice layer. This layer may be originated by surface melting, in a period with high temperature or insolation, causing percolation of water in the snowpack, and a freezing to ice within the snowpack at lower temperatures.
At a depth from 250 to 429 cm the snow facies were clearly similar to the top part, with some thin ice layers, or lenses, along the core. Some differences could be observed in snow density, mainly due to the snow grain size.
At a 492 cm-depth a transition to the ice facies was visible. After 4-5 cm of dense firn, the bottommost 81 cm of the core were represented by glacier ice, probably superimposed ice, related mostly to the complete percolation and melting of snowpack. The dense firn, in transition to ice, showed vertical and elongated and/or dendritic bubbles that represent the loss of porosity due to the formation of ice crystals. Below this small layer, there were 30 cm of ice with ice crystals of few mm in diameter and few bubbles, consistent with a strong recrystallization of the snow/firn packs after strong melting. In-fact, at 522 cm-depth, a dark layer was present, related to high accumulation of dust, representing probably a previous summer/spring glacier surface (Fig. 1b).
Below the dark level till the bottom of the core, glacier ice was present, with small crystals and more dispersed bubbles. This facies of the glacier seems to be related to the ice metamorphism, transforming snow/firn to glacier ice, even if we cannot exclude melting processes at this level (Fig. 1b).
Stable isotopes. Figure 2 shows the δ 2 H and δ 18 O values of melted ice. The δ 18 O values ranged between around − 17 and − 8‰ and δ 2 H values between around − 140 and − 50‰.
By evaluating the data along the ice core depth, it was possible to observe some general trends and several oscillations. Both δ 2 H and δ 18  Microscopical pollen analyses. Palynological analyses on three random samples confirmed the presence of pollen grains within the ice core ( Fig. 3; Table 1). In addition to pollen, various other debris were present in the samples (Fig. 3 Table S1). Samples were subjected to rarefaction (without replacement) at 16,948 reads per sample, corresponding to the minimum number of reads per sample (sample PI2). Reads were clustered in a total of 2339 OTUs, with a maximum of 1047 OTUs per sample (corresponding to the sample PI5). Sequencing results are detailed in the Supplementary Table S1.
The obtained dataset was the basis for the following taxonomic assignment.
Taxonomic assignment of environmental DNA sequences. The data analysis of OTUs was performed against the local reference database at a 99% similarity, resulting in 2164 OTUs (92.5%) taxonomically assigned to a single plant taxa. Overall, Fabaceae, Pinaceae and Poaceae were the most represented plant families, representing 69% of the total reads with a taxonomic assignment in the reference database. The taxonomic distribution of the OTUs spanned the whole spermatophyta clade from gymnosperms to angiosperms (Fig. 4), representing a total of 10 orders. In gymnosperms only the Pinaceae family was represented, while in angiosperms a higher number of families known to encompass species with anemophilous pollination was identified, spanning both monocotyledons and dicotyledons. In dicotyledons, both rosids and fabids clades were represented. While Pinaceae and Poaceae were identified also by light microscopy analyses, Fabaceae were identified exclusively by DNA metabarcoding. The two families with the highest numbers of reads, Pinaceae and Fabaceae, encompassed alone > 60% of the reads that could be identified. While in general the number of reads was often higher in families with higher numbers of taxa (see, e.g. Poaceae, Fabaceae and Pinaceae), the Caprifoliaceae were represented by only one genus (Lonicera; Fig. 4; Supplementary Table S2). In the case of sequences identified as Ficus carica, a curated BlastN search was carried out to the NCBI nucleotide database, resulting in the high-confidence identification of the sequences as Broussonetia papyrifera (L.) Vent., an invasive species of Asian origin.
The overall diversity is shown in the doughnut chart of Fig. 5, where each sample is symbolized by a circle, with the arch width proportional to the frequency of plant families. Plant biodiversity is clearly changing from one sample to another, supporting the estimated different flowering periods and years, as suggested by conventional light microscopy pollen analyses too.
When analyzed by Chao1 and Shannon diversity indices, the alpha diversity estimated for the samples resulted in a large range of values. Both indices assigned a very high alpha diversity to samples PI3 and PI5 and low alpha diversity for sample PI4 (Fig. 6).   Table S2) described 36 different plant taxa, with different taxonomic resolution. eDNA metabarcoding, associated with high throughput sequencing, assigned taxonomy mostly to the species level (21; 58.3%). In only one case the taxonomy was assigned to more than one species of the same genus (2.8%); in three cases to the genus level (8.3%) and in other three to more than one genus in the same family (8.3%). Assignments at the family level covered 8 cases (22.2%).

Discussion
In this study we tested whether metabarcoding of eDNA using the P6 loop of the plastid DNA trnL (UAA) intron 14 could be used as a proxy of vascular plant diversity when analyzing surface layers of the largest Italian glacier, the Adamello Glacier. The high sequence diversity of the trnL p6 loop region 14 allowed us to use a stringent cutoff of 99% sequence identity for the taxonomic assignment of the OTUs, higher than the 97% adopted by similar studies (e.g. 15,52 ). Despite the reduced length of the barcode marker sequence (138-152 bp, depending on the taxa), such approach resulted in a higher accuracy of the taxonomic assignment, as compared with that achievable by conventional palynology 3 . As expected, the two methods did not provide fully overlapping taxonomic representations, as indicated by the method-specific identification of some of the taxa. Palynological analyses by conventional light microscopy, however, were performed on a low number of samples, since they aimed mainly at confirming the presence of pollen rather than to compare the two approaches, as this has been extensively done in previous studies (reviewed in 53 ). Additionally, as shown by the conventional light microscopy analysis, glacier ice trapped, besides pollen grains, also plant remains, likely blown by wind during storms 54 , which can further contribute to assess plant species composition of the surrounding area, and to consider palynological and molecular analyses as complementary tools for biodiversity assessment 53 . Also the taphonomic processes acting on the pollen grains in the ice archive can contribute, at least in part, to the observed differences between the taxa identified by the two methods, although taphonomy of eDNA from glaciers remains largely to be investigated. Compared to eDNA from lake sediments, where extracellular DNA has been suggested to constitute the major component 18 , extracellular DNA in the ice samples of a temperate glacier as Adamello is expected to be a minor fraction for the following reasons: (1) repeated cycles of freezing and thawing can mechanically degrade exposed DNA through cryolysis 55 ; (2) extracellular DNA in the ice cannot complex with soil particles, which have been demonstrated to protect it from nuclease action 17,56 ; (3) If dissolved in melting water, despite the low temperature, the DNA is exposed to hydrolysis, due to the low ionic strength and to physical degradation due to UV exposure 18,57,58 ; (4) once dissolved in melting water extracellular DNA can be diluted below the detection limit or completely removed through percolation. As the visual stratigraphy and the stable isotope analyses we carried out suggest significant processes of meltwater percolation and refreezing of the ice core, we suggest that the majority of DNA analyzed by metabarcoding is of intracellular origin, both from pollen grain and plant debris. The lack of Cupressaceae and Taxaceae sequences from the metabarcoding record is in line with this scenario, as for these families, the exine layer of pollen walls is known to be very thin and easily broken 59,60 , causing the release of the pollen protoplast due to pollen bursting 61 . Once in the low ionic strength meltwater, the protoplasts can undergo rupture by hypo-osmotic stress and release their DNA extracellularly. Analogous taphonomic processes may differentially affect DNA preservation as well as transport processes in different species, possibly causing a partially biased representation of the plant biodiversity of the palynological catchment area.
In this study, eDNA was analyzed at six depth ranges, for a total length of about 6 m. The marked variations among ice core fractions detected by both the palynological and the isotope ratio analyses suggested that the investigated ice layers correspond to multiple seasons, possibly in the range between 4 to 6 years. The progressive reduction of signal in the stable isotope ratio data suggested a possible modification in the solid precipitation after the deposition, by processes of meltwater percolation and refreezing, which can reduce seasonal isotopic signals, induce isotopic enrichment and introduce time gaps 62 . This could be the case for the sampled ice core described by these analyses; deeper drilling in additional locations will allow to ascertain whether melt-affected ice is a general feature of the Adamello glacier or if it is a local peculiarity at the drilling site. Due to these processes, chemical signals are generally susceptible to various modifications which can affect the quality and quantity of information when analyzed on temperate glaciers. Palynological analyses, on the contrary, demonstrated to provide valuable information even when in temperate glaciers, where low latitude and climatic conditions do not guarantee freezing temperatures all year long 31 , due to the stability of these biological particles in the ice layers.
In general, the taxonomic representation resulting from eDNA metabarcoding was surprisingly broad, reflecting many plant species from the families of Betulaceae, Fagaceae, Oleaceae, Pinaceae, Poaceae and Urticaceae, typical of the alpine habitats at different elevations within the palynological catchment area of the Adamello glacier 63 . In our study, species accumulation curve reached saturation (see Supplementary Fig. S1), indicating the exhaustive sampling of almost all species present in DNA samples. Moreover, taxonomic assignment after eDNA metabarcoding was performed against a custom-made reference database 16 , ensuring data quality by a robust set of reference DNA sequences, as suggested by 64 . The only OTU identified by eDNA metabarcoding which could not directly be explained by the custom reference database used in this study, was the OTU initially identified as Ficus carica, a mediterranean species, which does not spontaneously grow in the immediate proximity of the glacier. Manual curation of the homology search, however, easily fixed the representativeness problem and identified the taxon as Broussonetia papyrifera, an invasive species, occasionally naturalized, and present in several Italian regions, including Trentino, Alto-Adige and Lombardy, which surround the Adamello glacier 65 . Recently, B. papyrifera has been recognized as a source of pollinoses in Asia 66 , and it is currently object of monitoring by PollNet, the Italian network of aerobiological monitoring centres (http://www.snpam bient e.it/2019/05/10/primo -studi o-in-itali a-sulla -distr ibuzi one-del-polli ne-di-brous sonet ia-papyr ifera / downloaded at February 28, 2020 www.nature.com/scientificreports/ detection of B. papyrifera pollen on the Adamello glacier indicates the high potential of glacier ice as a proxy of vascular plant diversity of the catchment area, including even alien/invasive species. Nowadays the integration of the DNA metabarcoding and the morphological approaches represent the best way to investigate the major part of the ecological sites 67 . Despite the accuracy of DNA metabarcoding analysis is strictly related to the reference database used and its quantification capability is still under study 68 , the advantages of the metabarcoding approach (overcoming the morphological approach limits, low cost and low time consumption) make it a fundamental method for ecological survey 6,53 .
In conclusion, we demonstrated for the first time that a biological signal is archived in the accumulation basin of the Adamello glacier and that eDNA metabarcoding is a feasible technique when applied to eDNA of firn and ice of Alpine glaciers. Moreover, eDNA metabarcoding gives a high-ranking taxonomic resolution to pollen signal, capable of disclosing plant biodiversity of the palynological catchment area, estimated by the traces of pollen and plant remains archived in this Alpine temperate glacier. The potential of sampling glacier DNA with eDNA metabarcoding will need to be further investigated in other glaciers and to different depths, to fully exploit the wealth of information still contained in these important archives of natural plant biodiversity. The relevance of this information on heavily endangered sites could improve existing knowledge on plant biodiversity changes during the past in large catchment areas and on the link of these changes to distinct climatic conditions, leading to results that could be applied on estimates of climate change effects on the future vegetational composition of the ecosystem.

Methods
Study site. The drilling site (coordinates UTMWGS84, E 617616-N 5111704) is located on Pian di Neve, the vast plateau that represents the accumulation basin of the Adamello Glacier which is situated between Brescia and Trento provinces in the Adamello-Presanella group (Fig. 7). Recent seismic geophysical surveys estimated the maximum thickness of the glacier (270 m; 69 ) and the main directions of the glacial flows.
Adamello Glacier has been classified as a unique glacial unit 70 and comprises the six units in which the glacier was divided in the Inventory of the Italian Glaciological Committee 71 : Miller Superiore (CGI n.600), Corno Salarno (CGI n.603), Salarno (CGI n. 604), Adamello or Pian di Neve (CGI n.608), Adamè (CGI n.609) and Mandrone (CGI n. 639). Considered in its entirety, Adamello Glacier can be classified as a Scandinavian summit glacier 72,73 . More precisely, it should be classified as a plateau glacier with radial tongues.
Pian di Neve is a large accumulation area from which the different tongues originate, the main of which is channeled towards the Val Genova (Mandrone Glacier). Pian di Neve plateau receives the contribution of other www.nature.com/scientificreports/ glacial flows originating from both the north (Corno Bianco) and the south-west (Corno di Miller, Corno di Salarno and Cornetto di Salarno). Adamello Glacier is fed by direct snowfall. In recent years the snow limit has often been placed at very high altitudes and the entire central glacier sector, including Pian di Neve, has been found in ablation conditions. The accumulation areas are reducing their size and are now concentrated in correspondence with the highest cirques, both in the eastern and western glacier sectors. Moreover, the open morphology reduces the shading effect due to the rocky ridges that delimit it.
Adamello Glacier is currently the largest glacier in the Italian Alps, with an area of 15.5 km 2 in 2015 (data downloadable at https ://webgi s.provi ncia.tn.it, update 2015 by Casarotto, C. and Trenti, A., accessed 28 April 2020). The total area of the glacier exceeded 30 km 2 in 1864 74  Ice core processing and visual stratigraphy. A 573 cm long core was drilled through the Adamello glacier (Trentino, Italy) in March 2015, by a mechanical drill system (SIPRE type, USA), with a core diameter of 100 mm and length of approximately 60 cm for each section. The sampled core was cut at Eurocold Lab, University of Milano Bicocca, in 5 cm sections for the stable isotope ratio analysis, and 10 cm sections for environmental DNA (eDNA) analyses. Samples were stored at − 20 °C. Thin sections, 1 mm thick, were cut from the ice core and placed between two crossed polarization filters. Size and orientation of the ice crystals were evaluated by the different color of the crystals, due to the polarization of the light.
Stable isotope ratio. Samples were thawed at room temperature in hermetically sealed tubes. As described in previous papers 78,79 , the ratio 18 O/ 16 O of water was analysed using an isotope ratio mass spectrometer (Isoprime, Manchester, UK) interfaced with an on-line automatic system that allows CO 2 /H 2 O equilibration (Multiflow, Isoprime, Manchester, UK). For the analysis of ratio 2 H/ 1 H of water, 200 μl of water together with Hokko bead platinum catalyst (Isoprime) were introduced into a reaction vessel, which was then attached to the online automatic equilibration system and filled with He containing 10% of H 2 . The ratios 18  Light microscopy pollen analysis. In order to verify the presence of pollen grains within the ice core, random samples, about 35 ml each were concentrated and processed with acetolysis 80 . Afterwards, pollen was mounted on slides, coloured with basic fuchsin and analysed by light microscopy (400 × magnification) for the identification of pollen grains by means of standard identification keys 81 , and comparison to a local pollen collection (Fondazione Mach pollen reference collection) and to pollen atlas (http://polli ni.iasma .it/polim age; 82 ).
Environmental DNA sample processing and extraction. Ten cm-long samples were thawed in sterile glass containers in a refrigerator for 24 h and subsampled for eDNA, as to be equally represented in the resulting composite samples. The six composite samples covered the entire length of the ice core, representing the following relative depth from the top: 0-172 cm (PI1); 173-275 cm (PI2); 276-449 cm (PI3); 450-489 cm (PI4); 490-509 cm (PI5); 510-573 cm (PI6).
A 36 ml aliquot of each composite sample, together with a negative control, was individually concentrated (Concentrator Plus, Eppendorf, Hamburg, Germany) to a final volume of 100 µl (6 to 7 h).
DNA extraction protocols were applied under a high performance Class II biosafety cabinet in a non-invasive DNA laboratory, free from contaminant small DNAs. Details on the procedures and controls used to prevent and identify contaminations are provided in Supplementary Methods. An automated high throughput (HTP), magnetic-bead based isolation of genomic DNA was performed by NucleoMag Plant, (Macherey-Nagel, Düren, Germany) and KingFisher Flex platform (Thermo Scientific, MA USA). Cell lysis was promoted by mechanical disruption (Mixer Mill MM200, Retsch GmbH, Germany). In order to control contamination during extraction, samples were loaded in plates, following a design where empty wells surrounded each sample (see Supplementary methods). The final elution volume for all samples was 100 µl. www.nature.com/scientificreports/ Each sample (negative controls and blank included) was amplified by a 1 st round of PCR using 50 µl reaction with 1 µM of each primer. More in detail, 10 µl of 5 × GO Taq Flexi Buffer and 2 µl forward and reverse primers, were used in combination with 5 µl of template DNA (5 ng/µl), 6 µl MgCl 2 25 mM, 0.25 µl dNTP 40 mM and 0.25 µl GO Taq Flexi [5U/µl]. PCR reactions were executed in triplicates by a Veriti 96 well thermal cycler (Applied Biosystems, Foster City, CA, USA) at the following cycling conditions: initial denaturation step at 95 °C for 2 min (one cycle); 40 cycles at 95 °C for 15 s, 52 °C for 15 s, 72 °C for 30 s; final extension step at 72 °C for 5 min (1 cycle).
The amplification products in triplicate were pooled, subsequently the amplicons library construction and the high throughput sequencing were performed at the FEM Sequencing Platform using an Illumina MiSeq (MiSeq Control Software 2.5.0.5 and Real-Time Analysis software 1.18.54.0).

Bioinformatics analysis.
Raw reads were processed with the MICrobial Community Analysis software MICCA v. 1.7.2 83 by setting the OTUs clustering at the 99% sequence similarity level and 75% query coverage, and, subsequently, OTU classification was performed against the BioAir reference database 16 . Chao1 and Shannon indices were estimated to calculate alpha diversity in terms of OTUs richness and diversity. Taxonomic tree for the OTUs was reconstructed using the NCBI taxonomy and visualized in

Data availability
Raw sequences were deposited in the European Nucleotide Archive ENA-EMBL, under the accession number PRJEB28043. www.nature.com/scientificreports/