Prof Jean-Michel Claverie
email: michel.claverie@univ-amu.fr

Keywords: metagenomics, Kamchatka, Yakutia, Siberia, permafrost thaw, antibiotic resistance genes

Abstract

In the context of global warming, the melting of arctic permafrost raises the threat of a reemergence of microorganisms some of which were shown to remain viable in ancient frozen soils for up to half a million years. In order to evaluate this risk, it is of interest to acquire a better knowledge of the composition of the microbial communities found in this understudied environment.

Here we present a metagenomics analysis of 12 soil samples from Russian Arctic and subarctic pristine areas: Chukotka, Yakutia, and Kamchatka, including 9 permafrost samples collected at various depths. These large datasets (9.2 10 11 total bp) were assembled (525,313 contigs > 5kb), their encoded protein contents predicted, then used to perform taxonomical assignments of bacterial, archaeal, and eukaryotic organisms, as well as DNA viruses.

The various samples exhibited variable DNA contents and highly diverse taxonomic profiles

showing no obvious relationship with their locations, depths or deposit ages. Bacteria represented the largely dominant DNA fraction (95%) in all samples, followed by archaea (3.2%), surprisingly little eukaryotes (0.5%), and viruses (0.4%).

Although no common taxonomic pattern was identified, the samples shared unexpected high frequencies of β-lactamase genes, almost 0.9 copy/bacterial genome. In addition of known environmental threats, the particularly intense warming of the Arctic might thus enhance the spread of bacterial antibiotic resistances, today's major challenge in public health. β-lactamases

were also observed at high frequency in other types of soils, suggesting their general role in the regulation of bacterial populations.

O R I G

I N A L U N E D I T E D

Introduction

Two of the front-page societal concerns are global warming and the increasing frequency of emerging infectious diseases, eventually turning into pandemics [START_REF] Morens | Pandemic COVID-19 Joins History's Pandemic Legion[END_REF][START_REF] Zhou | A pneumonia outbreak associated with a new coronavirus of probable bat origin[END_REF]. It turns out that both predicaments are partly linked. Through its influence on the ecology of living organisms vectoring infectious pathogens, the warming climate contributes to spreading endemic diseases from tropical regions into temperate areas of the globe, such as Western Europe [START_REF] Hertig | Distribution of Anopheles vectors and potential malaria transmission stability in Europe and the Mediterranean area under future climate change[END_REF][START_REF] Fuehrer | Monitoring of alien mosquitoes in Western Austria (Tyrol, Austria[END_REF][START_REF] Liu | Reviewing estimates of the basic reproduction number for dengue, Zika and chikungunya across global climate zones[END_REF]). This phenomenon is further amplified by the encroachment into pristine areas (e.g. tropical forests) where expanding human activities generate new risks through imponderable contacts with (mostly uncharted) microbial environments and their associated wild life hosts [START_REF] Bradley | Urbanization and the ecology of wildlife diseases[END_REF][START_REF] Keita | Looking in apes as a source of human pathogens[END_REF][START_REF] Plowright | Pathways to zoonotic spillover[END_REF][START_REF] Valentine | Sylvatic cycles of arboviruses in non-human primates[END_REF].

While such dangers are mostly pointed out as coming from the South, more recent concerns have been raised that new plagues could also come from the Arctic, through the release of infectious agents until now trapped in perennial frozen soils (i.e. permafrost) up to 1.5 km deep and 2-3 million years old [START_REF] Revich | Thawing of permafrost may disturb historic cattle burial grounds in East Siberia[END_REF][START_REF] Revich | Climate change and zoonotic infections in the Russian Arctic[END_REF][START_REF] Parkinson | Climate change and infectious diseases in the Arctic: establishment of a circumpolar working group[END_REF][START_REF] Huber | Symposium report: emerging threats for human health -impact of socioeconomic and climate change on zoonotic diseases in the Republic of Sakha (Yakutia), Russia[END_REF].

Climate warming is particularly noticeable in the Arctic where average temperatures are increasing twice as fast as in temperate regions [START_REF] Cohen | Recent Arctic amplification and extreme mid-latitude weather[END_REF]. One of the most visible consequences is the widespread thawing of permafrost at increasing depths [START_REF] Biskaborn | Permafrost is warming at a global scale[END_REF][START_REF] Turetsky | Permafrost collapse is accelerating carbon release[END_REF]) and the rapid erosion of permafrost bluffs [START_REF] Fuchs | Rapid Fluvio-Thermal Erosion of a Yedoma Permafrost Cliff in the Lena River Delta[END_REF], a phenomenon most visible in Siberia where deep continuous permafrost underlays most of the North Eastern territories.

The "ice age bug" concern has been periodically brought back to the public attention. For instance, when an exceptionally hot summer triggered local outbreaks of anthrax on Yamal Peninsula, Northwest Siberia, in 2016, as a deeper than usual summer season thaw of soils above the permafrost layer (i.e., the "active layer") exhumed infectious B. anthracis endospores buried in the frozen ground for 75 years [START_REF] Timofeev | Insights from Bacillus anthracis strains isolated from permafrost in the tundra zone of Russia[END_REF]. Historically frequent outbreaks of anthrax killed 1.5 million reindeer in Russian North between 1897 and 1925 and human cases of the disease have occurred in thousands of settlements across the Russian North [START_REF] Hueffer | Factors Contributing to Anthrax Outbreaks in the Circumpolar North[END_REF].

The capacity of deeper (hence much older) permafrost layers to preserve the integrity of "live" (i.e.

infectious) viral particles, was demonstrated by the isolation of two previously unknown giant DNA viruses: Pithovirus [START_REF] Legendre | Thirty-thousand-year-old distant relative of giant icosahedral DNA viruses with a pandoravirus morphology[END_REF] and Mollivirus [START_REF] Legendre | In-depth study of Mollivirus sibericum, a new 30,000-y-old giant virus infecting Acanthamoeba[END_REF]. Those 1A) radiocarbon dated to 35,000 years before present (BP) [START_REF] Legendre | Thirty-thousand-year-old distant relative of giant icosahedral DNA viruses with a pandoravirus morphology[END_REF]. This first demonstration that eukaryotic DNA viruses could still be infectious after staying dormant since the late Pleistocene came after many other studies [START_REF] Vishnivetskaya | Bacterial community in ancient Siberian permafrost as characterized by culture and culture-independent methods[END_REF][START_REF] Hinsa-Leasure | Characterization of a bacterial community from a Northeast Siberian seacoast permafrost sample[END_REF][START_REF] Graham | Microbes in thawing permafrost: the unknown variable in the climate change equation[END_REF]) going back to 1911 (reviewed by [START_REF] Gilichinsky | Microbial life in permafrost: A historical review[END_REF]). It was shown that ancient Bacteria could be revived from even much older permafrost layers, up to half a million years [START_REF] Johnson | Ancient bacteria show evidence of DNA repair[END_REF].

Contamination by modern bacteria have been suspected in the case of older isolates [START_REF] Willerslev | Isolation of nucleic acids and cultures from fossil ice and permafrost[END_REF]. A flower plant was also revived from 30,000 year old frozen tissues [START_REF] Yashina | Regeneration of whole fertile plants from 30,000-yold fruit tissue buried in Siberian permafrost[END_REF], as well as various fungi (from Antarctica) [START_REF] Kochkina | Ancient fungi in Antarctic permafrost environments[END_REF]) and amoebal protozoans [START_REF] Shmakova | Viable Species of Flamella (Amoebozoa: Variosea) Isolated from Ancient Arctic Permafrost Sediments[END_REF][START_REF] Malavin | Isolates from ancient permafrost help to elucidate species boundaries in Acanthamoeba castellanii complex (Amoebozoa: Discosea)[END_REF]. There is now little doubt that permafrost is home to a large diversity of ancient microorganisms that can potentially be revived upon thawing. Many of them, in particular among viruses, are unknown. If detected, the follow-up questions are: Are those microorganisms a threat for today's society?

