Microbial and Geochemical Evidence of Permafrost Formation at Mamontova Gora and Syrdakh, Central Yakutia

1 Introduction

Permafrost soils contain around 1,000 Pg of organic carbon, close to 50% of its total terrestrial storage (Schuur et al., 2015). Microbial activity stimulated by increasing air and ground temperatures and associated permafrost degradation is likely to increase trace gas release and accelerate climate change, threatening global carbon goals (Natali et al., 2021). The microbial community structure of permafrost soils draws substantial attention as a potential control in carbon release to the atmosphere (McCalley et al., 2014; Hultman et al., 2015; Tveit et al., 2015; Brouillette, 2021; Emerson et al., 2021). Permafrost microbiota are frozen in permafrost, and its biological performance upon permafrost thaw corresponds to their strategies and metabolic versatility (Ernakovich et al., 2015; Mackelprang et al., 2017; Zhou et al., 2020; Perez-Mon et al., 2021). Metagenomics and 16S RNA sequencing are widely used in permafrost microbial community studies, although 40–50% of the permafrost-derived DNA is relict (Carini et al., 2017; Liang et al., 2019). Viability or metabolic activity of the derived cells cannot be concluded from DNA-based methods, it is shown that dead cell’s DNA has a minor effect on the community structure (Carini et al., 2017; Burkert et al., 2019). The Northern Hemisphere permafrost was subject to significant alterations in the last 650 ka (Konischev, 2011; Vaks et al., 2013; Murton et al., 2021a), potentially affecting the microbial communities. Taxonomic diversity and ecological functionality of microbiota in various permafrost strata could have either remained unchanged and reflect the environmental conditions around freezing time or undergone adaptation to potential thawing and refreezing.

Recent studies suggest that both composition and diversity of permafrost microbiota vary with the age of deposits (Mackelprang et al., 2017; Burkert et al., 2019; Liang et al., 2019), the ice content (Burkert et al., 2019), and are subject to dispersal limitation and both homogenous and variable selection (Bottos et al., 2018). The permafrost origin is additionally important in driving the community structure, since significant differences in community composition were reported from late Pleistocene lacustrine-alluvial and Ice Complex (Yedoma) sediments (Rivkina et al., 2016). Recently, the transition from Pleistocene to Holocene was shown to initiate a major threshold-type shift in the composition and structure of permafrost microbial community in Central Yukon (Saidi-Mehrabad et al., 2020). These results draw particular attention to the potential effect of climate change on microbial activity and trace gas production. Our study is a follow-up to numerous preceding research efforts in Central Yakutia, Northern Eurasia, focusing on microbiota and climate-relevant gas production (Cherbunina et al., 2018; Kim et al., 2019; Hughes-Allen et al., 2021). This region is underlain by thick ice-rich horizon, known as Ice Complex, a series of fine-grained Pleistocene deposits of disputed origin (Solovev, 1959; Péwé and Journaux, 1983; Ivanov, 1984; Sukhodrovsky, 2002; Vasil’chuk et al., 2004; Schirrmeister et al., 2020). Under changing climate, landscapes on Ice Complex deposits are subject to intensified thermal disturbance and thermokarst development (Sejourné et al., 2015; Saito et al., 2018; Zakharova et al., 2018).

Previous studies on Central Yakutia ICs using 16S rRNA gene sequencing from ice wedge material have shown lower diversity and dominance of anaerobic species and psychrophilic bacteria with older age (Rakitin et al., 2020). Electron microscopy of ice wedge-derived material combined with X-ray microanalysis evidence high diversity and organic origin of bacteriomorphic particles, including those with signs of low-temperature damage to the cellular structures (Filippova et al., 2014). Recently, bacterial communities from both Pleistocene ice wedges and Miocene alluvial sands have been described (Brouchkov et al., 2017; Ivanova et al., 2017; Filippova et al., 2019).

The microbial community structure in permafrost soils is supposed to be relatively stable over time (Shade et al., 2013), but rapidly shift to a new stable state under changing conditions, experiencing threshold behavior (Saidi-Mehrabad et al., 2020). Ground freezing and permafrost development are natural thresholds, but they can occur either simultaneously with sediment accumulation, or significantly later. Moreover, frozen ground can be subject to thermal degradation and subsequent refreezing, in which case the structure of microbial communities is expected to reflect natural conditions of the post-degradation period. Massive ice wedges are known to be formed predominantly from snow meltwater, but, other surface water sources might interfere as well, i.e., shallow groundwater or thermokarst lake water. New species associated with such water sources can be introduced to the initial microbial community and could serve, if identified, as markers of both thermal degradation and/or water origin. Ice-rich deposits of central Yakutia accumulated during a prolonged cold epoch from marine isotope stages (MIS) 4 to MIS 2 and were subject to partial thermal degradation during MIS 3 and during the Holocene, which should have been reflected in the microbial composition of these deposits and associated ice wedges.

In our study, our aim was to describe and quantify the microbial community composition and structure with 16S rRNA gene sequencing and bioinformatic metrics, and to relate these community features to the geological history of the “Ice Complex” (IC) deposits and the origin of massive ice wedges at two key locations in Central Yakutia, Northern Eurasia: Mamontova Gora and Syrdakh, both well described in the preceding literature (Ravsky et al., 1960; Baranova and Biske, 1964; Markov, 1973; Svitoch, 1983; Baranova et al., 1976; Péwé et al., 1977; Ershov, 1989; Lazukov, 1989).

Three hypotheses were tested in the present study: 1) the composition and structure of microbial communities are similar in host deposits and ice wedges that formed syngenetically, and strongly differ when ice wedges formed epigenetically; 2) the composition of microbial communities can be related to the water origin of the ice wedges, and indicate if the IC strata were subject to thawing and subsequent refreezing (epigenetic freezing); 3) the similarity in microbial community composition and structure can be indicative of similar age of the IC strata in different locations, or their contemporaneous accumulation and freezing.

2 Study Sites

Central Yakutia occupies the interfluve surface of the Lena and Aldan Rivers in North-East Siberia, Russia, underlain by thick continuous permafrost. Its thickness varies from 150 to 200 m in the middle Lena River valley, from 200 to 400 m under the accumulative terrace sequence on the right (east) bank of the Lena River, and up to 500 m under the left-bank (west) denudation plain (Ivanov, 1984; Nikolaev et al., 2011). The IC deposits, or ice-rich late Pleistocene sediments of disputed origin hosting thick ice wedges down to 50 m depth, and are widespread in the region (Fedorov and Konstantinov, 2009; Konischev, 2011; Schirrmeister et al., 2013). Numerous episodes of the IC thermal disturbance are reported at the Pleistocene-Holocene transition and during the Holocene optimum, leading to widespread thermokarst development (Fradkina et al., 2005; Katamura et al., 2006; Katamura et al., 2009; Pestryakova et al., 2012; Nazarova et al., 2013; Ulrich et al., 2017a; Ulrich et al., 2017b; Ulrich et al., 2019).

Data from the two major ground ice exposures are presented in this study: Mamontova Gora and Syrdakh (Figure 1; Supplementary Figures S1–S5). The well-described Mamontova Gora exposure (N63°01.169′, E133°55.787′) serves as a reference section for the Miocene-Pleistocene unconformity in East Siberian stratigraphy (Baranova and Biske, 1964). The Aldan River creates this exposure by undercutting the Central Yakutian terrace sequence at its northern margin, where six main terrace levels and a contemporary floodplain are observed (Markov, 1973). This study investigates the deposits of the 50 m terrace (V terrace), and to a lesser extent the 80-m terrace (VI terrace) and the right-bank floodplain of the Aldan River. The lowest part of the outcrop is the Miocene (N₂) sandy stratum which is exposed slightly above the shoreline, overlain by Middle Pleistocene (Q₂) sands (Figure 2). These middle Pleistocene sands are overlain by the late Pleistocene loams (Q₃) hosting two generations of thick ice wedges (Kuznetsov, 1976).

FIGURE 1. Location of the study sites (A) Study region (B) Schematic geomorphological map of the Lena-Aldan interfluve with main genetic types of deposits and corresponding terrains (based on Ivanov, 1984, Solovev, 1959; Pravkin et al., 2018; Fedorov et al., 2018). Legend: the types of terrain are indicated by color, geomorphological levels by hatching: 1) denudation plain, 2) glacial-fluvioglacial plain, 3) erosion-accumulative Abalakh plain, 4) low terraces of Lena and Aldan, 5) Kerdem terrace (III), 6) Bestyakhskaya terrace (IV), 7) Tyungyulun terrace (V), 8) Magan terrace.

