The Zackenberg Valley, Northeast Greenland (74°30′N, 20°30′W), is a wide lowland valley dominated by continuous permafrost. According to Elberling et al. (2010), annual air temperatures averaged −9.5°C between 1996 and 2007 and summer temperatures varied between 3 and 7°C from 1996 to 2007 (Hansen et al., 2008; Hollesen et al., 2011). The vegetation in the valley consists primarily of wet hummocky fens and low shrub and graminoid species (Elberling et al., 2008), while the lowlands east of the Zackenberg river are Cassiope tetragona and grass heathland (Bay, 1998). Permafrost up to 1 m deep had an average temperature of −2°C in summer and −14°C in winter between 1997 and 2006 (Christiansen et al., 2008). The average maximum active layer thickness ranged from 40 cm up to 2 m in depth and increased by 0.8–1.5 cm per year between 1996 and 2012 (Elberling et al., 2013). Long-term soil temperature data from 1995 to 2020 from the monitoring site SAL-1 at depths of 0, 20, 40, 60, 80, and 100 cm were retrieved from the Greenland Ecosystem Monitoring database^{1 ,2} and averaged for August and September.

After an intense snowmelt event in the 2018 summer season, a formerly existing minor thermokarst, a depression of the eroding permafrost soil surface, collapsed 400 m north-east of the Zackenberg Research Station. The incision into the soil profile was about 1 m deep and created an erosion gully that led toward the Zackenberg River. An ice lens at 40–60 cm depth was visible, which had melted in 2019 before sampling took place (Supplementary Figure 1). The thaw depth was measured with a thin steel pole pushed into the soil vertically three times until resistance indicated the frozen permafrost table had been reached (Christiansen, 1999). Based on these measurements, the layers were defined as the active layer (the long-term seasonally thawing top 40 cm), the transition zone of newly thawed permafrost material (40–70 cm in 2019 and 40–90 cm in 2020), and permafrost as the continuously frozen layer below these (below 70 cm in 2019 and 90 cm in 2020) as visualized in Supplementary Figure 2. In 2019 and 2020, three biological replicates were taken aseptically for every 10 cm interval until a depth of 90 or 100 cm, respectively. Due to the long transport chain without stable frozen transport technology and the bias, inconsistent freeze-thawing processes were shown to inflict on cold-adapted microbiomes (Lim et al., 2020), all samples were stored at 4°C at the Research Station before controlled, cooled transportation to Denmark, then stored at 4°C until further laboratory analysis took place as reported earlier (Gittel et al., 2014b).

For radiocarbon dating, field-wet soil per 10-cm horizon for 2019 was sifted in technical triplicated with a 0.5-mm sieve in deionized water to retrieve macro plant residues for further analysis, while the roots were removed under a stereomicroscope. To counter the low biomass, the triplicates were pooled per depth, before chemical pre-treatment with HCl and NaOH. The samples were graphitized before measuring the ¹⁴C isotope activity using an accelerator mass spectrometer (Radiocarbon Dating Laboratory, Lund University, Lund, Sweden). The obtained peaks were compared to established ¹⁴C calibration curves, considering anthropogenic atmospheric ¹⁴C activity changes. The results were reported as age in ¹⁴C years in BP (before present = AD 1950) from 1650 to 1950 using the terrestrial calibration curve IntCal13 (Reimer et al., 2013) and for younger datings than 1963 in fM (fraction modern) using the Levin post-Bomb calibration (Levin and Kromer, 2004).

Loss on ignition was performed in technical triplicate for each 10-cm depth interval per year on wet soil, by air-drying at 70°C for 48 h, then weighing to determine the weight-based relative soil water (H₂O) loss. The samples were burned at 450°C for 2 h in ceramic cups and weighed to obtain the weight-based relative organic carbon content (SOM). To measure the pH, 10 ml of air-dried soil in triplicates were added to 50 ml of 1 M KCl in Falcon Tubes in technical triplicate. After shaking for 1 h and resting for 1.5 h, the pH was measured with a Mettler Toledo FiveEasy Plus™ pH Meter (Mettler Toledo GmbH, Gieβen, Germany).

Based on the radiocarbon age and organic matter content of the soil, we defined each a surface (1) and buried (2) organic (O) and mineral (M) horizon, resulting in two organic (O1, O2) and mineral (M1, M2) horizons (Table 1 and Supplementary Figure 2). In contrast, layers were defined based on the thermal state of the soil with a local active layer (AL) designated as the seasonally thawing material until 40 cm. The transition zone (TZ) was designated as former permafrost material which had thawed since the collapse to 70 cm in 2019 and 90 cm in 2020. The underlying permafrost layer (PF) remained frozen, indicated by visible ice crystals in the material (Supplementary Figure 3).