Could some of them represent an infectious hazard for humans, animals or plants? Can we identify microbes the metabolisms of which would accelerate the emission of greenhouse gases (e.g. methane, carbon dioxide) upon thawing?

Beyond the direct infectious risk represented by viable microorganisms released from thawing ancient permafrost layers, the persistence of a multitude of genomic DNA fragments from all the organisms that were present in the soil at various times is also to be considered. Although it originated from dead (or even extinct) organisms, the permafrost DNA content constitutes a historical library of genetic resources, from which useful elements (genes) can be reintroduced into modern organisms through bacterial transformation and/or lateral gene transfer in protozoa or viruses. Of particular interest are DNA fragments coding for bacterial virulence factors or antibiotic resistance genes the origin of which predates by far human history [START_REF] Finley | The scourge of antibiotic resistance: the important role of the environment[END_REF][START_REF] Perron | Functional Characterization of Bacteria Isolated from Ancient Arctic Soil Exposes Diverse Resistance Mechanisms to Modern Antibiotics[END_REF][START_REF] Potron | Emerging broad-spectrum resistance in Pseudomonas aeruginosa and Acinetobacter baumannii: Mechanisms and epidemiology[END_REF][START_REF] Wright | Environmental and clinical antibiotic resistomes, same only different[END_REF]).

To address these questions, a first step is to evaluate the DNA in frozen soils. The technique of choice (which avoids the risk of reactivating unknown dangerous pathogens) for this task is metagenomics. In contrast with targeted approaches, such as metabarcoding (i.e. PCR-amplified subset of environmental DNA), it offers the possibility of unexpected discoveries in all domains of life and provides an estimate of the relative abundances of the various organisms, although at the cost of a lower sensitivity. Finally, the metagenomic analysis of DNA obtained from up to sixty grams of frozen soil samples (Table S1) is much less sensitive to eventual contaminations by

O R I G I N A L U N E D I T E D M A N U S C R I P T
minuscule amounts of ambient DNA or contemporary microorganisms and more suitable for siteto-site comparisons than culture-based approaches.

Here we report the analysis of 9 sub-surface Siberian permafrost samples, and of 3 modern surface soils from pristine area in Kamchatka. All the permafrost samples were taken from carbon-rich, icerich frozen soils called Yedoma deposits [START_REF] Schirrmeister | Yedoma: Late Pleistocene ice-rich syngenetic permafrost of Beringia[END_REF][START_REF] Strauss | Deep Yedoma permafrost: A synthesis of depositional characteristics and carbon vulnerability[END_REF] found in vast regions of northeast Siberia (Lena and Kolyma river basins in Yakutia)(Fig. S1) and known for wellpreserved Late Pleistocene megafauna remains (mammoths, wholly rhinoceros). Yedoma permafrost is of special interest as it is prone to rapid thaw processes such as thermokarst and thermo-erosion [START_REF] Nitzbon | Fast response of cold ice-rich permafrost in northeast Siberia to a warming climate[END_REF] releasing not only carbon as greenhouse gases [START_REF] Schneider Von Deimling | Observation-based modelling of permafrost carbon fluxes with accounting for deep carbon deposits and thermokarst activity[END_REF][START_REF] Turetsky | Carbon release through abrupt permafrost thaw[END_REF] but also its formerly freeze-locked microbial content, in which we detected a high abundance of β-lactamases genes carried in the genomes (or plasmids) of a large phyletic diversity of bacteria.

Methods

Sample collection and preparation

We sampled Yedoma and thermokarst sediments from Central Yakutia from the Yukechi Alas landscape about 50 km southeast of Yakutsk (Table 1A). The Yukechi landscape is characterized by Yedoma uplands and drained as well as extant lake basins, so-called Alases, indicating active permafrost degradation processes [START_REF] Fedorov | Observations of surface dynamics with thermokarst initiation, Yukechi site, Central Yakutia[END_REF]. Field work took place in March 2015 during a joint Russian-German drilling expedition. We used a Russian URB2-47 drilling rig with a steel core barrel without a core catching system. Core barrels had different diameters, starting with 15.7 cm for the upper core segments and narrowing down to 8 cm for lower core segments, a setup preventing to get stuck in the deep borehole. To avoid contaminations, coring was conducted as dry drilling and no drilling fluid was used. The retrieved core segments were pushed out from the core barrel with air pressure only. The core segments were visually described and immediately stored in plastic bags or core foil and sealed tightly with duct tape. The wrapped and sealed cores were put into opaque and thermally insulated hard case expanded polypropylene (EPP) boxes (thermoboxes) and stored outside, where they were kept frozen (or refroze slowly if the sediment material was originally unfrozen) at outside temperature around -10 to -20°C. During the entire transportation, the core material was kept in the dark and frozen in the thermo-boxes.

Temperatures below freezing during the entire storage and transport chain were recorded with a temperature data logger in one of these boxes from the drilling site to the laboratories in Potsdam, Germany.

In total, three long (all ~20m below surface/lake ice) permafrost sediment and unfrozen thaw bulb 1A, Fig. S1). Biogeochemical characterization [START_REF] Windirsch | Organic carbon characteristics in ice-rich permafrost in alas and Yedoma deposits, central Yakutia, Siberia[END_REF][START_REF] Ulrich | Geochemistry and Weathering Indices of Yedoma and Alas Deposits beneath Thermokarst Lakes in Central Yakutia[END_REF] as well as greenhouse gas production [START_REF] Jongejans | Greenhouse gas production and lipid biomarker distribution in Yedoma and Alas thermokarst lake sediments in Eastern Siberia[END_REF]) are published. In Potsdam, the frozen cores were handled in a freezer laboratory at -10°C.

They were split lengthwise using a band saw and were subsequently subsampled. After cleaning the splitting surfaces by scratching with cleaned knifes, the core stratigraphies were described visually and photographed. The sampling core halves were cut into subsamples and divided for further analysis, while the archive core halves were wrapped and sealed again for potential future additional analysis. For the samples used in this study, we took the inner part of the cores only to avoid any other sources of contamination on the outer core margins. After extraction, the subsamples were send frozen with dry ice to Marseille, France, for further analyses.

The Stanchikovsky Yar sample (Sample B) from which two live viruses were previously isolated [START_REF] Legendre | Thirty-thousand-year-old distant relative of giant icosahedral DNA viruses with a pandoravirus morphology[END_REF][START_REF] Legendre | In-depth study of Mollivirus sibericum, a new 30,000-y-old giant virus infecting Acanthamoeba[END_REF], and three modern surface soils from pristine cold regions (Kamchatka), were also analyzed for comparison (Table 1B). Fig. S1 indicates the locations of these sampling sites.

DNA extraction

Samples were processed in BSL2 conditions using either the Miniprep (0.25 g of starting material)

or Maxiprep (20 g of starting material) version of the the Dneasy PowerSoil Kit from Qiagen.

Extractions were repeated until about 1 μg (or more) of raw DNA was obtained for each sample (Table S1). Samples B, C, D, E were processed using the miniprep kit, following the manufacturer protocol and by adding DTT to the C1 buffer at a 10 mM final concentration. All other samples (K, L, M, N, O, P, Q, R) were processed using the Maxiprep kit, including the addition of DTT to the C1 buffer at a 10mM final. Samples were ground for 20 seconds in an MP FastPrep homogenizer at a speed of 4 m/s then incubated for 30 minutes at 65 ° C and finally ground again for 20 seconds at 4 m/s. After elution in 5 ml of elution buffer the extracted raw DNA was concentrated on a silica column from the Monarch Genomic DNA purification Kit from NEB. The DNA contents were then quantified using Qu-bit (Thermo Fisher Scientific) and NanoDrop (Thermo Fisher Scientific) assays (Table S1) prior sequencing.