FIGURE 2. Mamontova Gora outcrop: (A) Geological section across the Aldan River in the area of the V terrace of the Mamontova Gora outcrop. (B) thermocirque, exposing two generations of thick ice wedges (photo by D.G. Shmelev). Legend: 1) active layer, 2) icy dark gray loams—cover layer, 3) icy dark gray, black loams with organic matter and thick-schlieren cryostructure, 4) thawed clays and silts; 5) wedge ice, milky white, 6) wedge ice, muddy, 7) loam with massive cryostructure, gray and dark gray, 8) water (here its the Aldan R and it is Syrdakh lake for Figure 3), 9) peat, 10) sands medium to coarse-grained, 11) gravel-pebble deposits with sandy interlayers, 12) organic inclusions (plant remains), 13) wood remains, 14) pebbles and boulders, 15) pseudomorphs, 16) lithological bedding.

Syrdakh alas, a local name for thermokarst depression, is located on the Lena-Aldan interfluve, about 85 km to NE from Yakutsk (Figure 3). This alas has an elongated form, around 2 km long and 1 km wide, and depression depth is approx. 30 m. Previous studies have observed massive ice wedges close to the shoreline of the Syrdakh Lake. Four boreholes have encountered ice wedges at depths from 2.2 to 16.6 m, and their toes occurred at deepest at several meters below the lake bottom (Figures 3A,B). On inter-alas surfaces, the ice-wedge depth reached 45 and 34 m near the Oner and Syrdakh Lakes, respectively (Bosikov, 1985). Our study on Syrdakh site is based on a ice wedge outcrop (N62°32.638′, E130°57.915′), exposed in a thermo-erosional gully, around 6 m deep and 10–15 m wide.

FIGURE 3. Geological section of the Syrdakh Ice Complex: (A) schematic section of the coastal slope with archival boreholes (Bosikov, 1985); (B) geological section of the outcrop exposed in a thermo-erosional gully (C) the photo of the outcrop (by M. Cherbunina). Legend for the figure, see above (Figure 2).

3 Materials and Methods

3.1 Field Sampling

The field work was carried out at 2016–2017. To determine the composition and properties of ice and sediments, at least six monoliths were taken from one genetic horizon and were transported to the laboratory in frozen state (n = 70 for sediments and n = 18 for ice). Ice wedges from both IC horizons, as well as modern ice veins from the active layer of the Aldan River floodplain, and ice wedge from Syrdakh exposure, were sampled in vertical and horizontal profiles for stable water isotope analyses (Supplementary Figure S7). Ice samples were melted in the field, collected in 15 ml plastic vials without headspace, sealed and stored at cool place at near zero temperature (total n = 75, the results for most of them (n = 65) were previously discussed in Vasil’chuk et al. (2017) , Vasil’chuk et al. (2019) and in Budantseva and Vasil’chuk (2017). Microbiological samples were taken from the cleaned exposure wall, heated with butane blowtorch flame, carried over to sterile Whirl–Pak^® Nasco bags (Nasco, Modesto), and kept frozen at −5°C (mimicking natural conditions) in 36-L Coleman^® isothermal containers with saturated NaCl solution as a cooling agent until delivered to the laboratory. In the lab, samples were stored at −20°C in a freezer. Sampling point locations are shown on Supplementary Figure S6 (n = 10 for sediments and n = 2 for ice are presented in this study and 3 samples (n = 2 for ice and 1 for the active layer soil is previously published in Rakitin et al., 2020).

3.2 Laboratory Analysis

3.2.1 Physical and Chemical Analysis

For sediments, wet bulk density, further referred to as bulk density, was determined by directly measuring the volume and weight of frozen samples. Soil dry weight was measured after oven-drying the samples at 105°C for 24 h (Yershov, 1998). Gravimetric ice content was recalculated from the frozen sample weight and dry soil weight and is given as weight percentage (wt%). For grain size analysis, the samples were treated with hydrogen peroxide to remove organic material. Subsequently, the organic-free samples were diluted and washed to neutral pH values. The grain size distribution (GSD) was analysed optically using Fritsch Laser Particle Sizer Analysette 22 and was displayed in 62 size classes between 0.15 and 1,027.24 μm. GSD parameters, including the mean grain size diameter (MGSD) were calculated using the Gradistat software (Blott and Pye, 2020). The total carbon (TC) in soils and soil particles of the ice (in wt%) was detected by an elemental analyzer (Vario EL CHNS analyser, Germany) The standard deviation for TC was ±0.1% for repeated measurements. Inorganic carbon content in Quaternary deposits of North-Eastern Russia is typically less than 10–15% of TC (Schirrmeister et al., 2011), while recent studies show larger values up to 22% at Batagay IC (Shepelev at al., 2020).

Ice wedge samples were melted at room temperature and filtered through pre-weighed filters. Filters with retained sediment were dried at 105°C for 24 h and weighted, sediment dry weight was recorded. The pH values were measured electrometrically (US EPA, 2017) using an Expert 001 ion meter (Econix-Expert Ltd., Russia). Total content of dissolved inorganic anions (Cl⁻, SO₄²⁻, HCO₃⁻), and cations (Ca²⁺, Mg⁺, Na⁺, K⁺) in ice filtrates and aqueous dry soil suspensions was measured using Russian standard methods. Hydrocarbonates content (GOST 31957-2012, 2019) was recalculated from total alkalinity, measured by titration with 20% accuracy. Major anions were measured by capillary ion electrophoresis (GOST 31867-2012, 2019; GOST 31869-2012, 2019), with accuracy ca. 15%, using ‘Kapel’ capillary ion electrophoresis system (Lumex, Russia). Analytical results were expressed in mg/L and recalculated to milliequivalents per liter (meq/L) as per standard practices. Then, equivalent masses were totalled separately across anions and cations and equalled 100% each, from where %-equivalent masses (%-eq.) were calculated for each anion and cation.

3.2.2 Radiocarbon Dating

Radiocarbon dating of the two samples from Mamontova Gora and Syrdakh IC was done by the Laboratory of Radiocarbon Dating and Electron Microscopy of the Institute of Geography of the Russian Academy of Sciences (IGAN) and the Center for Applied Isotope Research, University of Georgia (Athens, United States) using accelerator mass spectrometry (AMS). The calibration programme CALIB REV7.1.0 using the IntCal13 curve (Reimer et al., 2013), was used to calibrate these radiocarbon dates. The dates obtained at various facilities in the 1980’s and 1990’s and published without calibration were calibrated in Calib 8.2 (Stuiver et al., 2021) using the IntCal20 calibration datasets (Reimer et al., 2020).

3.2.3 Stable Water Isotopic Composition

Stable water isotope analysis was done by isotope ratio mass spectrometry with constant helium flux (CF-IRMS), with Finnigan Delta-V Plus mass spectrometer using gas-bench device, at the isotope laboratory of the Department of Geography, Moscow State University (MSU). Calibration was made using the V-SMOW and SLAP international standards, as well as the MSU internal standard, snow from the Garabashi glacier (Mount Elbrus area), with δ¹⁸O = −15.60‰, δ²H = −110.0‰. The measurement errors were ±0.6 and ±0.1‰ for δ²H and δ¹⁸O, respectively.

3.2.4 Profiling of Prokaryotic Communities Based on 16S rRNA Gene

gDNA (genomic DNA) isolation of samples C-3–C-12 (Table 1; Supplementary Figure S6) and primary amplification of the studied samples, depending on the type and age of the samples, was carried out in the “clean” rooms of the Institute of Geophysics of the Environment (IGOS) of the University of Grenoble-Alpes in Grenoble, France (IGE, CNRS-UGA, Grenoble) (Bulat et al., 2004). Cell disruption and gDNA isolation were performed mechanically using the FastPrep instrument (MPBiomedicals, United States) and the PowerSoilDNAisolationKIT (MoBioLabs, United States) with E matrix (beads) (MPBiomedicals, United States), as well as using the FastDNASpinKitforSoil kit according to the manufacturer’s method (MPBiomedicals, United States). The concentration was measured using a Qubit 2.0 fluorometer with a dsDNAHSreagentKIT (Invitrogen™, United States).

TABLE 1. The list of samples for microbiological analysis of the outcrops of Mamontova Gora and Syrdakh.

DNA for other samples (led-S, C-12a showed at Table 1; Supplementary Figure S6) was isolated using the FastDNA SpinKit and FastPrep-24 bead beating grinder (MP Bio, United States) and according to the manufacturer’s instructions. The libraries of the V4 region of the 16S rRNA gene for Illumina MiSeq high-throughput sequencing were prepared using PCR with following primer system: forward (5′-CAA​GCA​GAA​GAC​GGC​ATA​CGA​GAT​GTG​ACT​GGA​GTT​CAG​ACG​TGT​GCT​CTT​CCG​ATC​T XXXXXX ZZZZ GTGBCAGCMGCCGCGGTAA-3′), containing, respectively, the 5′ Illumina Linker Sequence, Index1, Heterogeneity Spacer (Fadrosh et al., 2014), and 515F primer sequence (Hugerth et al., 2014); and the reverse primer (5′-AAT​GAT​ACG​GCG​ACC​ACC​GAG​ATC​TAC​ACT​CTT​TCC​CTA​CAC​GAC​GCT​CTT​CCG​ATC​T XXXXXX ZZZZ GACTACNVGGGTMTCTAATCC-3′), containing the 3′ Illumina Linker Sequence, Index 2, Heterogeneity Space, and the Pro-mod-805R primer sequence, respectively. This primer set covers 86% of 16S rRNA gene sequence diversity of the Archaea and 84.3% of the Bacteria (Merkel et al., 2019). For each DNA sample, two libraries were prepared, which were sequenced in parallel using the MiSeq Reagent Micro Kit v2 (300-cycles) MS-103-1002 (Illumina, San Diego, CA, United States) on a MiSeq sequencer (Illumina, San Diego, CA, United States) according to the manufacturer’s instructions. The primary processing of raw reads was carried out as described earlier (Kallistova et al., 2020).

3.3 Data Treatment

All 16S rRNA gene sequence reads were processed by the SILVAngs 1.4 pipeline (Quast et al., 2013) using the default settings: 98% similarity threshold was used for creating operational taxonomic units (OTUs) tables, and 93% was the minimal similarity to the closest relative that was used for classification (other reads were assigned as “No Relative”). Abundance data were further treated in RStudio software (RStudio Team, 2021), an open-source GUI to R programming language (RCore Team, 2021). Statistical and graphical analyses relied on packages “tidyverse” (Wickham et al., 2019), “ggplot2” (Wickham, 2016). Dissimilarity and similarity-based distance metrics, function vegdist(), package “vegan” (Oksanen et al., 2020). Heatmaps were prepared using function pheatmap(), package “pheatmap” (Kolde, 2018). UniFrac and Principle Coordinate Analysis (PCoA) were also carried out in QIIME2 via q2-diversity function (Lozupone and Knight, 2005; Lozupone et al., 2007; Lozupone et al., 2011). The significance of weighted and unweighted unifrac matrices was calculated with Permanova test (Anderson, 2001).

4 Results

4.1 Stratigraphy

4.1.1 Mamontova Gora

At the Mamontova Gora site, an outcrop of the V terrace of Aldan River was studied in one location where numerous thermal denudation features, thaw slump scars, up to 20 m deep and about 200 m in diameter, exposed a two-horizons ice wedge structure (Figure 2). The complete outcrop profile including the underlying sediments is shown in Figure 4. In this profile, the following strata were described, top to bottom [as meters below surface (m bs)]:

FIGURE 4. Characteristics of the Mamontova Gora outcrop (Vterrace). Legend: 1)pebbles, 2)sand, 3)silt, 4)loam, 5)boulders and gravels, 6)peat and organic inclusions, 7)ice layers, 8)lithological layering, 9)pseudomorphs, 10)ice wedges, 11)total carbon content is the sediments, 12)total carbon content in soil particles of the ice wedges (relative to the dried soil particles).