DNA was extracted from 2 g field-wet soil per year and 10-cm depth interval for each biological replicate no later than 67 days after sampling, using the DNAeasy^® PowerSoil^® Kit (Qiagen, Hilden, Germany) according to the manufacturer’s instructions. Given the low biomass and DNA concentration, these extraction triplicates were pooled for sequencing in order to maximize coverage. The resulting extracts were quantified with a Qubit^® 2.0 Fluorometer (Thermo Fisher Scientific, Life Technologies, Roskilde, Denmark). For prokaryotic DNA, the 16S rRNA V3-V4 region was amplified with the primer set F341 (5′-CCTAYGGGRBGCASCAG-3′) and R806 (5′-GGACTACNNGGGTATCTAAT-3′) (Takahashi et al., 2014), while for fungal sequences the ITS2 region was targeted with the primers ITS1F2 (5′-GAACCWGCGGARGGATCA-3′) and ITS2 (5′-GCTGCGTTCTTCATCGATGC-3′) (Gaylarde et al., 2017). For the PCR, negative controls confirmed absence of contamination. The final PCR products were cleaned with magnetic beads (MagBio Genomics Inc., Gaithersburg, Maryland, US) and fresh 80% ethanol, and stored in 1x TE buffer at −20°C. The quality and size of the amplicons were verified by running an agarose gel and screening with Tapestation 4150 (Agilent Technologies, Santa Clara, California, US). Paired Illumina MiSeq 2 × 300 bp sequencing was performed at Teknologisk Institut, Taastrup, Denmark.

All of the following analyses were performed using the QIIME 2 pipeline and plugins version 2019.10.0 (Bolyen et al., 2019). Raw sequence reads were filtered by removing the multiplexing and gene region primer sequences with Cutadapt (Martin, 2011). The 50–60 cm 16S samples from 2020 and 0–10 cm and 60–70 cm ITS sequences from 2019 were removed due to an insufficient quality of reads after filtering. Noisy and redundant sequences were removed, to reduce sequencing errors and chimeras with dada2 denoise (Callahan et al., 2017). Unique filtered sequences were defined as amplicon sequencing variants (ASVs) and used for downstream analysis (Supplementary Table 1).

The prokaryotic ASVs were classified with the closest taxonomic affiliation with the classifier sklearn (Pedregosa et al., 2011) against the SILVA database (version 138_99) and the UNITE database version 8.2 (version 2020-02-20, Quast et al., 2013) with global and 97% singletons for fungi (Nilsson et al., 2019). The phylogeny was created with fasttree on classified sequences (Version 2.1.10, Price et al., 2010). For each sample, alpha phylogenetic diversity (PD) was determined with Faith’s PD index. Several other diversity indices were calculated but not used for downstream analysis (Supplementary Table 1). Beta diversity, as the variation of diversity across several samples, was calculated as Bray–Curtis dissimilarities (BC) as a non-phylogenetic index, and weighted Unifrac as pairwise phylogenetic distances.

On ASV level, phylogenetic alpha diversity correlation was tested with environmental metadata for variation, including pairwise Kruskal–Wallis covariation (non-parametric Mann–Whitney U-test) in QIIME2. Permutational multivariate analysis of variance (PERMANOVA) on 999 permutations was performed using the “adonis” function to determine significant drivers of community composition. A Principal coordinate analysis (PCoA) was performed on the weighted Unifrac distances and plotted as biplots, utilizing the Emperor QIIME2 plugin. Additional analyses were conducted using the “vegan” packages (Oksanen et al., 2019) in the R studio environment and R version 4.1.2 (R Core Team, 2021; R Studio Team, 2021). Correlation-based indicator species analysis (Dufrêne and Legendre, 1997) was performed on phyla abundances for each categorical parameter using the “indval” function. Non-Metric Multidimensional Scaling (NMDS) biplots with fitted environmental parameters were performed on BC dissimilarities based on the relative abundance of phyla.