DNA Sequencing, assembly, and annotation

Each metagenomic sample was sequenced with DNA-seq paired-end protocol on Illumina HiSeq4000 platform at the French National Sequencing Center "Genoscope" (https://jacob.cea.fr) producing datasets of 2x150bp read length, except for sample B (2x100bp read length) (Table S1).

Sequencing libraries were prepared as previously described in [START_REF] Alberti | Viral to metazoan marine plankton nucleotide sequences from the Tara Oceans expedition[END_REF] using an 'on

O R I G I N A L U N E D I T E D M A N U S C R I P T
beads' protocol. Briefly, when available, 250 ng of genomic DNA were sonicated using the E210

Covaris instrument (Covaris, Inc., USA) and the NEBNext DNA Modules Products (New England Biolabs, MA, USA) were used for end-repair, 3`-adenylation and ligation of NextFlex DNA barcodes (Bioo Scientific Corporation). After two consecutive 1x AMPure XP clean ups, the ligated product were amplified by 12 cycles PCR using Kapa Hifi Hotstart NGS library Amplification kit (Kapa Biosystems, Wilmington, MA), followed by 0.6x AMPure XP purification. Finally, libraries were quantified by qPCR (MxPro, Agilent Technologies, USA) and libraries profiles were evaluated using an Agilent 2100 Bioanalyzer (Agilent Technologies, USA).

Reads were quality checked with FASTQC (v0.11.7) (Andrews 2014). Identified contaminants were removed and remaining reads were trimmed on the 3' end using 30 as quality threshold with BBTools (v38.07) (BBMap). Assemblies of filtered datasets were performed using MEGAHIT (v1.1.3) [START_REF] Li | MEGAHIT: an ultra-fast single-node solution for large and complex metagenomics assembly via succinct de Bruijn graph[END_REF] with the following options: --k-list 33,55,77,99,127 --min-contig-len 1000. All filtered reads were then mapped to the generated scaffolds using bowtie2 [START_REF] Langmead | Fast gapped-read alignment with Bowtie 2[END_REF] with the --very-sensitive option. The mapping of the contigs was used to calculate the mean coverage of each contig. To ensure a reliable taxonomic assignation, only scaffolds longer than 5kb

were considered and submitted to metaGeneMark (v3.38) [START_REF] Zhu | Ab initio gene identification in metagenomic sequences[END_REF] with the default metagenome-style model. The predicted genes were then compared to the Refseq database [START_REF] Li | RefSeq: expanding the Prokaryotic Genome Annotation Pipeline reach with protein family model curation[END_REF] using the Diamond version of blastP (v0.9.31) [START_REF] Buchfink | Fast and sensitive protein alignment using DIAMOND[END_REF] with an E-value < 10 -5 , identity percentage threshold of 35 and retaining only the best hits.

The taxonomy of each scaffold was then inferred with a custom-made script applying a Last

Common Ancestor (LCA) method. At each rank, recursively, the taxonomy was conserved only if over half of the annotated genes presented the same keyword. Unclassified and ambiguous contigs were screened against Refseq with Diamond blastX (--sensitive option and E-value < 10 -5 ), in order to detect exons. Non-overlapping best hits were used. The above taxonomical protocol was applied and only eukaryotic scaffolds were conserved. An alternative taxonomical annotation at the single read level was computed by Kraken (v2) (Wood, Lu and Langmead 2019) using its standard database. This annotation was found to be much less comprehensive and was not used further.

The diversity and sequencing effort of each sample was estimated by Nonpareil (v3.304) [START_REF] Rodriguez-R Lm | Nonpareil 3: Fast Estimation of Metagenomic Coverage and Sequence Diversity[END_REF] on forward reads over 35 bp using the kmer algorithm.

Contig taxonomy and coverage data was used for beta-diversity estimation together with the kmer frequency of reads. To do so, 15 bp kmer frequency was retrieved for forward and reverse reads by Gerbil (v1.12) [START_REF] Erbert | Gerbil: a fast and memory-efficient k-mer counter with GPU-support[END_REF]). An in-house script was used to merge the produced files and compute the cosine distance between samples.

O R I G

I N A L U N E D I T E D M A N U S C R I P T
The Pfam domains (Nov. 2021 database version, [START_REF] Mistry | Pfam: The protein families database in 2021[END_REF]) of all ORFs were searched with Interproscan (v.5.39-77.0, Jones et al. 2014). The non-overlapping matches with E-value < 10 -5

were retrieved with an in-house script.

β-lactamase gene detection and analyses

ORFs with a best Diamond blastP match to a β-lactamase and β-lactamase-related proteins in the NCBI Refseq database were retrieved and analysed (Table S2). Functional domains were inferred by Interproscan (v5) against Pfam, ProSitePatterns, and by Batch-CDsearch (Marchler-Bauer and Bryant 2004) against the COG database [START_REF] Galperin | COG database update: focus on microbial diversity, model organisms, and widespread pathogens[END_REF]) (E-value < 10 -5 ) (Table S3). ORFs at least 100 amino acids in length exhibiting a β-lactamase domain according to Pfam and COG (except for class D solely defined by COG domain) were kept for further analyses. We attempted to identify the β-lactamase-encoding contigs as potential plasmids using PlasClass [START_REF] Pellow | PlasClass improves plasmid sequence classification[END_REF] and PlasFlow [START_REF] Krawczyk | PlasFlow: predicting plasmid sequences in metagenomic data using genome signatures[END_REF]. However, the results were found to be unreliable and not used further.

To normalize the β-lactamase count by a biologically relevant metric, we estimated the number of bacteria in our metagenomic dataset using 120 bacterial single copy genes from the Genome Taxonomy database using GTDB-Tk [START_REF] Chaumeil | GTDB-Tk: a toolkit to classify genomes with the Genome Taxonomy Database[END_REF] as reference (Table S4).

For comparison, the abundance of β-lactamases was also computed in all of the terrestrial metagenomics datasets available from the IMG/M database from June 2021 [START_REF] Chen | The IMG/M data management and analysis system v.6.0: new tools and advanced capabilities[END_REF].

1802 datasets were initially downloaded, then checked for contig length, coverage and quality. 434 datasets from known biomes and of sufficient sizes (at least one contig over 10kb and at least 3 bacteria equivalent sequenced) were retained and analysed with the exact same protocol used to detect ORFs, β-lactamases and bacterial single-copy genes in our own dataset (Table S5).

All statistical analyses were done using v4.1.1 of R (R Core Team 2021). Bray-Curtis dissimilarities were calculated with package vegan (v2.5-7) [START_REF] Oksanen | Vegan: Community Ecology Package[END_REF] and principal coordinate analysis performed with package stats4 (v4.1.1).

O R I G

I N A L U N E D I T E D M A N U S C R I P T
O R I G I N A L U N E D I T E D M A N U S C R I P T

Results

DNA extraction yield and quality

Large variations were observed in the DNA contents of the different soil samples, despite their similar macroscopic appearances (black/brownish fine silt and sandy compact soils), including between surface samples (taken from vegetation-free spots). Using the same extraction protocol, up to 75 times more DNA could be recovered from similarly dated ancient soils from different sites (Table S1, for instance sample B vs. P). A lesser range of variation (10 times) was observed between surface samples (Sample C to D) already suggesting that the sample age was not the main cause of these differences. Interestingly, very similar total numbers of usable base pairs and good quality reads were determined from the 250 ng of purified DNA used for sequencing) (Table S1). Except for sample B (sequenced as shorter reads on a different platform), we observed no correlation between the initial DNA content of the sample and the number of usable good quality reads (Table S1),

indicating that, once purified, all DNAs exhibited a similar quality.

Sample diversity

In contrast, the assembly of a similar number of sequenced reads from the different datasets, resulted in a highly variable number of contigs (≥5 kb). As these values are linked to the complexity of the sample (species richness-number of different species-and their relative abundance -species evenness-), this suggested that the samples exhibit globally different microbial population compositions and structures. As our samples are expected to be variable in richness and evenness while the number of reads is quite similar (except for sample B), the sequencing effort is expected to be different for each sample. This was confirmed by the Nonpareil tool that calculates the sequencing effort and an alternative to the Shannon diversity index computed from the kmer diversity of reads (Rodriguez-R et al.

2018).

Prior to any taxonomical annotation we could thus estimate the sample diversity as well as the estimated sample coverage (i.e. the proportion of the sample diversity actually sequenced) that ranges from 0.39 (sample E) to 0.9 (sample L) (Fig. S2). Given the same number of reads for two samples, a larger coverage corresponds to a lesser diversity.

According to Figure 1, the richest is sample E and the poorest is sample L. There are large differences between samples, including among the Kamchatka surface samples which appear the most diverse. The Nonpareil diversity values do not significantly correlate with the amounts of DNA recovered from one gram of soil (Table S1). In contrast and as expected, the

O R I G I N A L U N E D I T E D M A N U S C R I P T
diversity values negatively correlate with the median coverage of the contigs in each sample (Pearson correlation of -0.75, p< 0.005).

When focusing on a given borehole, we found that the diversity appeared to increase with sample depth in a given borehole (Fig. 1). However, diversity values were quite variable for samples taken at similar depths across different boreholes, in line with the lack of a consistent relationship between the diversity and the age (from radiocarbon dating) of the sample (Table 1A).

2018).

Global species contents

As depicted in Figure S3.A, bacteria account for over 90% of the estimated abundance in most samples and 87% in average but they represent 95% of contigs (> 5kb). This is true for all modern cryosols and most ancient permafrost layers despite their different origins and DNA contents. Even higher proportions of bacterial sequences were computed by Kraken2 [START_REF] Wood | Improved metagenomic analysis with Kraken 2[END_REF]. However, we found this method to be less sensitive and was not used further in this study. Four samples, L, N, R and Q appear to be poorer in bacterial

O R I G I N A L U N E D I T E D M A N U S C R I P T
sequences (<= 77%) due to their higher proportions of unclassified and Archaea sequences (Fig. S3.A).

In all samples combined, Archaea represent 3.2% of the abundance but only 1.8% of the total contig number. In most samples they are very rare except in samples Q (unfrozen talik) and 3 samples from different depths in the same borehole (L, R, N; respectively 12 m, 16 m, 19 m) where they represent 6.5 to 11.3% of the total abundance (Fig. S3.A). In the sample where they are in sizable amounts, the archaeal population is dominated by Bathyarcheota Eukaryotic DNA represents less than 1% of contigs and 0.5% of the total abundance, except for two sub-surface samples from the same borehole (R, N) with 1% and 2% of eukaryotes respectively (Fig. S3). The Viridiplantae (land plants and green algae) and Unikont (mostly Fungi, Metazoan and Amoebozoa) constitute most of the represented clades except for sample E and O (Fig. S5.A). The latter only exhibits 68 total eukaryotic contigs, 21 of which were attributed to the Bacillariophyceae (the class of diatoms). With 30 other contigs most likely originating from green algae, this suggests that this 40,000-year-old layer originated from a lake bottom deposit (Table 1A) (Fig. S5.B).

Finally, according to the LCA assignment method, viruses were predicted to only represent 0.4% of the global coverage of all samples (Fig. S3), with one sample (R) standing out by reaching 2.3%. Such low representation prompted us to reanalyze our datasets using Virsorter, a leading software tool specifically dedicated to the identification of viruses in metagenomics data [START_REF] Roux | VirSorter: mining viral signal from microbial genomic data[END_REF]. This new protocol increased the average abundance of virus to 2.68%, with only a handful of samples exhibiting values higher than 2%, including the R sample at 8.6% (Fig. S6). families known to gather giant or large DNA viruses. A detailed analysis of these virus populations will be published elsewhere.

A high diversity that points out syngenetically trapped communities

The global community analysis and an analysis of the kmer frequencies of reads in each dataset do not consistently put samples from the same borehole close to each other (L-N-R, P-K, Q-M-O) (Fig. S7). These analyses do not allow a meaningful clustering of samples originating from the different boreholes, or even distinguish surface samples from permafrost samples (Fig. S7).

Focus on the bacterial populations

The relative abundances of phyla representing more than 0.5% of the total sum of coverages are depicted in At lower taxonomical ranks the distributions become scattered, as well as less informative

given the smaller proportion of non-ambiguously annotated contigs (66.3% at the phylum level, 43.6 % at the order level, but only 16.4% at the genus level). The most abundant order overall is Rhizobiales that is present in all samples (Fig. 3). Unclassified Actinobacteria are the Actinobacteria represent the major group. In some samples (C, D, E, B, L, O), the high abundance of Proteobacteria is divided between multiples orders in the phylum rather than due to a predominant one. In samples R, N, P, K, M on the other hand, Proteobacteria are mostly represented by Rhizobiales (Fig. 3). Some of the orders found also in deep samples are known to be aerobic. This is the case for Vicinamibacteraceae (Huber and Overmann The abundance is estimated by the sum of coverages (normalized to 1 in each sample).

The analysis suggests a lot of variability across samples, making it illusory to expect that a tractable number of bacterial orders could be used to classify the various types of cryosol and permafrost into recurrent taxonomical types. While surface samples and most deep samples form separate clusters (Fig. 3) we should be careful drawing conclusions as the surface samples come from separate locations. When samples from the same borehole do cluster together (L-N-R and P-K) their dominant bacterial communities still exhibit some

O R I G I N A L U N E D I T E D M A N U S C R I P T
variations (Fig. 3). Interestingly, the large distance separating sample Q from the same borehole samples M and O is consistent with its unfrozen state as part of an under lake talik (Table 1A). Statistical analysis using ANOSIM confirms these observations (Fig. S7) showing that the distances between samples from the same borehole are not significantly lower than the distances between samples from different boreholes (R=0.388, p<0.069).

As a complimentary attempt to interpret the microbiome differences between samples, we screened the Pfam database, identified a total of 9481 motifs associated to specific protein families and/or molecular functions, and computed their distribution across samples. We expected that this type of analysis could reveal ecological/environmental constraints shared by various samples in spite of their different microbial communities. A list of the fifty most abundant Pfam motifs globally found in our dataset is given in Table S6. This list includes seven of the twenty most abundant motifs found in the whole Pfam database, supporting the validity of our annotation protocol. Wright 2019) we specifically investigated the presence and diversity of these proteins/enzymes in our dataset. Although the ancient microorganisms trapped in the pristine permafrost layers studied here have little probability to come in direct contact with human-beings, the risk remains that they might contribute to new β-lactam resistances once in contact with modern bacterial species.

A total of 1071 β-lactamase genes were detected by diamond blastP against Refseq and confirmed by a domain search (Fig. S9). Among them, 1066 were found in bacterial contigs, 4 in an archaeal contigs and 1 in an ambiguous contig. Figure 5 presents the numbers, origins and classes of bacterial β-lactamases. Although these numbers are small compared to the total number of contigs (Table S1), they nevertheless suggest that, in average, most (87%) of the sampled Bacteria harbour one copy of a β-lactamase gene, once normalized by the frequency of 120 known single copy genes detected in each sample (Fig. 5, Table S4). The normalized β-lactamase copy-number ranges from 0.39 *confidence interval: 0.22 -0.55+ per bacteria (sample N) to 1.75 *confidence interval: 1.54 -1.97+ (sample K) (Fig. 5). We noticed that the lowest prevalence of β-lactamase-encoding bacteria correspond to samples from the same boreholes (L-R-N, and Q-M), suggesting an actual ecological constraint. In most cases, we identified a single β-lactamase copy per contig, except for 17 contigs carrying several occurrences. In 8 of these, two β-lactamase genes were located next to each other and in 14 contigs, the two or three β-lactamase genes were found less than 10 ORFs away from each other.

Despite the large fluctuation in the average β-lactamase copy number, the relative proportions of the 4 enzyme classes remains fairly stable across samples. Class C *52-83%+

was always the most represented, followed by class B (14% in average), then Class A and D (10%) (Fig. 5). All classes were present in every sample except for sample E (no class A βlactamase).

Finally, we investigated the distribution of β-lactamase genes among the various bacterial phyla, as one given group of bacteria could concentrate most of β-lactamase genes. We thus simply compared the relative abundance of each bacterial phylum within the set of all cognate contigs with that of those exhibiting a β-lactamase gene. Bacteroidetes appeared to be the most enriched group in β-lactamase genes, while Planctomycetes exhibited the lowest proportion. Although the two distributions (Fig. 6) are different, they remain significantly correlated (Pearson coefficient > 0.8) pointing out that the occurrences of β-lactamase genes in these cryosols were shared among the diverse phyla constituting these microbiomes.

O R I G

I N A L U N E D I T E D M A N U S C R I P T

Comparing β-lactamase gene occurrences in the permafrost versus other soils types

To compare the β-lactamase gene proportion found in the above Russian cryosols to other types of soils, we extended our analysis to over 1500 terrestrial metagenomic assemblies publicly available from the JGI. As soil communities are usually very complex (e.g. (Rodriguez-R et al. 2018)), their assembly into large contigs is usually challenging [START_REF] Alteio | Complementary Metagenomic Approaches Improve Reconstruction of Microbial Diversity in a Forest Soil[END_REF].

Accordingly, only 434 data sets were found suitable to be further analysed by our protocol. If the average number of β-lactamase genes per bacteria was again found to be close to one, it was quite variable (mean: 1.09 ± 0.77) (Fig. 7). Interestingly, if our samples did not stand out, other permafrost samples from Alaska and north Sweden ranked high in the list (Fig. 7). The samples with most occurrences of β-lactamase genes per bacteria were agricultural soils followed by riparian soils. There is significantly more β-lactamase copy number in agricultural

O R I G I N A L U N E D I T E D M A N U S C R I P T
samples than in other terrestrial soils except for the JGI permafrost samples and riparian soils (Fig. 7). Bog and forest soils also exhibited β-lactamase gene copy numbers close to one.

The environments exhibiting the smallest β-lactamase gene copy numbers are shale carbon reservoirs, uranium contaminated soils and mire (a kind of wetland) (Fig. 7). We found the quantities of DNA recovered from the various cryosols to be highly variable, from more than 10 µg/g to 0.02 µg/g, thus a factor of 500 (Table S1). Most of the low yielding samples are from ancient (subsurface) permafrost layers, albeit without a strict correlation with age. The Stanchikovsky Yar sample (B) exhibited a particularly high DNA content. However, once purified, the recovered DNA exhibited similar sequencing yields (i.e. number of read/ng) and assembled into contigs of sizes similar to that of other samples (Table S1). This suggests that the variation in DNA content is not due to its degradation, but reflects real differences in the global amount of microorganisms per gram of sample. Such a result is significant because it supports the hypothesis that genes encoded by the ancient DNA could be recycled within contemporary microorganisms through transformation.

An unexpected increase of diversity with sample depth

The nonpareil diversity (a good correlate of α-diversity and of the Shannon index (Rodriguez-R et al. 2018) was found to increase with depth for each of the three boreholes. This is in contrast to the trend previously reported for culturable permafrost microorganisms [START_REF] Gilichinsky | Long-term preservation of microbial ecosystems in permafrost[END_REF]). As argued by Abramov et al. (2020) [START_REF] Abramov | Are permafrost microorganisms as old as permafrost?[END_REF], microorganisms trapped in permafrost do not multiply and, because of gradual death, their diversity is expected to decrease over time. They claimed that, although growth as subzero temperature remains possible for some microorganisms (for instance [START_REF] Panikov | Growth kinetics of microorganisms isolated from Alaskan soil and permafrost in solid media frozen down to -35°C: Growth kinetics in frozen media[END_REF][START_REF] Mykytczuk | Bacterial growth at -15 °C; molecular insights from the permafrost bacterium Planococcus halocryophilus Or1[END_REF]), the nutrient level and diffusion is too low to sustain growth as seen in the laboratory. The trend exhibited in our data between age and diversity is thus probably due to actual variations in diversity between syngenetically trapped microbiomes.

Scarcity of Eukaryotes

A striking finding of our taxonomic analyses is the very low representation of eukaryotic organisms, only present in trace amounts in all our samples except for two deep samples (R, N) (Fig. S3, Fig. S5). This result is unexpected given the fact that the warm spring-summer period sees the growth of a dense vegetation cover, accompanied by a variety of mosses and fungi, and a proliferation of insects, which should leave lasting traces in the metagenome of permafrost, once buried. Furthermore, the most abundant bacteria we found belong to the order Rhizobiales which are often in symbiosis with plants. The Siberian permafrost harbored

O R I G I N A L U N E D I T E D M A N U S C R I P T
also members of Burkholderiales (6% in average, up to 27%), Chitinophagales (3%, and up to 9%), and Cytophagales (0.7%), all of which are enriched in fungal-bacterial networks [START_REF] Bonito | Fungal-Bacterial Networks in the Populus Rhizobiome Are Impacted by Soil Properties and Host Genotype[END_REF]. A previous culture-based study of Siberian, Canadian and Antarctic permafrosts estimated the count of fungal cells around 10 2 per gram of dry soil, distributed among 49 species (Ivanushkina, Galina and Ozerskaya 2005). As fungi are in form of spores in permafrost [START_REF] Vorobyova | Micro-organisms and biomarkers in permafrost[END_REF], part of their diversity might remain unseen because of difficulties to extract their DNA. The proportion of eukaryotic cells might also be much lower than that of bacteria. In such case, their sequencing coverage might be insufficient to assemble contigs of a least 5kb in size. The larger size of eukaryotic genomes, the presence of introns, and their lower coding density due to large intergenic regions also strongly decrease the likelyhood of their identification by similarity searches in protein databases.

However, kmer-based read classification (Fig. S3.B) is consistent with an overall low proportion of eukaryotes. The remaining possibility is thus that eukaryotic DNA is rapidly recycled into microbial DNA during the seasonal decomposition process that takes place in the superficial active layer before being perennially frozen in deeper permafrost.

Variable bacterial community composition

As expected, bacteria constituted most of the DNA sources (Fig. S3). Compared to bar-coding methods, metagenomics allows an exploration of the sample diversity that is, in principle, only limited by a threshold of minimal abundance for a given species. This threshold results from 3 main factors: 1-the proportion of cells (or viral particles) of a given species (mostly unicellular) in each sample, 2-the efficiency with which their DNA can be extracted (e.g.

sporulating vs non-sporulating species), 3-our capacity to identify them by sequence similarities searches in reference databases. The last factor involves variations in proteincoding density (high in viruses and prokaryotes, much lower in most eukaryotes except for some fungi and protozoans), the fraction of the database corresponding to known members of the various phyla (high for Bacteria, much lower for viruses and Archaea), the level of sequence divergence among homologues within different domains (very high in viruses, lower in eukaryotes). Thus, even for a similar coverage (i.e. number of redundant reads per genome position) in a metagenomics mixture, bacteria are more likely to be identified than archaeal and eukaryotic organisms, viruses being the least likely to be identified (i.e.

taxonomically unclassified). The comparisons of organism frequencies across domains are thus much less reliable than within a given domain. Yet, the predominance of the bacterial domain is unlikely to be due to a bioinformatic/annotation bias. To estimate the abundance

O R I G I N A L U N E D I T E D M A N U S C R I P T
of a given group we used the sum of coverages of the contigs attributed to the group. In theory, a better estimation would have been the number of attributed reads relative to the genome size [START_REF] Nayfach | Toward Accurate and Quantitative Comparative Metagenomics[END_REF]. In the absence of information on the genome sizes, the sum of coverages thus remains a good approximation.

Globally, the phyla found to dominate the microbiomes in our cryosol and permafrost samples do not differ from those already reported as most abundant in previous analyses of temperate soils [START_REF] Janssen | Identifying the dominant soil bacterial taxa in libraries of 16S rRNA and 16S rRNA genes[END_REF] or permafrosts [START_REF] Tripathi | Distinct Taxonomic and Functional Profiles of the Microbiome Associated With Different Soil Horizons of a Moist Tussock Tundra in Alaska[END_REF][START_REF] Xue | Bacterial and Archaeal Metagenome-Assembled Genome Sequences from Svalbard Permafrost[END_REF]. Those are Proteobacteria (40%), Actinobacteria (22%), Bacteroidetes (10%), Firmicutes (9%), Acidobacteria (5%), Plantomycetes (4%) and Chloroflexi (3%) (Fig. 2). However, even the most common phyla do exhibit very different proportions between samples: 76% to 9% for Proteobacteria, 59.7% to 1.4% for Actinobacteria (Fig. 2). We noticed that the same phyla dominate both in the surface cryosols and the permafrost sample except for members of the Firmicutes that are mostly found below the surface (<1.3% at the surface, up to 28.5% in permafrost). In contrast to a previous study [START_REF] Tripathi | Distinct Taxonomic and Functional Profiles of the Microbiome Associated With Different Soil Horizons of a Moist Tussock Tundra in Alaska[END_REF] we did not detect a significant amount of Crenarchaeota. On the opposite, we detected Bathyarcheota, absent from this same study (Fig. S4).

The high taxonomical variability depicted here confirms that no significant mixing occurred between the distinct microbiomes of the permafrost layers, in particular for the Q, M, O samples extracted from the same borehole (Fig. 2 &4). The weak clustering of samples L, N, R and P, K using the Bray-Curtis dissimilarity might solely reflect the influence of the soil type or environmental conditions shaping bacterial communities.

Although dead cells would result in altered DNA and thus a lower efficiency of sequencing and assembling large enough contigs, it remains possible that most of the bacteria identified in our study might be dead. A previous metagenomic study of ancient Alaskan permafrost estimated that only 15 to 18% of cells were alive. However, the DNA depletion to remove dead cells did not significantly changed the relative proportion of bacterial groups (from phylum to family) or the phylogenetic diversity [START_REF] Burkert | Changes in the Active, Dead, and Dormant Microbial Community Structure across a Pleistocene Permafrost Chronosequence[END_REF].

There is a growing concern that taxonomy alone is insufficient to interpret the complexity of soil processes [START_REF] Baldrian | The known and the unknown in soil microbial ecology[END_REF] given the already high diversity of metabolisms and life-styles encountered in a given bacterial group. A good example of such a difficulty is the presence of Anaerolineales species in all our samples (Fig. 3). Although their anaerobic metabolism is consistent with life in deep permafrost, all known isolates exhibit optimal growth

O R I G I N A L U N E D I T E D M A N U S C R I P T
temperature in the 37°C to 55°C range [START_REF] Sekiguchi | Anaerolinea thermophila gen. nov., sp. nov. and Caldilinea aerophila gen. nov., sp. nov., novel filamentous thermophiles that represent a previously uncultured lineage of the domain Bacteria at the subphylum level[END_REF][START_REF] Yamada | Anaerolinea thermolimosa sp. nov., Levilinea saccharolytica gen. nov., sp. nov. and Leptolinea tardivitalis gen. nov., sp. nov., novel filamentous anaerobes, and description of the new classes Anaerolineae classis nov. and Caldilineae classis nov. in the bacterial phylum Chloroflexi[END_REF] or come from hot environments [START_REF] Konishi | Unique substrate specificity of a thermostable glycosyl hydrolase from an uncultured Anaerolinea, derived from bacterial mat on a subsurface geothermal water stream[END_REF]. On the other hand, these detected

Anaerolineales species appear to cohabit with close relatives to known psychrophilic bacteria such as Devosia psychrophila and Dyadobacter psychrophilus the presence of which is not unexpected in cryosols.

Thawing permafrost as a source of antibiotic resistance genes

Our study indicated the presence of the four main classes of β-lactamases [START_REF] Ambler | The structure of beta-lactamases[END_REF], the evolution of which can rapidly lead to ESBL in the context of strong antibiotic selection in clinical settings. Class A EBSL includes cephalosporinases and six types of dreaded carbapenemases [START_REF] Sawa | Molecular diversity of extended-spectrum β-lactamases and carbapenemases, and antimicrobial resistance[END_REF]. Class B β-lactamases are metalloenzymes using zinc at their active center compared with a serine residue for other classes.

They hydrolyse carbapenems and degrade all β-lactam agents except monobactams [START_REF] Sawa | Molecular diversity of extended-spectrum β-lactamases and carbapenemases, and antimicrobial resistance[END_REF]. Class C β-lactamases derives from the ampC gene carried on the genome of many Enterobacteriaceae. Variants of these enzymes are known to reduce sensitivity to carbapenems [START_REF] Sawa | Molecular diversity of extended-spectrum β-lactamases and carbapenemases, and antimicrobial resistance[END_REF]. Finally, class D βlactamases, also known as oxacillinases (OXAs) evolved from degrading an extended spectrum of cephalosporins to hydrolyse carbapenems [START_REF] Sawa | Molecular diversity of extended-spectrum β-lactamases and carbapenemases, and antimicrobial resistance[END_REF]. However, from their sequences alone (except for near 100% identical residues) it is impossible to predict the substrate range and the potential clinical risk associated to the diverse β-lactamases identified in our dataset.

The presence of antibiotic-resistant bacteria carrying various types of β-lactamases has been previously reported in pristine Arctic and Antarctic surface cryosols and ancient permafrost [START_REF] Allen | Functional metagenomics reveals diverse betalactamases in a remote Alaskan soil[END_REF][START_REF] Costa | Antibiotic resistance is ancient[END_REF][START_REF] Petrova | Genetic structure and biological properties of the first ancient multiresistance plasmid pKLH80 isolated from a permafrost bacterium[END_REF][START_REF] Perron | Functional Characterization of Bacteria Isolated from Ancient Arctic Soil Exposes Diverse Resistance Mechanisms to Modern Antibiotics[END_REF][START_REF] Kashuba | Ancient permafrost staphylococci carry antibiotic resistance genes[END_REF][START_REF] Van Goethem | A reservoir of 'historical' antibiotic resistance genes in remote pristine Antarctic soils[END_REF][START_REF] Haan | Unearthing Antibiotic Resistance Associated with Disturbance-Induced Permafrost Thaw in Interior Alaska[END_REF]. However, these previous studies did not evaluate the proportion of bacteria carrying a β-lactamase gene or exhibiting a β-lactam resistance in the global population of soil bacteria, most of which are unculturable (or eventually dead).

The number of identified β-lactamase genes was divided by the estimated number of distinct bacteria to obtain a normalized gene count. This original approach provides a biologically relevant quantification that can be generalized to target other enzymes and detect specifically enhanced metabolic features in prokaryotic communities. The usual approach

O R I G I N A L U N E D I T E D M A N U S C R I P T
normalizes the ORF count by the number of reads, and is less biologically relevant than the ORF count per bacteria.

We found few previous studies against which comparing our results about the presence of βlactamases in a quantitative manner. Almost all investigations of environmental "resistomes" were performed on populations of culturable bacteria sampled near the surface and selected for their antibiotic resistance (e.g. [START_REF] Haan | Unearthing Antibiotic Resistance Associated with Disturbance-Induced Permafrost Thaw in Interior Alaska[END_REF]. A few others were restricted to a given bacterial species (e.g. Escherichia coli) [START_REF] Pormohammad | Prevalence of antibiotic resistance in Escherichia coli strains simultaneously isolated from humans, animals, food, and the environment: a systematic review and meta-analysis[END_REF]. We The high proportion of β-lactamase-encoding bacteria found in our 12 cryosol samples initially appeared significantly greater than previously reported in the few publications concerning other pristine environments, undisturbed by anthropogenic activities. However, following the analysis of 434 supplementary data sets, the proportion of β-lactamases in our 12 samples came out near the average of its value for a variety of other soils, and for instance twice less than in some American agricultural lands (Fig. 7). Other permafrost samples turned out to be richer in β-lactamase than ours, some even overlapping with the level of cultivated soils (Fig. 7). The comparison of β-lactamase copy numbers between biomes revealed statistically significant differences (Fig. 7). Thus, β-lactamases are not distributed evenly among diverse environments under the influence of unknown factors that remain to be identified in future studies.

The function of β-lactamases in soil microbial communities is not clear, but our results might including members of the Actinobacteria (e.g. streptomyces species) [START_REF] Ogawara | Self-resistance in Streptomyces, with Special Reference to β-Lactam Antibiotics[END_REF]. It turns out that these bacteria are among the most abundant in the cryosol microbiomes (Fig. 2). Among eukaryotes, fungi (also abundant in our samples, Fig. S5) can also produce many complex β-lactam-related compounds. They are the original source of two foundational βlactam antibiotics: penicillin and cephalosporin [START_REF] Brakhage | Aspects on evolution of fungal beta-lactam biosynthesis gene clusters and recruitment of trans-acting factors[END_REF]). The high β-lactamase abundance might thus result from the selective advantage they confer to their bearers in the middle of a biochemical war between soil microbes.

Alternatively, β-lactamase may also play a role in a less violent scenario, outside of the simplistic "war" paradigm. Below the minimal inhibitory concentration, antibiotics may modulate bacterial gene expression, thus acting as signalling molecules rather than lethal weapons. In such context, β-lactamases might interfere with quorum sensing [START_REF] Yim | Antibiotics as signalling molecules[END_REF]. Another possibility is that some of these genes encode a broader function (such as carboxylesterases or DD-carboxypeptidases active on β-lactams [START_REF] Higgins | Characterization, cloning and sequence analysis of the inducible Ochrobactrum anthropi AmpC β-lactamase[END_REF][START_REF] Nan | A Novel VIII Carboxylesterase with High Hydrolytic Activity Against Ampicillin from a Soil Metagenomic Library[END_REF][START_REF] Pandey | MSMEG_2432 of Mycobacterium smegmatis mc2155 is a dual function enzyme that exhibits DD-carboxypeptidase and β-lactamase activities[END_REF]) providing a fitness gain in addition to protecting against soil antibiotics).

  Acanthamoeba-infecting viruses were isolated (and cultivated) from a wall-sampled layer of a well-(the Stanchikovsky Yar) (Sample B, Table

(

  talik) cores were retrieved. One from Yedoma deposits (61.75967°N, 130.47438°E), one below a Yedoma upland (61.76086°N, 130.47466°E), one from the adjacent drained Yukechi alas basin (61.76490°N, 130.46503°E) (Table

Fig 1 :

 1 Fig 1: Nonpareil diversity compared to sample depth. Nonpareil diversity was computed on reads over 35 bp in size. This metric is correlated to the Shannon index (Rodriguez-R et al.

(

  a recently described clade of methanogens belonging to the TACK phylum)[START_REF] Evans | Methane metabolism in the archaeal phylum Bathyarchaeota revealed by genome-centric metagenomics[END_REF], Nitrososphaeria (a class of chemolithoautotrophic ammonia oxidizing Archaea belonging to the Thaumarchaea and globally distributed in soils)[START_REF] Hatzenpichler | Physiology, and Niche Differentiation of Ammonia-Oxidizing Archaea[END_REF], and Methanomicrobia (a class of methanogens belonging to the Euryarchaeota and previously noticed to be enriched in the sediment microbiomes associated with the rhizosphere of macrophytes)(Behera et al. 2020) (Fig S4).

  The newly annotated viral (mostly phage) contigs originated from contigs attributed to Bacteria by the LCA method (Fig.S6.B). The family distributions between theses samples is quite variable but, in general, the most abundant viruses are Mimiviridae, Pithoviridae (the prototype of which was previously isolated from ancient permafrost[START_REF] Legendre | Thirty-thousand-year-old distant relative of giant icosahedral DNA viruses with a pandoravirus morphology[END_REF])), Caudovirales, Phycodnaviridae and Ascoviridae, all

Figure 2 .Fig 2 :

 22 Fig 2: Bacterial phyla abundances across samples. The abundance is estimated as the fraction of the total sum of coverages.

  second most abundant group suggesting that the cryosol Actinobacteria are very different from their cultured relatives. In samples taken individually, Rhizobiales is the most abundant order in samples N and R where they represent 53 and 23% of the bacterial abundance respectively. Burkholderiales are the most abundant order in sample E, D and O and Clostridiales are the most represented in samples Q and O. In the other samples, unclassified

Fig 3 :

 3 Fig 3: Bacterial abundance across samples. Samples and most abundant Bacterial orders (over 0.5% of the total normalized coverage) were clustered using the Bray-Curtis distance.

  Figure4.A presents a Bray-Curtis dissimilarity map computed from the Pfam counts in the samples. This map confirmed the functional similarity of the L-N-R samples, of the P-K samples (with the insertion of M from a different borehole, but from a similar depth and age as K), and to a lesser extent of two surface soil samples (D, E). The rest of the samples exhibited unrelated functional patterns, with the unfrozen talikderived sample Q as the most dissimilar.Table S6 also indicated the presence of 3 Pfam motifs, found in proteins involved in bacterial antibiotic resistance: TetR-like transcription regulators (PF00440, ranked 29 th /9481), Metallobeta-lactamases (PF00753, ranked 44 th ), and glyoxalases (PF00903, ranked 46 th ). This prompted us to focus our next analysis on Pfam motifs putatively associated with virulence factors, the result of which are presented as a heatmap in Figure4.B. In this graph, three different protein families associated to bacterial antibiotic resistance appear to be largely dominant. The top one corresponds to a large family of enzymes catalyzing different reactions, including glyoxalases, dioxygenases, and bleomycin/fosfomycin resistance proteins. Unfortunately, this Pfam motif is not specific enough to distinguish the minority of proteins truly involved in antibiotic inactivation from the rest of the family. Fortunately, the next most abundant functional families (β-lactamases and Metallo-β-lactamases) are much better defined and for the most part correspond to enzymes active against various types of β-lactam antibiotics. Their distribution is further analyzed in the next section. Finally, the less abundant TetR (tetracyclin resistance) family in the heatmap is defined by an ambiguous a large diversity of DNA-binding proteins most of which are not involved in antibiotic resistance.

Fig 4 :

 4 Fig 4: Analysis of the functional similarity between samples. (A) Bray-Curtis dissimilarity matrix computed from all Pfam counts in bacterial contigs. (B) Heatmap of the abundance of Pfam motifs potentially associated to bacterial virulence factors.

OFig. 5 :

 5 Fig. 5: Bacterial β-lactamases classes and copy numbers. The distribution of the 4 main Ambler's classes of β-lactamases in the various samples is shown by the left bars. "Class Crelated serine hydrolase" stands for proteins with a β-lactamase AmpC domain and whose best diamond blastP matches were annotated "β-lactamase related serine hydrolase". The number of bacterial single copy genes is shown by the adjacent grey right bar (mean count out of 120 reference single copy genes, with error bars corresponding to standard deviation). The boxes on top of each barplot indicate the average β-lactamase copy number per bacterial cell.

Fig. 6 :

 6 Fig. 6: Side-by-side comparison of the proportion of bacterial phyla in all samples with the distribution of β-lactamase carrying bacteria. These proportions are globally correlated (Pearson correlation coefficient > 0.8), indicating that β-lactamases genes do not originate from a single dominant phylum. Bacteroidetes appears the most enriched phylum in βlactamases.

Fig 7 :O

 7 Fig 7: β-lactamase gene copy numbers in available terrestrial metagenomes. The values are computed as the ratios between the β-lactamase gene counts and the average number of single copy genes identified in each sample by GTDB-tk. Boxes are ranked according to the median of the β-lactamase gene copy number. Only datasets with 3 or more estimated entire bacteria sequenced were included. The brown box corresponds to the 12 new metagenomes analysed in this work. The p-values of the pairwise comparison (Wilcoxon test) between our samples, agricultural lands and permafrost to other datasets are shown in the above heat map. Empty grey boxes correspond to non-significant differences.

  found one study analysing the proportion of β-lactam-resistant bacteria within different agricultural soils, untreated or amended with manure (Udikovic-[START_REF] Udikovic-Kolic | Bloom of resident antibiotic-resistant bacteria in soil following manure fertilization[END_REF]. This proportion was found to vary from 0.67% to 7.4% in untreated vs. manure-amended soils.Most of these resistant bacteria presumably encoded a β-lactamase, thus in a much lesser proportion than found in our cryosol samples. Another metagenomics analysis of the distribution of antibiotic resistance genes (ARG) in 17 pristine Antarctic surface soils estimated the relative frequency of ARG to all genes in the *1.3-4.4 10 -5 + range, which approximately correspond to 17% of soil bacterial encoding at least one ARG[START_REF] Van Goethem | A reservoir of 'historical' antibiotic resistance genes in remote pristine Antarctic soils[END_REF]) (assuming an average content of 3850 gene per bacteria (diCenzo and Finan 2017)). Moreover, this fraction includes all ARG (e.g. multidrug resistance efflux pump, aminoglycoside acetyltransferase/nucleotidylyltransferase, aminocoumarin resistant alanyl-tRNA synthetase, etc) of which β-lactamases only constitute a small proportion.

  help building new hypotheses. Several soil microorganisms produce antibiotic β-lactams,

  Unexpectedly, this study revealed a new potential danger, due to the abundance of βlactamases found encoded in a large phyletic diversity of cryosol bacteria. The DNA from these permafrost bacteria, dead or alive, thus constitutes an immense reservoir of historical antibiotic resistance genes. Their transfer to contemporary pathogenic bacterial in a clinical setting might further contribute to the antibiotic resistance crisis, now considered one of the biggest public health challenges of our time.In a context of global warming, ancient permafrost becomes an antibiotic resistance flavored ice cream ready to be consumed without moderation by all passing bacteria.

Table 1A : Deep permafrost samples analyzed in this study

 1A

	Sample	IGS Code	Latitude°	Longitude°	Locality	Depth (m)	Dating Years BP	Type
	P-1084T	B	68.37016	161.41555	Stanchikov-skiy Yar	-	35,000	Wall sampling, Yedoma bluff
	Yuk15-Yed-1 /Yed-6 Yuk15-Yed-1 /Yed-7	K P	61.75967	130.47438		16 11	>49,000 >36,000	Core, Yedoma upland
	Yuk15-Alas-1 / Ala-11	L				12	>28,000	Core, from a
	Yuk15-Alas-1 / Ala-10 Yuk15-Alas-1 / Ala-9	N R	61.76490	130.46503	Yukechi Alas landscape	16 19	>45,000 >42,000	drained thermokarst lake basin (Alas)
	Yuk15-YuL-15/Y1	O				19	>40,000	Core, Yedoma
	Yuk15-YuL-15/Y2	M	61.76086	130.47466		16	>48,500	upland, under lake (4 m water
	Yuk15-YuL-15/Y4	Q				6	53	depth)

# # : this soil sample was not frozen (talik). O, M, Q sample data from : Jongejans LL, Grosse G,

Ulrich M, et al. (2019)

Radiocarbon ages of talik sediments of an alas lake and a yedoma lake in the Yukechi Alas, Siberia. PANGAEA, https://doi.org/10.1594/PANGAEA.904738

Table 1B : Cold surface soils samples analyzed in this study

 1B

	Sample	IGS Code	Latitude°	Longitude°	Locality	Type
	P2	C	54.54972	160.58194	Kamchatka, Kronotsky river bank	Vegetation-free soil
	P4	D	55.09861	160.34944	Kamchatka, Kizimen volcano	Tundra soil
	P5	E	55.11500	159.96333	Kamchatka, Shapina river bank	Vegetation-free soil

  Downloaded from https://academic.oup.com/microlife/advance-article/doi/10.1093/femsml/uqac003/6564729 by SCDU Mediterranee user on 13 April 2022

Downloaded from https://academic.oup.com/microlife/advance-article/doi/10.1093/femsml/uqac003/6564729 by SCDU Mediterranee user on 13 April 2022

Acknowledgements

We are deeply indebted to our volunteer collaborator, Alexander Morawitz, for collecting the Kamchatka soil samples. We thank Dr. Audrey Lartigue for technical advices about DNA extraction, and acknowledge the computing support of the PACA Bioinfo platform. We thank M. Ulrich (DFG project #UL426/1-1) and P. Konstantinov for helping with fieldwork at the Yukechi site as well as the Alfred Wegener Institute and Melnikov Permafrost Institute logistics for field support and sample acquisition.

Funding

This work was supported by France Génomique (Grant #ANR-10-INBS-09-08) to J-M C and CNRS PRC research grant (#PRC1484-2018) to CA. E C-F was supported by a PhD grant (DGA/DS/MRIS #2017 60 0004). GG and JS were funded by ERC PETA-CARB (#338335) and the HGF Impulse and Networking Fund (ERC-0013).

Data availability

All raw reads are available through the EMBL-EBI under the study number: PRJEB47746 (secondary accession: ERP132049). Samples accession numbers: ERS7649018, ERS7649019, ERS7649020, ERS7649021, ERS7649022, ERS7649023, ERS7649024, ERS7649025, ERS7649026, ERS7649027, ERS7649028. Downloaded from https://academic.oup.com/microlife/advance-article/doi/10.1093/femsml/uqac003/6564729 by SCDU Mediterranee user on 13 April 2022

Competing Interests

The authors declare no competing financial interests

Contributions

JMC and CA initiated the project and designed the study, SR and SS performed the bioinformatics analyses, JMC wrote the initial version of the manuscript, EC processed the sample and performed DNA extraction, KL supervised DNA sequencing, GG, JS, and ANF collected the deep core samples from the Yukechi site. All authors contributed to the writing of the final manuscript.