Unit 1 (0–1.5 m bs): Silty brown loams, thawed from the surface down to 0.9 m (active layer depth at the time of description). The ice content in the frozen part of the active layer ranges between 27 and 31 wt%, decreasing towards the surface. The stratum is heterogenous, and there is clear separation in color and TC content between the active layer, and transient layer (Shur et al., 2005). The TC content is between 2.5 and 3.5 wt% in the transient layer and varies from 1.4 to 2.0 wt% in the active layer.

Unit 2 (1.5–3.7 m bs): “Ice Complex I” (ice wedges and host deposits, upper horizon). The host deposits are dark grey, heavy, silty loams with thick-schlieren cryostructure and with peat interbeds locally. The ice content in loams is from 45 to 50 wt%, TC from 4.1 to 4.2 wt%; peat layers have higher ice content, from 129 to 160 wt% and higher TC, from 6.6 to 6.8 wt%. Particulate TC content in soil particles of the ice wedge ranges from 3.7 to 4.2 wt%, relative to the dried soil particles. This unit gave a median ¹⁴C date of 43.3 cal ka BP (this study, Supplementary Table S1). The 0.5–2 m wide and up to 3 m long wedges penetrate into the lower horizon.

Unit 3 (3.7–15.7 m bs): “Ice Complex II” (ice wedges and host deposits, lower horizon). Host deposits are dark grey, ice-rich, silty loams and grey sandy loams. Ice wedges are up 5 m thick, the size of the polygonal grid is about 10–15 m. Syngenetic origin is indicated by layered cryostructures in the host sediments that turn upward adjacent to the wedge owing to the ice growth which was also noted earlier by Popp et al. (2006). Ice content varies from 49 to 66 wt%, TC from 2.1 to 4.2 wt%. The soil particles fraction content changes from 1.9 to 2.4 vol%, particulate TC content in soil particles of the ice wedge ranges from 2.1 to 3.2 wt% and decreases down the profile.

Unit 4 (15.7–19.0 m bs): Blueish-grey heavy ferruginous loams with high peat content, covered with a film of 1–2 mm similar to biogenic mats when interacting with air. Ice content is from 28 to 35 wt%, TC from 3.4 to 4.0 wt%. The lattice-like cryostructures are observed here with cell size of 3–4 by 7–10 cm.

Unit 5 (19.0–38.4 m bs): Light grey sands, mainly coarse-grained, with horizontal layers and cross-bedding structures, ferruginous nodules, and interbedded with dark grey sandy loams. Ice content is low, varying from 18 to 24 wt%. Cryostructure is massive for the sands, and horizontally layered for the loams, with a width of 1–3 mm and a distance between them of 5–8 mm.

Unit 6 (38.4–46.2 m bs): Well-rounded gravel-pebble deposits with sandy interbeds. Ice content is from 18 to 20 wt%, TC from 0.9 to 1.4 wt%.

Unit 7 (below 46.2 m bs): Yellowish-grey mixed sand, from fine to coarse, with massive cryostructure containing wood residues, pebbles, and ferruginous interlayers. Ice content varies from 19 to 32 wt%, TC ranges from 0.7 to 4.3%. Ice wedge casts were described previously in this stratum (Markov, 1973).

The highest VI Aldan River terrace was not studied in detail during our field expedition. Only basal sandy deposits were sampled for microbiological analyses, as their appearance is close to the lowest layer (Unit 7 above) of the V terrace. Its ice content varies from 19 to 31 wt%, TC from 0.9 to 3.5 wt%. Therefore, samples were taken here to allow inter-comparison of microbial composition between two terrace levels and similar deposits.

4.1.2 Syrdakh

At the Syrdakh site, two main strata were exposed in the visible section of the profile (Figure 5):

FIGURE 5. Characteristics of the Syrdakh outcrop. Legend for the figure see above (Figure 4).

Unit 1 (0–2.7 m bs): Interbedded loams and sandy loams. The active layer at the time of observations was 1.2 m deep with inclusions of poorly decomposed organic matter in the form of a bluish felt. The ice content in the frozen part of the active layer varies from 25 to 29 wt%, TC is 2.0 wt%. The transient layer has abundant peat inclusions, and TC content around 2.1 to 2.2 wt%. Cryostructure is massive, with thick ice layers toward the layer base. This layer gave a median ¹⁴C date of 10.6 cal ka BP (this study, Supplementary Table S1).

Unit 2 (2.7–5.1 m bs): “Ice Complex” (ice wedges and host sediments). Ice wedges have apparent vertical layering. Host sediments are loams with schlieren and the lattice-like cryostructures. Ice content ranges between 46 and 65 wt%, the TC content ranges from 1.2 to 2.1 wt%.

4.2 Water and Sediment Chemistry