At our study site, the weight-based relative soil organic matter (SOM) content indicated a mainly organic top 10 cm horizon (O1). It was underlain by the silt and clay mineral horizon (M1), which reached until the active layer limit at 30–40 cm depth with an average SOM content < 10%. Both were dated to be younger than 58 years old. Formerly separated from these by ground ice, the 20-cm deep organic horizon (O2) from 50–60 cm depth indicated up to 3,841 years old underlying material and after the melting of the ice lens in 2019 highest relative, weight-based soil water content (H₂O) of 28.27%. The deepest horizon (M2) from 60 cm depth and deeper consisted of sand and pebbles and up to 26,571-year-old material. This horizon included the upper permafrost table at 70 cm in 2019 and 90 cm in 2020. In 2019, the pH ranged from 3.87 of the top 10 cm to 7.01 at 70–80 cm depth, while the pH range in 2020 was stable throughout the profile from 4.01 to 4.91 (Table 1).

Amplicon sequencing revealed a frequency of 10,256,910 counts across 18 samples for raw 16S sequences. After filtering, 25,471 unique reads were defined as prokaryotic amplicon sequence variants (ASVs, for sequencing statistics see Supplementary Table 1). Overall, four ASVs remained taxonomically unclassified and 75 ASVs were archaeal (0.19% relative abundance). Ubiquitous phyla included Proteobacteria (29%), Acidobacteriota (19%), Verrucomicrobiota (13%), Bacteroidota (11%), Actinobacteriota (8%), Chloroflexi (6%), Gemmatimonadota (4%), and Patescibacteria (3%), as well as rarer taxa such as Firmicutes (0.7%), Bdellovibrionata (0.3%), and Cyanobacteria (0.1%), relative abundance across all samples indicated in brackets (Figure 1). While no single ASV was shared by all samples, the 11 most abundant ASVs with overall > 1% of all counts depicted 20% of total relative abundance (Figure 2). The overall most dominant identified genera belonged to the order Burkholderiales (Comamonadaceae family: Polaromonas and Rhodoferax, Gallionellaceae family: Sideroxydans, and Hydrogenophilaceae family: Thiobacillus), while Chthoniobacterales were dominated by Candidatus Udaeobacter, and Pyrinomonadales by RB41.

Specifically in permafrost samples (2019: 70–90 cm, 2020: 90–100 cm), Proteobacteria increased with depth, with the order Burkholderiales accounting for up to 60% of reads in the permafrost samples. Within this order, the genus Polaromonas was the most abundant ASV in the deepest permafrost sample of 2019 (80–90 cm), where also Firmicutes mainly occurred. Rhodoferax had the highest counts in the deepest and only permafrost sample of 2020 (90–100 cm). Sideroxydans occurred most in permafrost samples in 2019 (60–90 cm). In 2019, Thiobacillus dominated the community most between 60 and 80 cm (Figure 2). Actinobacteriota abundance increased with depth but decreased from 2019 to 2020.

Within the freshly thawing permafrost material of the transition zone (2019: 40–70 cm, 2020: 40–90 cm), relative abundance of Bacteroidota increased between 2019 and 2020, driven mainly by the Sphingobacteriales order in 60–90 cm deep samples. The phyla MBNT15 and Patescibacteria also increased in the transition zone samples of 40–80 cm in 2019 and 60–80 cm in 2020, respectively. Verrucomicrobiota were present throughout the thawed samples of active layer and transition zone. In both years, the genus Cand. Udaeobacter was found both in the thawed mineral top and buried organic horizons until 60 cm depth (M1 and O2). Here, Acidobacteriota were also abundant and most represented by the orders 11–24 and Pyrinomonadales, which were most abundant between 10 and 50 cm depth, decreasing toward 2020. For the latter order, the genus RB41 was among the highest counting ASVs in thawed soils in 2019, only seconded by Cand. Udaeobacter. Similarly, Chloroflexi also decreased in relative abundance with depth and time. Gemmatimonadota, Desulfobacterota, and Myxococcota abundance was rather constant throughout the samples, while the latter two had their maximal abundance in the surface organic horizon.

The indval analysis indicated the correlation of phyla with categorical environmental parameters (Supplementary Figure 6). Significantly more taxa were dominant in 2019, and only a few increased toward 2020, including Bacteroidota, Bdellovibrionota, and Fibrobacterota. These also correlated with increased abundance in the transition zone, together with MBNT15 and Patescibacteria, compared to other layers. While most taxa dominated the active layer (AL), Firmicutes indicated the strongest signal for significant abundance in permafrost samples. In contrast, each horizon had numerous relatively abundant taxa, without clearly outstanding trends (although Fibrobacterota and Firmicutes had their highest values in the deepest mineral horizon compared to other horizons; Supplementary Figure 6).