At Mamontova Gora, samples were taken from two Ice Complex (IC) horizons, and from ice veins in the contemporary Aldan River floodplain ca. 4–5 m above the low-flow water stage. At the Syrdakh site, a single sample was collected from the exposed ice wedge. The two IC horizons differ in their chemical composition (Figure 6), as previously discussed by Kuznetsov (1976). The total dissolved solids (TDS) content varies from 393 to 687 mg L⁻¹ and is slightly higher in the upper IC layer. Hydrochemical type of the “Ice Complex I” wedge ice is HCO₃-Mg (n = 1), “Ice Complex II” ice, HCO₃-Ca (n = 2) or HCO₃-Na-K (n = 3). Contemporary ice vein of the Aldan River floodplain is HCO₃-Na-K (n = 1). We have also observed high Cl⁻ content in the ‘Ice Complex II’ ice, up to 40%-eq.

FIGURE 6. Piper diagram of the chemical composition of ice wedges of the Mamontova Gora outcrop according to the authors (this article) and Kuznetsov (1976), %-equiv. 1) ice wedge of the upper tier (our results), 2) ice wedge of the upper tier (Kuznetsov, 1976), 3) ice wedge of the lower tier (our results), 4) ice wedge of the lower tier (Kuznetsov, 1976), 5) modern ice of the floodplain (our results), 6) Holocene ice of a 6–9 m terrace (Kuznetsov, 1976).

Chemistry of aqueous soil extracts prepared from sediments is uniform across the studied profiles. TDS content ranges from 90 to 390 mg per 100 g dried soil in the “Ice Complex I” host deposits, from 60 to 350 mg in the “Ice Complex II” deposits, and from 100 to 180 mg in the I (9 m) Aldan River terrace sediments (Kuznetsov, 1976). The water type is HCO₃-Ca, rarely HCO₃-Mg “Ice Complex I” or HCO₃-Na-K (“Ice Complex II”).

Water type of the aqueous soil extract in samples collected from Syrdakh exposure is HCO₃-Na-K in the upper horizons, and HCO₃-Ca in the lower horizons, surrounding the ice wedges, and with lower TDS compared to the Mamontova Gora samples, from 51 to 121 mg per 100 g dried soil.

Increased sodium content is typical for ice-containing sediments of the region (Anisimova and Pavlova, 2014). Water and sediment chemistry of the Mamontova Gora exposure was studied previously (Kuznetsov, 1976; Vasil’chuk et al., 2004; Vasil’chuk et al., 2017) as one of the proxies of ice origin and tracers of landscape evolution. According to Kuznetsov (1976), ice wedges of the V Aldan River terrace, classified as late Pleistocene, have higher TDS content, from 230 to 640 mg L⁻¹, than lower-lying wedges of the I terrace, or Holocene wedges, from 71 to 219 mg L⁻¹. Dominant water type in late Pleistocene ice wedges was HCO₃-Mg, in Holocene ice wedges, HCO₃-Ca. Total dissolved solids content in samples collected from an “Ice Complex II” late Pleistocene ice wedge (Vasil’chuk et al., 2017), varied from 80 to 476 mg L⁻¹. Hydrochemical type was HCO₃-Ca, switching to HCO₃-Mg between 7.0 and 7.8 m from the surface. The pH value varied over a wide range, from 4.4 to 7.6.

4.3 Water Stable Isotope Composition

Both upper and lower ice wedge generations of the V Aldan River terrace, a contemporary ice wedge of the Aldan River floodplain, as well as ice wedge of the Syrdakh outcrop were sampled. The results are partially presented and discussed in (Vasil’chuk, 1988; Vasil’chuk, 1992; Budantseva and Vasil’chuk, 2017; Vasil’chuk et al., 2017; Vasil’chuk et al., 2019), therefore here only a brief overview is given adding data on upper ice wedge horizon.

The upper “Ice Complex I” ice wedges have δ¹⁸O values from −29.6 to −24.4‰ and −28.8 to −27.4‰ for two profiles δ²H values from −227 to −192‰, −219 to −207‰ with a slight difference in variance between the profiles (Figure 7; Supplementary Figure S7B). The atmospheric origin of ice with insignificant participation of non-atmospheric waters was previously suggested based on the one of the profiles (Vasil’chuk et al., 2019). Nonetheless, in several samples highly enriched in heavy isotopes, deuterium excess, or d-excess, descends below 3.0, which may indicate participation of surface water from active layer or mixed snowmelt and shallow marshes waters that underwent sufficiently strong evaporative loss.

FIGURE 7. Ratio of δ¹⁸O and δ²H values in late Pleistocene and Holocene ice wedges of an outcrop of a V terrace on the right bank of the Aldan River (Mamontova Gora outcrop). Sampling profiles 2 and 3 refer to the lower tier of the ice complex, 6 and 5 to the upper tier. GMWL, global meteoric water line.

The lower “Ice Complex II” ice wedges show highly variable stable isotope composition (Figure 7; Supplementary Figure S7B). The δ¹⁸O values from −29.58 to −24.69‰ were observed in one location, while from −30.89 to −27.89‰ in the other location about 3 m apart, close to δ¹⁸O from −31.5 to −28.5‰, previously reported in (Popp et al., 2006). Isotopically heavier values were obtained for the upper parts of the ice wedge, conflicting with previously published results, where δ¹⁸O values from −29 to −25.9‰ were obtained in the upper part of the wedge, and from −22.7 to −16.5‰ in the high-TDS bottom toes of an ice wedge (Vasil’chuk and Vasil’chuk, 1998).

Holocene and contemporary ice wedges from the Aldan River floodplain are generally isotopically heavier, with δ¹⁸O values from −25.9 to −23.2‰, and δ²H values from −196 to −178‰ (Figure 7; Supplementary Figure S7C).

At Syrdakh, along the vertical profile, δ¹⁸O and δ²H values varied narrowly in a range from −31.8 to −30‰ and from −251 to −231‰, respectively (Figure 8; Supplementary Figure S7D). Horizontal δ¹⁸O variations were in the range from −32.5 to −29.2% at a 0.5 m depth, and from −31.5 to −30.5 at a 1.2 m depth. Isotopic composition of Syrdakh ice suggest its atmospheric origin and is in the typical regional range for snow. At the same time, d-excess values vary from 3 to 12‰ and fall below 5‰ on multiple occasions, which reflects the relatively frequent participation of non-meteoric water in ice formation, for example, evaporated water from interpolygonal ponds.

FIGURE 8. Ratio of δ¹⁸O and δ²H values in the late Pleistocene ice wedges of the Syrdakh outcrop. GMWL, global meteoric water line.

4.4 Prokaryotic Communities

A sample list with additional information is given in Table 1. Characteristics of amplicon libraries and alpha diversity of prokaryotic communities are summarized in Table 2.

TABLE 2. Characteristics of the amplicon libraries of the Mamontova Gora and Syrdakh.

Microbial communities of the studied samples are mainly represented by bacteria. Archaea are minor and constitute no more than 0.6% in all samples. The bacterial communities differ both in the composition of dominant groups and in their ratio. Each of the nine phyla of bacteria exceeded 3% abundance in one or several studied geological units: Acidobacteria, Actinobacteria, Bacteroidetes, Chloroflexi, Firmicutes, Gemmatimonadetes, Proteobacteria, Parcubacteria, Verrucomicrobia (Figure 9).

FIGURE 9. Relative abundance of prokaryotic phyla (A) and classes (B) in different horizons of Mammoth Mountain and Syrdakh (the samples AL, MMG and MMP correspond to previously studied samples of the ice complex of Mamontova Gora:an active layer (AL), an ice wedge of the upper (MMG) and lower horizons (MMP) (Rakitin et al., 2020).

A number of reads per sample varied from ∼4,200 and ∼14,300 after using all quality filters. The coverage for the studied microbial communities ranged from 0.58 to 0.80, estimated by Good’s coverage index, meaning that between 20 and 42% of reads were from OTUs that appear only once in the samples. Such coverage is insufficient for a complete description of the phylogenetic diversity of studied microbial communities. However, all major existing groups have been identified. The diversity of prokaryotes was estimated with the Shannon index, which varied within a wide range from very moderate to relatively high values, from 2 to over 5.2. Sampling points are shown in Supplementary Figure S6.

4.4.1 Mamontova Gora

Miocene Alluvial Sands (Samples С-3, С-5 From the “Unit 7” of V Terrace, С-6 From the VI Terrace)

The dominant phylum in all three samples was Proteobacteria (58–61%), represented mainly by the class Gammaroteoabcteria (44–58%), the other predominant phyla were Bacteroidota (3–24%) and Actinobacteriota (7–24%). Firmicutes constitute a significant part only in C-6 community (9%) whereas C-5 was characterized by a high presence of the Patescibacteria phylum (14%). At the genus level, the main feature of these three samples is high abundance of Gallionella (6–19%) and Sideroxydans (3–6%), which are almost completely absent in other samples. A distinctive feature of samples C-3 and C-5 is the presence of Salinibacterium (5–6%).

Middle Pleistocene Alluvial Sands, Overlying Miocene Sands (Samples C-7, C-8 From Unit 5)

There are three dominant phyla in these samples: Firmicutes (31–43%), Bacteroidota (34–35%) and Proteobacteria (17–31%), in both cases almost equally represented by Gammaproteobacteria and Alphaproteobacteria. Phylum Firmicutes is almost entirely represented by Desulfosporosinus genus and this is the characteristic feature of these samples. Phylum Bacteroidota is exclusively represented by uncultivated microorganisms of env.OPS 17 group which is typical for other samples where this phylum is highly present. Another distinctive feature of these samples is the relatively high abundance of Diaphorobacter (7–10%).

“Ice Complex II” Late Pleistocene Deposits (Samples C-4, C-11) From Unit 3

The late Pleistocene samples are characterized by a high degree of community complexity, with Shannon Index from 4.6 to 4.8, mainly represented by Actinobacteriota (20–30%), Proteobacteria (22–32%), Firmicutes (21–38%), and Bacteroidota (7–8%). At the genus level, these two samples are very different from each other. For example, the most common genus in sample С-4 is Clostridium (16%), whereas in sample С11, Bacillus (21%).

“Ice Complex I” Deposits of Uncertain Age (Late Pleistocene or Holocene) (Samples С-9, С-10 From Unit 2)

The sample C-10 was taken from the same ice wedge of the upper IC horizon as the sample “MMG” described in Rakitin et al. (2020), but closer to the base of the ice wedge. The community of this sample is the most complex among studied, Shannon Index value is 5.2. The phyla Actinobacteriota (39%), Proteobacteria (23%), Chloroflexi (11%), Bacteroidota (7%), and Firmicutes (5%) form the backbone of the community. Unlike sample C-9, there are no clearly dominant genera of microorganisms here. Each genus makes up no more than 4% of all microorganisms. The sample C-9, obtained from hosti loams, has the same dominant phyla but in different proportions: Actinobacteria (7%), Proteobacteria (6%), Chloroflexi (30%), Bacteroidota (7%), and Firmicutes (27%). At the genus level, there are two clearly dominant groups: uncultured Chloroflexi of the Gitt-GS-136 cluster (24%) and Bacillus genus (12%).

4.4.2 Syrdakh

Ice Complex Late Pleistocene Ice Wedge (Sample Led-S From Unit 2)

The Ice Complex sample from Syrdakh exposure was dominated by the phyla Actinobacteriota (30%), Firmicutes (23%), Chloroflexi (21%), and Proteobacteria (19%). In contrast to other studied samples, phylum Bacteroidota is almost absent from this sample. Another distinctive feature of this community is the abundance of phylogenetically deep cluster of uncultured Actinobacteriota named MB-A2-108 (8%). Besides, at the genus level there are several widely represented groups: Bacillus (9%), uncultured Chloroflexi of the Gitt-GS-136 cluster (13%), uncultured Chloroflexi of the KD4-96 cluster (6%).

Overlying Loams (Samples С-12, С-12a From Unit 1)

The phyla Chloroflexi (15–30%), Actinobacteriota (26–28%), Proteobacteria (7–22%), Acidobacteriota (4–6%), and Bacteroidetes (2–3%) are dominant in these samples. Also, phylum Gemmatimonadota (12–19%) is abundant in these samples, unlike other studied communities. This phylum is almost entirely represented by uncultivated microorganisms. Among Actinobacteria, a significant part of the community in both samples is represented by uncultivated groups of Actinobacteria, Thermoleophilia and Acidimicrobiia classes. At the genus level, uncultured Chloroflexi of the Gitt-GS-136 cluster (7–19%) are also highly abundant.

5 Discussion

5.1 Age and Origin of the Mamontova Gora and Syrdakh Deposits

5.1.1 Geological Evidence

Starting from Miocene, the Mamontova Gora region experienced tectonic subsidence compensated by sedimentation, and a monotonous sandy stratum accumulated since that time (Ivanova et al., 2015). The sands of the VI Aldan River terrace do not contain cryogenic formations. Their age is considered pre-Quaternary (Baranova and Biske, 1964), while others subdivide this stratum into three sub-horizons: Miocene, transitional and lower Pliocene (Ravsky, 1960). These same sands are exposed at the base of the V Aldan River terrace. Miocene flora from the Mamontova Gora basal sands dates to Middle Miocene (Nikitin, 2007). During that time, according to paleobotanical data, the climate was warm, with mean annual temperature about +12°C, mild winters with frosts not lower than −20°C, and precipitation from 1,000 to 1,500 mm (Biske and Baranova, 1976).

In the late Pliocene, the climate was still humid, but winters were getting colder, it is assumed that mean annual temperatures were below 3°С. Permafrost was most likely absent in the region during this time (Fradkina et al., 2005; see Supplementary Figures S8, S9). Stable depositional environment shifted to tectonic stability or, possibly, uplift and erosion. Monotonous cross-bedded sands are overlain by a pebble horizon, with erosional base contact and multiple unconformities. The similar pebble horizon is dated back to Middle Neopleistocene in the close Chuyskaya Gora exposure on the right bank of the Aldan River ca. 400 km downstream from Mamontova Gora (Ivanova et al., 2015). In the Batagay megaslump exposure at least one episode of permafrost thaw and erosion occurred sometime between MIS 16 and 6 (Murton et al., 2021b).

Markov (1973) attributes the accumulation of sediments of the upper part of the VI terrace level to late Pliocene, while the alluvium of the V terrace accumulated during the Middle Pleistocene. These alluvial layers contain traces of syngenetic frost cracking (Katasonov and Ivanov, 1973; Ershov, 1989), as well as wedge ice in the upper part of these sands, directly below the overlying thermokarst basin deposits (Markov, 1973). A single available U/Th date, 300 ± 5.7 ka BP, was obtained from loamy interbeds in this pebble layer (Katasonov and Ivanov, 1973). The middle Pleistocene environment during MIS 6, 170–130 ka BP, was characterized by a cold and severe climate in southern and northern Siberia (Chlachula, 2003; Andreev et al., 2004). The lacustrine loamy deposits which overlay the alluvial sands yield the U/Th date of 176 ± 2 ka BP (Katasonov and Ivanov, 1973). Ice Complex deposits of Central Yakutia accumulated during the MIS 4-2, a prolonged cryochron comprising the Zyryan stadial (MIS 4), the Kargin interstadial (MIS 3) and the Sartan stadial (MIS 2). Despite numerous detailed studies and the abundance of published radiocarbon dates (Tananaev, 2021), the age and origin of the Mamontova Gora Ice Complex deposits are still debated. Radiocarbon dates from host deposits enclosing the ice range from 36.7 cal ka BP (IM-155) to 47.4 cal ka BP (SI-1972) at depths from 3 to 8.8 m (Péwé et al., 1977; Kostyukevich et al., 1984). A series of AMS ¹⁴C dates from organic matter dispersed in wedge ice of “Ice Complex II,” from 14.9 to 21.9 cal ka BP, suggest their epigenetic origin (Vasil’chuk et al., 2004). Palaeoclimatic data obtained from the Central Yakutian Ice Complex evidence severe winters, and overall cold and dry environment for time intervals around 41, 21, and 13 cal ka BP (Popp et al., 2006).

The upper IC horizon, “Ice Complex I,” remains poorly studied, because it was only rarely found exposed in the thermal denudation scars of the Mamontova Gora. The age and genesis of these deposits remain unclear. On one hand, according to visual descriptions, host deposits of the upper IC horizon are, from their appearance lenses of lacustrine loams of the lower IC horizon subject to thawing during the Holocene and subsequently frozen with the development of ice wedges of epigenetic origin. The AMS ¹⁴C date 43.3 cal ka BP from the host loams, and similar stable water isotopic composition to the lower IC horizons, oppose this interpretation.

The age of Ice Complex at the Syrdakh site can be estimated from 16.6 cal ka BP (IM-433) to 23.6 cal ka BP (IM-433), from Ener and Syrdakh lakes, respectively (Tananaev, 2021). Dispersed organic material from Syrdakh ice yields AMS ¹⁴C date of 21.7 cal ka BP (KIA 26367) at 2 m depth, contemporaneous to the enclosing loams and supporting its syngenetic origin (Popp et al., 2006). Age distribution in an ice wedge in a nearby Ulakhan Syrdakh lake, 13.1 cal ka BP (KIA 26364) in the ice wedge margin and 3.8 cal ka BP (KIA 26365) in the ice wedge center, suggests that even if initial freezing and massive ice formation proceeded syngenetically, at least two later episodes of ice wedge development occurred in the region. These cold spells correspond respectively to Younger Dryas and Subboreal stage of the Holocene, signals of which are widely present in published regional paleoclimate proxies (Fradkina et al., 2005; Katamura et al., 2006; Nazarova et al., 2013).

The late Pleistocene-Holocene transition, between 11 and 9 ka BP, led to the widespread development of thermokarst lakes in Central Yakutia (Ulrich et al., 2019). The obtained (in this study) date of 10.6 cal ka BP years for the deposits overlying the ice complex is logical and associated with the accumulation of precipitation in the conditions of climatic warming of the Holocene. The Holocene climatic optimum in Central Yakutia falls on the period between 6.7 and 5.0 ka BP years, as evidenced by the data of the paleoarchives on diatom records (Pestryakova et al., 2012; Ulrich et al., 2017a).

5.1.2 Soil and Water Chemistry Evidence

Soil Carbon

The content of ТС is quite high, with a maximum in the loams of the upper horizon, reaching 6.6–6.8 wt% in peat interlayers (for unpeat interlayers, the value is about 4.0 wt%). In the lower horizon, these values are slightly lower, 3.4–4.0 wt%. High values of the content of organic matter are typical for the Yedoma sediments, the values in the active layer correspond to those previously described for Central Yakutia (Shepelev et al., 2016), but significantly higher than the data obtained for the Yukechi ice complex, which is very poor in terms of the content of organic matter (less than 1 wt%), dated 18–49 thousand years ago (Windirsch et al., 2020). Higher carbon contents in the upper horizon can be associated with more convenient conditions for the accumulation of organic matter, for example, due to the thawing of the lower IC horizon with the formation of lakes, the re-deposition of matter from other places, or with its better preservation. In Strauss et al. (2015), a higher total organic content was also observed in thermokarst deposits formed during the thawing of ice complexes. The available dating of 43.6 cal ka BP years, from the host deposits of the upper horizon, similar to the lower IC horizon (see above) may indicate the absence of new accumulation of organic matter at the time of the formation of the upper IC horizon (otherwise there would be younger dating). In this case, higher TC values are most likely associated with its better preservation. The more intensive decomposition or leaching of soil carbon could significantly reduce its content in the underlying layers. In the ice complex of the Syrdakh outcrop we observe lower values of carbon content, 2.1–2.2 wt%. The TC content is low in syngenetic ice complex deposits of Syrdakh, which were not subjected to melting and re-freezing.

Stable Water Isotope Evidence and Hydrochemical Evidence

Stable water isotopic composition of ice wedges is used to provide proxy-data for paleoclimatic and paleoenvironmental reconstructions (Rozanski et al., 1997; Sturm et al., 2010). In North-East Siberia, contemporary ice wedges are reported to have consistently heavier composition, getting significantly lighter toward late Pleistocene (Vasil’chuk, 1991; Vasil’chuk, 1992; Meyer et al., 2002; Wetterich et al., 2008; Meyer et al., 2010; Opel et al., 2011; Boereboom et al., 2013). Under warmer Holocene conditions and in different depositional environments, i.e., under riverine influence on flood plains, other mechanisms could have been involved, as suggested by the data from Mamontova Gora floodplain wedges (Opel et al., 2018).

Our stable water isotope data on Mamontova Gora wedge ice shows a pronounced difference between Holocene floodplain wedges and both Ice Complex horizons. Isotopic composition becomes increasingly lighter with depth and, therefore, age of the horizon. Water sources in ground ice can be deduced from both isotopic composition and basic water chemistry.

The slope of the regression line, in coordinates δ²H-δ¹⁸О, reflects not only the nature of water (atmospheric nature according to the GMWL correspondence, isotopically transformed water according to the displacement of values relative to the GLMW—this is mainly due to evaporation processes that reduce the slope of the line), but also ice formation processes in a closed or an open system. So, we get the highest ratio coefficient (8.2) for the lower generation wedge of Mamontova Gora which we interpret forming from winter snow, and a low slope coefficient (6.6) for the other profile from this wedge, which we can interpret as water source form surface or evaporated waters (seasonal thawed layer, or swamps with polygonal ponds). In principle, such low coefficients are also characteristic of ice formation in a closed system. The upper horizon of wedges is characterized by co-isotope slopes 7–8, which corresponds to the atmospheric origin of ice and insignificant participation of waters of non-atmospheric origin. The lowest values of the slope of the regression line, equal to 6.3 are obtained for Syrdakh ice, which probably reflect the relatively constant participation in the ice wedge formation of non-meteoric water. Ice segregation in a closed system, i.e., freezing lake basin (Fotiev, 2015), can also be involved in ice accumulation. High abundance of alkalis, Na⁺ and K⁺, shows important cryogenic metamorphization effect. Increased Cl⁻ content positively related with total dissolved solids suggests closed system conditions in the subsurface compartment. Otherwise, in open system chlorides are expected to leach from soils during the repetitive freeze/thaw cycles (Ivanov and Vlasov, 1974).

In general, high TDS content suggests dominant water origin of wedge ice from the subsurface compartment, while variations in major ion composition and ratios might relate to alternating ice segregation mechanisms or variable degree of cryogenic metamorphization. The latter leads to an increase in Na⁺ content, and alters the ratios between major ions, especially cations, and between alkaline Earth metals, Ca²⁺ and Mg²⁺ being the most affected (Ivanov and Vlasov, 1974). Point grouping in both the Na⁺ and Cl⁻-normalized molar ratios of Ca²⁺ and Mg²⁺ plots (Figure 6) suggest higher degree of cryogenic metamorphization in contemporary ice and several samples of the “Ice Complex II” wedge ice. High pH variability, along with extremely high Fe and DOC concentrations, above 40 and 53 mg L⁻¹, respectively, suggest an important impact of soil solutions in ground ice chemistry through the fixation of colloidal complexes produced through Al-Fe-humus release process (Pereverzev, 2009). High sodium carbonate content in the underlying lacustrine deposits, uncommon in loess layers, is seen by Kuznetsov (1976) as evidence for the subaerial loess origin, while for us it may also explain high Na content in the upper horizons of the “lower-generation” wedge ice through pore water filtration to the freezing front under epigenetic freezing.

Our results suggest, that at the Mamontova Gora exposure, the lower IC horizon, “Ice Complex II,” is a non-homogenous and polygenetic feature, where the isotopically lighter ice wedges were formed from snow meltwater with insignificant input from other atmospheric waters, while heavier ice wedges correspond to water origin from evaporated surface and shallow groundwater sources, i.e., lacustrine waters or soil moisture in the sub-lacustrine talik, interpolygonal ponds and shallow marshes. The isotopic composition of the Late Pleistocene ice wedges in Central Yakutia shows that snowmelt water was the dominant water source for ice wedge growth and yields the isotopic signal of winter precipitation. Surface waters, especially from alas lakes and shallow marshes, as suggested by their significant evaporative transformation, could also have a local effect on the isotopic composition of wedge ice.

5.2 Microbial Communities vs. Deposits

The bacterial communities of the studied horizons differ both in the composition of the dominant groups and in their ratio. Similarity patterns assessed using Bray-Curtis dissimilarity and phylogenetic (weighted UniFrac) metrics are not entirely identical at the class and phylum levels (Figures 10, 11; Supplementary Tables S2, S3).

FIGURE 10. Bray-Curtis dissimilarity heatmap for phyla (A) and classes (B). The samples AL, MMG and MMP correspond to previously studied samples of the ice complex of Mamontova Gora: an active layer (AL), an ice wedge of the upper (MMG) and lower horizons (MMP). For numerical values see Supplementary Tables S2, S3

FIGURE 11. Principal coordinate analysis (PCoA) scores plot of beta diversity based on UniFrac distance metric. MM for Mamontova Gora, S for Syrdakh. The samples AL, MMG and MMP correspond to previously studied samples of the ice complex of Mamontova Gora: an active layer (AL), an ice wedge of the upper (MMG) and lower horizons (MMP).

5.2.1 Neogene-Middle Pleistocene Alluvial Sands

The most ancient sediments, Miocene alluvial sands (С-3, С-5, С-6) yield one of the lowest Shannon Index values (3.0–3.3), and rather high similarity coefficients at the phylum level. The greatest similarity is between samples C-3 and C-6, belonging to the same VI terrace (dissimilarity 15%), sample C-5 of a V terrace is characterized by a lower similarity (dissimilarity 26–32%) with them. These sediments were deposited under warm and humid conditions of the Miocene and underwent epigenetic freezing only in the middle Pleistocene (see above, Section 5.1). According to the UniFrac distance metric, they are also combined with the syngenetic riverbed sands deposits (С-7, С-8) of V terrace. It is interesting that at the same time they show significant differences at the phyla and class levels, reaching a discrepancy of 39–57% with samples C-3 and C-5, and up to 78% with sample C-5. In middle Pleistocene sands (C-7, C-8), the Shannon Index values are lower (1.99–2.29), and less OTUs are allocated—this seems to be quite natural considering their synchronous formation and freezing. In addition, the samples of these sands have the highest similarity among the samples at the phylum level (dissimilarity 14%), which is explained by a fact that in fact these are the sub-samples from a single monolith. The profiles of the main communities remain the same for the Miocene alluvial sands and Pleistocene ones, but their ratio changes dramatically; in the sands formed in the much more severe climatic conditions of the middle Pleistocene (C-7, C-8), which were subjected to freezing synchronously with their accumulation, representatives of the phylum Firmicutes (with dominance of class Clostridia, genera Desulfosporosinus) begin to dominate (up to 43%), which either were completely absent, or were present in a minor amount in the Neogene sands formed in warm conditions (maximum 8%).

Published results on the Neogene sand community of a V Aldan River terrace (Brouchkov et al., 2017) show the same basic profiles at the phylum level, with a minor number of Firmicutes, the absence of Actinobacteritota, but with a predominance of Bacteroidota. In the Middle Pleistocene sands, Firmicutes phyla increases in abundance (up to 43%), while it was absent, or present in only a minor amount, in the Neogene sands formed under warmer conditions (maximum 8%). Actinobacteriota also remain in trace amounts, but the proportion of the phylum Bacteroidota slightly increases from Neogene to Middle Pleistocene sands. A large proportion of Betaproteobacteria, detected predominantly in samples C-3, C-5 and C-6, indicates the adaptability of communities to life at low temperatures and low nutrient content (Johnson et al., 2007). Actinobacteria and cold-loving obligate anaerobic, capsule-forming Bacteroidetes are also mainly distributed in strata of the Arctic permafrost (Jansson, Taş, 2014; Taş et al., 2018).

5.2.2 Ice Complexes: Ice Wedges

The studied samples from ice wedges of the Mamontova Gora exposure demonstrate significant differences in the community composition of the lower and upper Ice Complex horizons, further confirming that these are different generations of ice wedges with different origin (Figures 10, 11). We assume that the observed differences are related more to the genesis than the differential preservation of DNA since the microbial community structure within an epoch is relatively stable over time (Shade et al., 2013), while once the system surpasses a threshold, microbial parameters rapidly shift to a new stable state (Saidi-Mehrabad et al., 2020).

The group of samples from the upper “Ice Complex I” ice wedges (С-10 from this study, and ice sample of upper IC horizon (MMG) from Rakitin et al. (2020), are clustered together with Syrdakh ice wedge (led-S), based on the UniFrac distance metric, and are significantly distanced from the lower “Ice Complex II” (MMP in Rakitin et al., 2020). Bray-Curtis dissimilarity index at the phylum level between the lower IC ice wedges and all other samples varies from 0.52 to 0.68, while intra-group dissimilarity values are from 0.30 to 0.32. The maximum dissimilarity of 0.40 is between Syrdakh and previously studied MMG samples, mostly owing to closer position of the MMG sample to the active layer, allowing infiltration of organisms from the active layer (Rakitin et al., 2020).

These two groups of samples differ in both diversity metrics and compositional profiles. The ice of the lower “Ice Complex II” horizon yields a poor prokaryotic community, with the dominance of Firmicutes (54.5%) followed by Proteobacteria (31.4%). The communities of other named samples more diverse: the highest Shannon Index value, 5.17, is observed in the upper “Ice Complex I” horizon, it is slightly lower, 4.37, for the Syrdakh IC sample. At the phylum level, in addition to their greater diversity, the main difference between these groups is the emergence of photosynthetic anoxygenic gram-negative bacteria Chloroflexia in samples C-10 and led-S Syrdakh). This phylum can be associated with photosynthetically active aquatic ecosystem (11–21%), with prevailing Gitt-GS-136 and KD4-96 classes, which are most often found in soil or river sediments, freshwater lakes (Mehrshad et al., 2018), and were also widespread in taberal sediments, including recently thawed ones, together with Actinobacteriota (Gaiellales, 0319-7L-14) (Winkel et al., 2019).

We find a certain correspondence between the water origin ant the degree of evaporative transformation in ice wedges and the microbial community composition. Isotopically lighter ice wedges of the lower “Ice Complex II” horizon originating from snow meltwater yield also a poor, slightly diverse communities. Ice wedges at Syrdakh are same meltwater additionally fractionated by evaporation in surface or subsurface compartments. The Chloroflexia bacteria is absent in the first case but are abundant (21%) in the second. Minor evaporative transformation of waters in the upper ‘Ice Complex I’ horizon, corresponds with moderate presence of Chloroflexia (11%) in its bacterial community. This inter-relation can be driven by higher water temperature associated with higher evaporation, and with this, higher photosynthetic activity in the water later stored in ice wedges.

Besides, we observe the high similarity between the microbial communities of Syrdakh and the upper “Ice Complex I” horizon at Mamontova Gora, that may reflect their similar genesis and/or age. The enclosing sediments of the Syrdakh ice wedges have ¹⁴C dates between 23.6 and 16.6 cal ka BP. The similar calibrated age interval, from 21 to 14.9 cal ka BP, is attributed to ice wedges of the lower “Ice Complex II” horizon of the Mamontova Gora exposure (Vasil’chuk et al., 2004). Host deposits of this “Ice Complex II” horizon date to 47–35 cal ka BP.

5.2.3 Ice Complex: Host Deposits and Cover Layer

The highest similarity coefficients at the phyla level, Bray-Curtis index 0.22, were obtained for loams of the lower horizon, samples С-4 and С-11 (Figure 10; Supplementary Tables S2, S3). They are mainly represented by the same phyla, differing in the contribution of the phylum Firmicutes: 38% for С-4 vs. 21% for С-11. Starting from sample c11 and further up the vertical profile, the ratio remains in which Bacilli (Bacillus sp.) prevail over Clostridia (Clostridium sensu stricto 13) several times. Chloroflexi in С-11 (up to 3.9%), the share of the Thermoleophila class among Actinobacteriota is higher. On the Unifrac cluster plot, samples С-4 and С-11 are plotted far from each other, in contrast to С-4 vs MMP (Figure 11).

The Firmicutes species in the studied samples were observed in significant quantities only in syncryogenically frozen deposits: Middle Pleistocene sands, samples С-7 and С-8; late Pleistocene loams С-4 and С-11. But they are almost absent in the epigenetically frozen sediments formed in warm and humid Miocene conditions (samples С-3, С-5, С-6). This observation leads us to the suggestion that representatives of this phylum may be a marker of syngenetic freezing of deposits. The prevalence of the phylum Firmicutes, detected using the same primers and technique in clay permafrost layer from Western Spitsbergen, corresponding to the end of The Last Glacial Maximum, supports our hypothesis of Firmicutes as a marker of syngenetic sediment formation (Karaevskaya et al., 2021). Other studies of syngenetically frozen deposits exposed in the CRREL Permafrost Tunnel, Alaska, United States, showed that the content of Firmicutes OTUs, primarily from the spore-forming Clostridia and Bacilli classes, increased from an average relative abundance of 13% in the youngest age category (19 cal ka BP) to 79% in the oldest samples (33 cal ka BP) (40–60% for 27 cal ka BP) (Kanevskiy et al., 2008; Mackerlprang et al., 2017). Microbial community structure from the CRREL tunnel is close to that of our samples from the lower “Ice Complex II” horizon of Mamontova Gora.

The host loams of the upper “Ice Complex I” horizon of Mamontova Gora, sample C-9, have a high similarity coefficient to late Pleistocene loams of the lower IC layer, notably, the С-11 loam, located above С-4 sample in the vertical profile, at the dissimilarity level of phyla 0.32 and 0.42, respectively. UniFrac metrics, however, plots them close to the active layer sample even though the Bray-Curtis index at the phylum level, 0.58, is low (Figure 11; Supplementary Tables S2, S3). The reason for discrepancy between these coefficients is not very clear and requires further study. On the one hand, the significant differences in the microbial community composition between the active layer and the underlying frozen strata are typical (Müller et al., 2018), but on the other hand the close distance on Unifrac metris is logical too, since it deals with phylogenetic distances between observed organisms as well.

At the phylum level, Cloroflexi and Firmicutes are almost completely absent from the active layer but comprise up to 57% of the community in the underlying loams (sample С-9), which mostly differs from loams of lower horizon at the phylum level by the presence of the Chloroflexi phylum up to 30% (Gitt-GS-136 and KD4-96 classes). The fact that phototrophic bacteria are typical aquatic microorganisms, rarely occurring in dry soil, but developing very actively when flooded with water, confirms the hypothesis about the partial thawing of the ice complex deposits, followed by freezing. The fact that this thawing occurred in the Holocene is contradicted by the available dating and isotopic oxygen-deuterium results, which show similar data for the upper and lower ice (see above), in addition. We assume the IC of Mamontova Gora thawed in MIS 3 and was followed by refreezing of the deposits forming the new unit. The redeposition of material could take place there as it is shown in Kaplina (2011) for North Yakutia. This explains our dating in 43.3 cal ka BP of the upper IC horizon which fits in the range of data of the lower IC horizon. Therefore, the presence of the phylum Firmicutes in the upper loam is logical, since they are inherited communities of syngenetically frozen deposits of the lower horizon, In addition, a large amount of Chloroflexi (15–30%) can be noted in loams dated of 10.6 cal ka BP (our study), overlying the Syrdakh ice complex (С-12, С-12a), and the complete absence of Firmicutes there. If we adhere to the hypotheses about the formation of overlying loams due to the thawing of the IC (Vtyurin, 1975; Kaplina, 2011) during the Holocene optimum, then the presence of Chloroflexi (Gitt-GS-136 and KD4-96 classes) can be considered as a biotracer of this process indicating a likely watering of the territory, waterlogging due to excessive moisture, which commonly leads to thawing.

Freezing of loam deposits of the upper stage in the Holocene was, accordingly, epicryogenic, when the reconstructed temperature is characterized by values approximately equal to the present or slightly lower.

5.2.4 Deposits of Ice Complexes. Soils vs. Ice Wedges

High similarity at the phyla level of host loam and ice wedge ice of the lower generation, as well a similar composition of the main bacterial profiles at this level (dissimilarity 0.26) and clustering together the C-4 and MMP samples according to Unifrac, confirm the simultaneous formation and freezing of sediments and ice wedge. The C-11 sample, located closer to the upper part of the lower tier, that showed some distance from C-4 along Unifrac (and is closer to the sediments of upper tier of IC), but having a sufficiently high coefficient of similarity with the upper loam at the phylum level (0.32) can most likely be interpreted as undergone short-term thawing and immediately frozen.

As for the ratio of the upper ice wedge to the host sediments, the low dissimilarity (0.45) most likely indicates the non-simultaneous formation of loams relative to the upper ice wedge, as the ice wedge was formed epigenetically.

6 Conclusion

Our conclusions reflect on the hypotheses suggested in the Introduction:

1) The composition and structure of microbial communities are similar in host deposits and ice wedges under syngenetic sedimentation, and strongly differ under epigenetic formation of ice wedges.