Archaeal relative abundance was marginal and even absent in permafrost samples. Phyla accounting for at least 1% of all archaeal counts included Crenarchaeota (59%), Nanoarchaeota (25%), Halobacterota (9%), Euryarchaeota (5%), and Thermoplasmatota (2% relative abundance within only archaeal ASVs across all samples). The two highest counting archaeal ASVs were assigned as Bathyarchaeia on 39% of all archaeal counts. The phyla Nanoarchaeota, Euryarchaeota, and Thermoplasmata each were driven mainly by one order, Woesearchaeales (25%), and the methanogenic Methanobacteriales (5%) and Methanomassiliicoccales (2%), respectively. Most of the Halobacterota abundance was also driven by other methanogenic taxa, such as the orders Methanomicrobiales and Methanosarcinales (7 and 2% total relative abundance).

Across 17 samples, raw ITS2 sequences had 1,543,256 counts, revealing 1,624 fungal ASVs after filtering, of which 63% remained unclassified. Only four phyla were classified as Ascomycota (21%), Basidiomycota (14%), Chytridiomycota (3%), and Kickxellomycota (< 0.005% relative abundance across all samples). No phylum appeared across all samples (Figure 1). For Basidiomycota, the genus Hebeloma was the most abundant in the first 20 cm in both years, accounting for 83% of the counts in the 10–20 cm deep sample in 2020 and 9% of all fungal counts. Ascomycota were particularly dominated between 0 and 50 cm depth by the genera Cystodendron, Blumeria, and Neobulgaria, accounting for 7% of all fungal counts (Figure 2). Most fungi within the permafrost samples remained unclassified.

Sequence statistics and alpha diversity indices were calculated and are summarized in Supplementary Table 1. As all alpha diversity indices indicated similar trends, we utilized Faith’s PD for all downstream analyses. The prokaryotic alpha diversity (overall Faith’s PD = 181 ± 91) consistently decreased with depth in 2019, while in 2020, very high diversity was maintained throughout the whole active layer and peaked in 40–50 cm deep samples before decreasing with depth (Figure 3). In contrast, fungal phylogenetic alpha diversity was comparatively lower (overall Faith’s PD = 54 ± 24), peaking in both years at 30–40 cm depth.

In our study, fungal diversity was generally lower than prokaryotic diversity and most driven by the difference of age and the thaw state of the soil. We could confirm representation of Ascomycota by Leotiomycetes and Basidiomycota by Agarimycetes in alignment with Gittel et al. (2014b), which might indicate a local pattern specific to Northeast Greenland, or even the Zackenberg valley. Ascomycota dominated our fungal community, in agreement with previous findings (Malard and Pearce, 2018; Margesin and Collins, 2019). The identified Ascomycete genera Blumeria, Cystodendron, and Neobulgaria of the Helotiales order are recognized as saprotrophic and plant pathogens, growing on decaying organic matter (Ekanayaka et al., 2019), as abundant in the buried organic horizon, where they occurred predominantly. Fungal abundance was commonly reported with decrease with depth in Arctic soils and is linked to a saprotrophic metabolism and the breakdown of labile carbon compounds (Gittel et al., 2014a).

Prokaryotic alpha diversity in our study was twice as high as reported in other studies (Chu et al., 2010; Tveit et al., 2013; Frank-Fahle et al., 2014; Gittel et al., 2014b; Schostag et al., 2019), which did not detect a deeper second peak, as we did at the 40–50 cm depth in 2020. Faith’s PD values both for prokaryotic and fungal communities were comparable to one of the few permafrost studies that previously employed this index (Feng et al., 2020). A decrease of alpha diversity with depth and age has been previously shown (Tripathi et al., 2018; Burkert et al., 2019). Interestingly, utilizing relic permafrost DNA, if amplifiable, was earlier found to lead to up to 45% increased diversity estimates (Fierer, 2017). Age was a significant driver of biodiversity for both kingdoms. The clustering of samples from the youngest samples opposed to higher dispersal in samples older than 3,500 years has also been described before (Tripathi et al., 2018), and aligns with the significant age effect on diversity (Mackelprang et al., 2017; Burkert et al., 2019).

For prokaryotic diversity, each layer depicted a significantly different community, but clearly only the most significant one, the oldest horizon, was visually split in the community, as indicated by the PCoA plots. For fungi, diversity changes were driven significantly similarly driven by the age effect and the transition to the permafrost layer, indicating a strong role for the permafrost table depth of the soil rather than organic or mineral content. Furthermore, current work has suggested that physical barriers to microbial dispersal were the main driver of changes in the permafrost soil microbiomes (Bottos et al., 2018). We think the separation of the community in the study by horizons most closely depicts these barriers of habitat, hence its significance for explaining changes in community composition.