High similarity coefficients were obtained between ice wedges and host sediments (Bray-Curtis dissimilarity 0.26 at the phylum level) of the lower IC unit of Mamontova Gora, while significant differences for observed for the wedges and host sediments of the upper IC generation (0.45), which is also confirmed by Beta diversity results based on the UniFrac distance metric, which considers phylogenetic distances between observed organisms in the computation.

2) The composition of microbial communities can be related to the water origin in the ice wedges and mark the IC strata subject to thawing and subsequent refreezing (epigenetic freezing).

We found a certain correspondence between the water origin and the degree of evaporative transformation in ice wedges and the microbial community composition. The Chloroflexia bacteria represented mostly by Gitt-GS-136, KD4-96 classes is absent in the lower generation ice wedges of Mamontova Gora originating from snow meltwater, and are abundant at the Syrdakh IC (21%), formed by the meltwater additionally fractionated by evaporation in surface or subsurface compartments. Minor evaporative transformation of waters in the upper “Ice Complex I” layer, corresponds with moderate abundance of Chloroflexia (11%) in its bacterial community. This inter-relation can be driven by higher water temperature associated with higher evaporative loss alongside with higher photosynthetic activity in the surficial water, later stored in ice wedges upon freezing.

For the soils, the abundance of the uncultured Chloroflexi of the Gitt-GS-136 and KD4-96 classes, is characteristic only to the layers that underwent thawing and refreezing. These bacteria were previously found both in freshwater bodies and in taberal sediments, including those recently thawed, and in our study, were abundantly present, up to 30% of the microbial community, in the loams from the upper generation IC Mamontova Gora, presumably thawed in MIS 3 interstadial, and in cover layer of Syrdakh IC, thawed during Holocene.

3) The similarity in microbial community composition and structure can be indicative of similar age of the IC strata in different locations, or their contemporaneous accumulation and freezing.

Different microbial communities of sediments and ice of the upper and lower generations of Mamontova Gora mark different stages of its formation. Based on the similarity of the communities for the upper IC layer of the Mamontova Gora and the Syrdakh IC, it can be concluded that these ice formations are somewhat synchronous. Still to confirm this we need data on the microbial communities of the host deposits of the Syrdakh IC.

In addition, the following interesting results were found that require further research:

1) Microbial diversity did not decrease monotonically down the profile according to sediment age. Thus, the lowest diversity with regard to soils was found in the Middle Pleistocene sands, and not Miocene sands. This may be caused by the deposition and long-term existence of the latter in a warm humid climate which creates favorable conditions for microbial diversity, rich species composition accumulated in the warm epoch was higher. Despite the destruction of cells later under harsh cold conditions, we suggest more DNA preserved from this stratum than from the later sediments formed in the more harsh conditions. These sediments were frozen only in the Middle Pleistocene, while the sands of the Middle Pleistocene were formed in much colder conditions synchronously with freezing (syngenetic freezing).

2) An abundant number of representatives of the phylum Firmicutes (the dominance of Clostridia is replaced by Bacilli at the beginning of the late Pleistocene) is characteristic only of sediments that were formed syngenetically: their accumulation and freezing took place at the same time in very harsh conditions.