Previous studies confirmed that prokaryotic survival in permafrost is due to the occupation of viable microhabitats, such as brine channels (Gilichinsky et al., 2003) or the ability to adherence to silt or clay particles, whereas the fungal mycelium often form larger networks across rather than within the soil horizons (Jansson and Hofmockel, 2019). Jansson and Hofmockel (2019) found bacteria usually outcompeted fungi as the temperature increased, due to higher growth rates and competition for newly available nutrients. However, for all microorganisms, our results indicate the strongest community shift from frozen to thawed ground.

While nutrient load and lability of organic carbon decreases within intact permafrost soils, the more stress-resilient oligotrophic (slow-growing specialist) microorganisms are usually favored over copiotrophic (fast-growing and generalist) taxa (Fierer et al., 2007). But these ecosystems are facing increased abrupt erosion until the end of this century, both globally (Turetsky et al., 2020) as well as in the Zackenberg Valley (Westermann et al., 2015; Rasmussen et al., 2018). Christensen et al. (2020) have already documented the connection between increasing events of extreme precipitation and summer temperatures, and the permafrost erosion site studied here. Habtewold et al. (2021) recently confirmed a correlation between higher carbon released from warmed Canadian soils and increased abundance of copiotrophic taxa. A former study on Alpine permafrost soil found soil warming to trigger copiotrophy (Perez-Mon et al., 2022), indicating elevated mineralization potential of increasingly bioavailable soil carbon and fresh nutrient inputs after thawing. The taxa that dominated this vulnerable zone of freshly thawed ancient carbon stocks in our study included both the copiotrophic Bacteroidia and MBNT15 phyla, as well as the oligotrophic Chthoniobacterales and Pyrinomonadales (Ho et al., 2017). But Bacteroidia were found as one of the clearest thaw zone indicator phyla and also dominated the 2020 thaw samples with more than half of all counts. They have higher carbon mineralization potential and can adapt their population size better to the available resources. In comparison, their oligotrophic counterpart Acidobacteria might be better adapted to the thermally variable active layer due to their higher stress resilience (Fierer et al., 2007). Interestingly, the abundance of copiotrophic taxa, such as the Comamonadaceae family and the often as copiotrophic classified Firmicutes even in our deepest permafrost samples (Ho et al., 2017), indicates that the short time period of storage at thawing temperatures initiated a microbial response to the abrupt thaw, both in situ and in vivo. This, clearly, indicates a need for controlled studies to verify and quantify the metabolic response of ancient permafrost taxa to thaw. Fierer et al. (2007) argued that knowledge about community composition and dominance of oligo- and copiotrophy is a crucial component for improving the prediction permafrost carbon loss due to microbial mineralization (Langille et al., 2013). Due to the quick response to readily available carbon, copiotrophic metabolism has been reviewed before to correlate with carbon release from soil ecosystems (Trivedi et al., 2013; Hurst, 2019).

Malard et al. (2019) and Tripathi et al. (2019) explored the description of stochastic (random and dispersal) and deterministic (biotic and abiotic selection) drivers of community composition in permafrost soils. They found top-soils to be dominated by deterministic processes, while deeper soils were controlled by stochastic processes. As permafrost soils are highly heterogenic, an assessment of these processes would help to confirm the significance of changes in this vulnerable ecosystem on its microbiome and should be considered in future efforts. Still, we can interpret from our community that higher abundance of specialist taxa in the permafrost samples indicates an importance of deterministic processes to govern, while these soils are intact.

In Arctic remote sampling sites, both core sampling techniques to reduce contamination (Bang-Andreasen et al., 2017) and temperature stable sampling and transport oppose great challenges. Due to the remote sampling site and long transport route, we decided to store the samples at 4°C and extract the DNA as fast as possible. Greater biological perturbation could have been induced with freezing and thawing. While DNA sequencing is a proxy for the total microbial community composition and not just the viable one—including taxa that might have ceased or increased abundance in response to the storage induced thaw—viability assays, incorporation, or transcriptomic studies are needed to quantify the copiotrophic taxa trends visible in our samples (Mackelprang et al., 2021). Still, taxa abundance in the deepest and most intact permafrost sample in 2019, as opposed to other depths and years, indicated the potential of depicting close to in situ conditions in ancient permafrost before and while thawing.