3) The composition of microbial communities of the host syncryogenic loam deposits of the lower layer of the Mamontova Gora IC (dated to 47.4–36,7 cal ka BP) is similar to the composition of microbial communities (phylum level) in the syncryogenic deposits of the CRREL Permafrost Tunnel in Alaska, dated to 33–27 cal ka BP which makes such marker as microbial community similarity very promising to paleoreconstructions.

4) Absence of archaea in all samples earlier than late Pleistocene, which is probably due to their poor preservation.

Data Availability Statement

Metagenomic data uploaded to Genbank (https://www.ncbi.nlm.nih.gov) as bioproject number PRJNA734507 (SRX11067028-SRX11067041). The data on isotopes that support the findings of this study are openly available in PANGAEA at https://doi.org/10.1594/PANGAEA.921024.

Author Contributions

MC conceived the project and wrote the original draft of the manuscript. MC and DS participated in the field work and collected the data. SB and AYM performed DNA isolation. AYM supervised and provided bioinformatics analysis. AYM, EK, and MC performed microbiological interpretations. YV and NB performed isotopic analysis and interpretations. NT performed water and sediment chemistry interpretations. AR and AVM participated in bioinformatics analysis. MC, EK, AYM, and AB supported overall proxy data interpretation. All authors discussed the results and contributed to the final manuscript.

Funding

The reported study was funded by RFBR according to the research project 21-55-75004. The research of YV and NB was financially supported by the Russian Scientific Foundation (grant No. 19-17-00126, isotope analysis). The work of AR and AVM was supported by Ministry of Education and Science of the Russian Federation. The work of BS was carried out with the partial support from the Act 211 Government of the Russian Federation, Agreement no. 02.A03.21.0006.

Conflict of Interest

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

Publisher’s Note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

Acknowledgments

We thank Samsonova Vera, Karzhavin Vladimir, Pankov Alexander, Andreevskaya Maya and Alexander Osipov for their for their invaluable assistance in field work.

Supplementary Material

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

References