The sampling sites were chosen based on their grazing intensity, which was identified by long-term monitoring of animal preferences for grazing sites and additional observations over several days at the start of the field campaign in July 2019. We applied a space-for-time approach to compare apparently similar sites selected within the same landscape units that ideally only differed by grazing intensity, instead of long-term monitoring of non-grazed sites along with gradual herbivore introduction. For this pilot study with limited field time, we decided to prioritize the number of herbivore treatments and landscape positions investigated over the number of replicates within a given treatment and site in order to allow a preliminary analysis of a larger range of site conditions. We implemented a nomenclature scheme for our study sites that provides information on the landscape type (B for thermokarst basin or U for Yedoma upland) and the grazing intensity (3: intensive grazing; 2: occasional grazing; 1: no grazing). Five sites were selected to cover these different site characteristics. We sampled one soil core per site during our field campaign (Figure 2): B3 was an intensively grazed site in a wet area of a thermokarst basin; B2 was an extensively grazed site within the thermokarst basin, close to the fence of Pleistocene Park; B1 was a non-grazed site within the thermokarst basin, just outside the park’s fence; U3 was an intensively grazed site on a Yedoma upland; and U1 was a non-grazed site on a Yedoma upland. “Intensive grazing” and “occasional grazing” are relative terms in this context, with “intensive” referring to a daily presence of feeding animals over several hours throughout each day and “occasional” describing an occasional animal presence with feeding occurring along animal tracks. Location details are provided in Table 1.

In this setup, the area outside the fence was defined as an exclosure with no access for the park’s animals. We chose our sampling sites so that the three B sites were located in the same drained thermokarst basin, and the two U sites were located on the same Yedoma upland complex. The goal was to choose sites with similar properties so the animal influence was the main variable and differences between all other characteristics were minimized. Unfortunately, site B1 was experimentally deforested in 2015, prior to our study, which resulted in an unavoidable environmental difference that may have impacted soil composition at this site. The main impact is most likely soil compaction of the active layer due to the use of heavy machines. Also, flooding regime (seasonal or occasional) is different between our sites and might have some effect on the soils. Since animal trampling as well as defecation and browsing always occur where animals graze, in the following we will use the term “grazing” to indicate all animal activity including trampling, defecation, and foraging. Each site is exposed to the same grazing pressure year-round, as there are no separate summer and winter ranges, enabling us to study the net effect of the animals’ winter and summer impacts.

Firstly, we prepared a visual description of the surroundings at each sampling site, including the main vegetation type. We identified and sampled the most abundant vegetation species and estimated their coverage based on one by 1 m plots (Figure 3).

Secondly, we removed the active layer using a spade, until we hit permafrost. We measured the thaw depth and sampled the soil profile by using fixed-volume steel cylinders with a volume of 250 cm³. Due to very wet ground conditions, we were not able to collect cylinder samples at the intensively grazed site B3, where the ground was saturated from the surface to the frozen ground (38 cm). Instead, we cut blocks 8–10 cm high (see also Windirsch et al. (2021) for exact sample dimensions) out of the soil profile using a knife. The organic top layer was sampled separately at all sites.

Thirdly, we used a Snow Ice and Permafrost Research Establishment (SIPRE) permafrost auger with an inner diameter of 7.6 cm in order to sample both the still-frozen parts of the active layer and the underlying permafrost. We reached maximum sampling depths of 110 cm below surface (bs) at B3, 108 cm bs at B2, 127 cm bs at B1, and 114 cm bs at U3 (Supplementary Figures S2, S3). After drilling, soil samples and cores were individually wrapped in sterile plastic bags. All samples were brought in a frozen state to our laboratories for further analysis.

Due to high and dense shrub vegetation and therefore inaccessibility, we were not able to sample the U1 location during our summer field campaign. The Yedoma upland site outside the fence of Pleistocene Park was therefore sampled in the following winter with the SIPRE auger to a depth of 72 cm bs and was completely frozen at the time of extraction.

In the laboratory the frozen cores were cut into approximately 5 cm samples according to stratigraphy using a band saw. Afterwards, all samples were freeze-dried (Zirbus Sublimator 15) and weighed pre- and post-drying for ice and water content determination (Mettler Toledo KERN FCB 8K0.1, accuracy ±0.1 g). The dry samples were split into subsamples for biogeochemical and sedimentological analysis. In the following, the mean sample depth will be used to describe the position of each sample within the soil column.

The subsamples for biogeochemical analysis were homogenized using a planetary mill (Fritsch Pulverisette 5) and weighed into tin capsules and steel crucibles for measurement. We determined total carbon (TC), total nitrogen (TN), and total organic carbon (TOC) by combustion analysis using a vario EL III and a soliTOC cube (both Elemantar Analysensysteme). Afterwards, the carbon-nitrogen ratio (TOC/TN) was calculated from TOC and TN, giving information about the state of degradation and the source of the OM. TOC/TN ratios could not be calculated for samples with TN or TOC below the detection limit of 0.1 wt%.

OM for radiocarbon dating was taken from the dried original samples. For dating, the Mini Carbon Dating System (MICADAS) at the Alfred Wegener Institute Bremerhaven was used. We calculated the results in calibrated years before present (cal yr BP) using the calibration software Calib 8.2 and applying the IntCal20 calibration curve (Reimer et al., 2020; Stuiver et al., 2021).

We analyzed the ratio of stable carbon isotopes and used it as a proxy for the degree of decomposition of the OM following Diochon and Kellman (2008). Soil δ¹³C values differ between material sources and can distinguish between vegetation communities (Malone et al., 2018). Samples for δ¹³C analysis were homogenized using a planetary mill, and subsequently treated with hydrochloric acid at 50°C to remove carbonates. Measurements were done using a Delta V Advantage Isotope Ratio mass spectrometer (MS) supplement equipped with a Flash 2000 Organic Elemental Analyzer. Results are given in ‰ compared to the Vienna Pee Dee Belemnite (VPDB) standard (Coplen et al., 2006).

Finally, we determined grain size distribution within samples to check sediment substrate similarity. Therefore, OM was removed from the samples using hydrogen peroxide. Grain size distributions were subsequently measured using a MasterSizer 3000 (Malvern Panalytical).

We determined the annual average number of freezing days (daily mean air temperature below 0°C) over a period of 9 years. For this we used data from 2009 to 2017 measured at the Chersky meteorological station (station RSM00025123). Data were obtained from the National Oceanic and Atmospheric Administration’s National Climatic Data Center (NOAA NCDC) database.

In order to check for differences between intensively grazed and non-grazed sites, we combined our TOC data for intensive (B3 and U3) and non-grazed (B1 and U1) sites. We did this for the minimum active layer depth (38 cm at B3) to ensure similar conditions for all sites used. This left us with n = 11 samples for intensively grazed, n = 5 samples for occasionally grazed, and n = 10 samples for non-grazed sites. Since grazing intensities were artificially altered during the last few decades, we expect the most pronounced differences in the seasonally thawed layer, which corresponds to more recent periods. We visualized the distribution of TOC contents in the active layer in boxplots for intensive, occasional, and non-grazed sites of our study to explore possible differences in TOC content in relation to grazing intensity.

To identify correlations between TOC content, water or ice content, and sediment type (via mean grain size), we used principal component analysis (PCA). Data were initially normalized to values between 0 and 1. The PCA was conducted in the R environment using the “stats” package (RCoreTeam 2021).

At site B3 we found Calamagrostis langsdorfii as the only species in our one by 1 m plot with a coverage of more than 95%. No other species was found underneath. C. langsdorfii reached a height of approximately 30 cm. At B2 Poaceae were dominant but could not be identified on a species level due to a lack of blossoms as well as herbivory damages (Figure 3). Therefore, height could also not be recorded accurately. Poaceae covered approximately 70% of our plot with small amounts of Salix ssp. in between (less than 10% coverage). Close to our plot, Salix ssp. and Larix ssp. shrubs were present with a height of up to 2.5 m. Contrary to the grazed sites, the predominant species at B1 was Beckmannia syzigachne with an approximate coverage of 65% and an approximate height of 25 cm (Figure 3). Underneath, low extents of Saxifraga ssp., Equisetum ssp. and low-growing Salix ssp. were found. After deforestation in 2015, no Larix ssp. were present. U3 featured a much lower-growing vegetation (up to 15 cm height), holding approximately 50% of unidentified graminoids which could not be identified due to herbivory damages (Figure 3). In between we found Rubus ssp. (15–20%), Saxifraga ssp. (less than 5%), Salix ssp. (less than 5%), and some individuals of Vaccinium ssp. and Arnica ssp., which all lacked flowers and could therefore not be identified on species level. Towards the less intensively grazed upland areas, Salix ssp. shrubs up to 2.5 m in height grew with increasing density.

We found that the seasonal thaw depth (July 2019) decreased with grazing intensity, giving a thaw depth of 38 cm below surface (bs) in B3, 59 cm bs in B2, and 85 cm bs in B1 with a semi-frozen zone between 80 and 85 cm bs. In U3, we measured a thaw depth of 53 cm bs. Since U1 was sampled in winter, we were not able to measure thaw depth at this site.

At B3, TOC generally decreased from top to bottom with the highest value of 25.66 wt% at 20.25 cm bs and the lowest value of 1.18 wt% at 92.5 cm bs (Figure 4). In contrast to this, TOC values for B2 peaked around the frozen-unfrozen interface with values between 4.89 wt% (21.25 cm bs) and 30.10 wt% (64 cm bs); the frozen part contained generally more TOC. The highest TOC values among all sites were found in a peat layer at B1 in the frozen core part, with a peak value of 52.80 wt% (92.5 cm bs). The unfrozen core part contained much less OC with values between 0.79 wt% (35.25 cm bs) and 9.65 wt% (3.75 cm bs).

At the U3 upland site, TOC was higher in the unfrozen core part with values between 3.52 wt% (11.25 cm bs) and 10.73 wt% (49 cm bs). In the frozen core part below, TOC values were rather homogeneous, varying between 1.01 wt% (103.5 cm bs) and 2.65 wt% (58.5 cm bs). In contrast to the previous sites, TOC at U1 was homogeneous throughout the core at values between 1.24 wt% (33 cm bs) and 2.54 wt% (27 cm bs) with a slight increase at the bottom to 4.50 wt% (71 cm bs).

While TOC/TN ratios for B3 (range 12.68–25.96) and B2 (range 14.83–18.48) were similar (Figure 5), values for B1 were higher, in the range of 16.94–24.90. In the upland sites, there was a strong contrast between U3 (11.40–29.19, mean of 17.75) and U1 (28.63 and 29.27) in TOC/TN ratios.

δ¹³C values were also similar between the basin sites (B3: −30.63‰ vs. VPDB (uppermost sample) to −26.43‰; B2: −29.60‰ to −28.13‰; B1: −30.17‰ to −27.89‰) (Figure 4). Upland sites showed higher δ¹³C values in the range of −28.06‰ to −23.49‰. For full δ¹³C values, please see Supplementary Figure S4 or the published data set (Windirsch et al., 2021).

Radiocarbon dates were measured from samples at greater depth than the expected grazing influence. These data are provided in the Supplementary Table S1.

The soil cores were silty overall and had relatively homogeneous grain size distribution (Figure 6) across sites (Supplementary Figures S2, S3). In B3, clay content decreased slowly with depth from 16.69 vol% at 76 cm bs to 7.50 vol% at 106 cm bs, while at the same time sand content increased from 3.46 vol% to 10.68 vol%. Clay content varied throughout B2 with a minimum of 14.16 vol% at 28.75 cm bs and a localized maximum of 22.17 vol% at 13.75 cm bs (Figure 6). Sand content decreased with depth, starting at 12.73 vol% (5 cm bs) and reaching 2.43 vol% at the bottom (106 cm bs). In B1, clay content showed a peak value of 21.96 vol% at 76 cm bs, followed by a peat layer down to 106 cm bs with no grain size data available. Mean clay content for this core was 10.84 vol%. Sand content varied between 5.37 vol% (35.25 cm bs) and 11.31 vol% (43.75 cm bs) with two low values, 2.54 vol% at 52.75 cm bs and 2.22 vol% at 124 cm bs.

U3 showed similar characteristics as B3 with clay content decreasing with depth (highest value 19.76 vol% at 49 cm bs, lowest value 6.34 vol% at 112 cm bs), while simultaneously the sand content increased to a peak value of 11.52 vol% at 103.5 cm bs (Figure 6). Grain size composition for U1 showed no peaks, with clay contents between 8.03 vol% (71 cm bs) and 12.60 vol% (33 cm bs) and sand contents between 8.20 vol% (44.5 cm bs) and 14.01 vol% (71 cm bs).

The absolute water and/or ice content showed generally drier conditions in the frozen core of B3, with values fluctuating between 27.67 wt% (34.75 cm bs) and 78.59 wt% (1.25 cm bs). Similar characteristics were present in B2 (34.26 wt% at 28.75 cm bs to 71.10 wt% at 82.5 cm bs). In contrast, B1 showed large differences between the unfrozen upper part (17.19 wt% at 35.25 cm bs, 39.36 wt% at 3.75 cm bs) and the frozen lower part (39.59 wt% at 124 cm bs, 85.74 wt% at 112.5 cm bs).

U3 showed similar water content values at the top (32.00 wt% at 11.25 cm bs) and bottom (31.93 wt% at 112 cm bs) with higher values around the freezing interface (55.01 wt% at 49 cm bs) and lower values above and below that (22.11 wt% at 3.75 cm bs and 15.59 wt% at 87.5 cm bs). Water content in U1 was more stable, ranging between 25.63 wt% (39.5 cm bs) and 43.88 wt% (71 cm bs) with a similar value at the top (43.24 wt% at 2.5 cm bs).

The stable isotope characteristics of the pore water are shown in Supplementary Figure S5.

The PCA revealed a positive correlation between water/ice content and TOC content (Figure 7). While B3, B2, U3, and U1 tended to form clusters, B1 values were spread across parameters. For full PCA scores, see Supplementary Table S2 as well as Supplementary Figure S1.

In Figure 8 we visualized the range of the TOC values in the minimum active layer (depth 38 cm), grouped by grazing intensity (Figure 8). TOC is clearly higher under intensive grazing compared to non-grazed sites, with no overlap except for outliers. Due to potential dependency of TOC values within the unfrozen part of a core, and due to the small sample sizes, we did not test for statistical significance.

For our studied cores we found that intensively grazed sites (B3, U3) were covered by generally shorter and sparser vegetation (Figure 3, B3: grassland; Figure 3, U3: shrubby grassland tundra), with taller grazing-resistant individuals in between, compared to occasionally or non-grazed sites (Figure 3, B2: grasses and herbs, and B1: grasses, herbs and low shrubs). This could likely be a result of the reduction in shrub expansion through large herbivore action (Suominen and Olofsson 2000), while the differences in the flooding regime between the sites promoting different vegetation types might also have an impact on this. The grazed sites located in the thermokarst basin (B3 and B2) were still flooded seasonally. The upland sites were generally better drained than the sites within the basin. We aimed to select representative sampling locations, based on overall vegetation and observed animal routes. However, we had limited information on representativeness in terms of soil type and wetness across this landscape. Still, we ensured that our study sites are comparable in terms of soil properties by measuring grain size composition of the soil material (Figure 6). This revealed similar characteristics across sites with slight variations in clay and sand content, and silt contents ranging between 68 and 87 vol%, which shows that our sites are comparable in terms of physical soil properties.

Our results showed that the changes in vegetation height and structure correlated with large herbivore activity. This is in agreement with former studies in similar Arctic settings (Sundqvist et al., 2019; Skarin et al., 2020; Verma et al., 2020). However, the effects of seasonal flooding on vegetation composition and structure (sites B3 and B2) cannot be neglected. These differences in flooding regime marked a disparity between our sampling sites, since flooding at the basin sites B3 and B2 occurred seasonally or infrequently, while there are no records of any flooding for the other sites.

In our study area, shrubs were established on the upland before the introduction of large numbers of animals, and are now retreating, probably due to grazing pressure. While there are now higher shrubs at U3, there is a gradient in density of Salix ssp. shrubs towards less intensively grazed areas from individual shrubs to dense, forest-like coverage. At the edges of the drained thermokarst basin, shrubification took place on spots that could only be grazed occasionally (Olofsson et al., 2009) after the water level had retreated; the water retreated as the lake drained, allowing for incremental vegetation establishment. This leaves us with high soil wetness as another explanation for the absence of shrub vegetation (Martin et al., 2017). The clearing of the larch forest at the non-grazed site B1 4 years prior to our study is also a plausible explanation for the limited shrub extent there.

The restricted grazing space within the fenced Pleistocene Park area led to a higher revisitation rate of animals to sites within the fence. This artificially increased grazing pressure led to the different vegetation structures observed inside and outside of the fence. Such shorter graminoid vegetation types are linked to grazing intensity (Forbes 2006). Monteath et al. (2021) found that the extinction of the Pleistocene megafauna was not responsible for causing an increase in shrubification, but rather that shrubification as a consequence of climate change during the early Holocene could have led to the extinction of megafauna as a result of habitat loss. However, modern rewilding approaches are planned and realized with recent herbivores capable of living in today’s climate, though species composition might have to be adjusted continuously with warming conditions. Also, the cause and consequence in this setup are still unclear, with opposite findings across studies, depending on the methods applied (Murchie et al., 2021). Numerous studies have found that contemporary herbivores were indeed capable of reducing shrub extent in Arctic environments on a local scale (Olofsson et al., 2009; Köster et al., 2015; Ylänne et al., 2018; Verma et al., 2020; Mekonnen et al., 2021). We found shallow thaw depths at intensively grazed sites and a much thicker active layer at the non-grazed B1 site. Generally, the increase in active layer thickness is associated with warming of the ground, influenced by summer temperatures, flooding, insulation from snow, and vegetation density and composition (Walker et al., 2003; Skarin et al., 2020; Magnússon et al., 2022). Therefore, as was stated by Zimov et al. (2012), the direct influence on these insulating factors could affect the active layer dynamics due to a stronger cooling of the ground in winter (Figure 1). However, to determine any clear trends, repeated annual measurements of the thaw depth are needed.

Data from our basin sites showed that the more intensively grazed area featured smaller seasonal thaw depths in July 2019, compared to the occasionally or non-grazed sites (Figures 4, 9A). Keeping in mind the seasonal flooding at B3 and B2, one would expect greater thaw depths at these sites, since higher soil moisture promotes heat flow, and the relatively warm flood water adds heat to the soil as well (Magnússon et al., 2022). Examining the assumed permafrost table depth based on cryostratigraphic characteristics in the cores, we observed a contrary effect with an active layer depth of 61 cm bs under intensive grazing at B3, 74 cm bs under occasional grazing at B2, and 115 cm bs for the non-grazed B1 site. This supports our original hypothesis of grazing impacts being the main driver behind thaw depth differences, likely via vegetation changes.

At B3, the graminoid vegetation lay flat beneath the snow cover, offering little or no insulation. The snow itself was compacted by animal trampling, which further reduced the insulation properties of the snow cover. On the other hand, the forb undergrowth of the grassland at B2 provided slightly more resistance against the weight of the snow cover. Combined with less intensive trampling and, thus, a less compacted snow cover, the B2 sampling site was characterized by better winter insulation properties compared to B3. The Salix shrubs at B1 promoted localized snow trapping and a loose snowpack, both of which have been shown to provide effective insulation against cold winter air temperatures (Myers-Smith et al., 2011). Trampling was also missing in this non-grazed area. In addition to those surface conditions, the composition of the sediment itself could have contributed to a deeper active layer at site B1. From the grain size distribution data (Figure 6) we saw that material among all studied sites was quite homogeneous (Table 1). Exceptions were the sampling sites B1 and U1, which had a lower proportion of clay and were therefore generally dominated by coarser-grained material. High ice contents in the frozen part of the B1 core support the idea of an insulating effect of the peat layer that prevented the ground from thawing any deeper, but also helped to maintain the thick active layer found here. This assumption implies that the peat layer worked as a heat flux barrier for the permafrost below, stopping summer thaw from the top, but also preventing bottom-up refreezing in winter. This could lead to present or future talik formation above the peat layer.

We were not able to compare thaw depths for the upland sites: the cryostratigraphically determined permafrost table depth of approximately 50 cm bs at U3 was, however, smaller than the active layer depths measured by Abramov et al. (2019) and Shmelev et al. (2021) for the close-by Mt. Rodinka site, located on a Yedoma upland with a mean active layer depth of 80 cm. Also evident in this study (Abramov et al., 2019) were smaller thaw depths, around 40 cm, for thermokarst basin sites featuring shorter vegetation in the Kolyma region, which is similar to our B3 site.

In our study, the reduced active layer thickness observed at the intensively grazed sites could also have originated from active layer compression by animal trampling, which would mean that even with shallower thaw depth the same amount of soil thawed. However, this last point cannot be tested with our data, since we were not able to measure bulk densities due to wet ground conditions and root structures. While such ground compaction is very common in peatlands (Batey 2009), the sandy silt soils that we found are only moderately compressible (Akayuli et al., 2013; Halcomb and Sjostedt 2019), especially as they were frozen throughout a major time of the year, making soil compaction a weak argument for the shallow active layer at intensively grazed sites. The regular flooding of the intensively grazed basin site, which usually leads to an increase in thaw depth (Magnússon et al., 2022), combined with the in fact shallower active layer found in our study when comparing to the non-grazed B1 site, suggests that the cooling effects induced by animal activity could be much larger than previously expected.

We identified several variables that may be responsible for explaining differences in TOC between sites, including herbivory. Besides shorter graminoid vegetation types, other effects were also aligned with cooler temperature and/or less carbon degradation, such as soil moisture, as visualized by our PCA (Figure 7). This coexisting influence of grazing and soil moisture should be tested in other locations. However, the contrast in TOC/TN ratios and δ¹³C ratios (Figure 5) of the studied deposits between intensive and non-grazed locations suggested that animal presence is affecting OM decomposition (if we assume a comparable OM source). But based on the differences in vegetation cover, we found that this contrast might also originate from different OM input as well as from different nitrogen input via the presence or absence of animal droppings. Since grazing alters the vegetation cover—assuming that in the absence of large herbivores the vegetation would be more homogeneous throughout the study area—the observed strong differences in TOC/TN and δ¹³C ratios have likely resulted from the differences in animal impact on the landscape. Grazing is, therefore, a reasonable explanation for the strong differences between sites, including TOC content.

Mean TOC was six times higher in the top 38 cm of B3 (minimum active layer depth) compared to B1, with intermediate values in B2 (Figures 4, 8). For the upland sites, mean TOC was twice as high under intensive grazing in those core sections lying within minimum active layer depth (38 cm, B3) compared to non-grazed cores. The effect was not visible in lower core parts, perhaps because of the relatively short time span of 23 years since intensified animal introduction into the area (Figure 4). Instead, the TOC pattern in the frozen core parts was opposite to the higher TOC found in the active layer under intensive grazing influence. This suggests different pre-existing site conditions prior to herbivore introduction. The contrasting TOC pattern in the active layer suggests that animal grazing can help increase carbon storage in a relatively short time span (Figure 9B) but grazing needs more time to affect deeper permafrost carbon storage apart from thermal stability, hence reduced thaw depth. Therefore, the TOC differences between sites were not significant when comparing the entire depth of the sampled cores. We explain this as a result of i) the TOC-rich peat layer found in B1 that shifts the TOC median value for non-grazed sites upwards, combined with ii) low TOC values in the frozen part of the intensively grazed sites which offset the high active layer TOC. Since the peat layer was found at 85 cm bs and dated to 4,327 cal yr BP, we expect no connection to herbivory influence in a 23-year timespan, which is why we tested for significance in the active layer.

The significant differences in upper soil TOC between sites of intensive large herbivore presence (B3 and U3) and non-grazed sites (B1 and U1) agreed with the expected effects of grazing (Figure 9B).

All organic carbon found in the upland areas showed higher δ¹³C values (Figure 5), which we interpreted as a result of colder, drier, and more compact ground conditions in the Yedoma deposits and therefore less OM decomposition, as well as differences in plant species, resulting in different δ¹³C signals compared to the seasonally flooded and less compact ground of the drained thermokarst basin. Higher TOC/TN ratios in non-grazed sites indicated a different source material, originating from shrubby tundra, in comparison to grassland vegetation in occasionally and intensively grazed sites. However, these vegetation differences might originate from animal grazing impacts.

Above-ground carbon storage is also an important aspect of ecosystem carbon storage. Although we were unable to address above-ground carbon storage in this study, grazing activities of large herbivores have likely had an impact on above-ground biomass storage in vegetation by providing nitrogen fertilization via droppings. Mekonnen et al. (2021) stated that under intensified shrubification, shrubs take up large amounts of nutrients from the soil into the vegetation biomass. It is possible that CO₂ emissions from the soil might also be partially captured by this boosted plant growth and incorporated into plant biomass. However, nutrient availability is much more important for plant growth and hence carbon accumulation. Such nutrients become available when OM is decomposed, but the increased carbon accumulation by plant growth can to some extent offset the carbon losses associated with permafrost-stored OM decomposition (Keuper et al., 2012; Turetsky et al., 2020). While shrub growth increases above-ground biomass, it generally leads to net ecosystem carbon losses via heterotrophic respiration (Phillips and Wurzburger 2019). The formation of grasslands, on the other hand, increases above-ground carbon cycling as well, showing high productivity (Gao et al., 2016). While showing the highest respiration rates, soil carbon accumulation in wet Arctic grasslands is greatest when compared to other Arctic vegetation types (Bradley-Cook and Virginia 2018). This might be linked to the relative ease of grassland litter transfer into the soil, e.g. via animal trampling, as well as deeper root penetration of graminoids, compared to shrubby and sturdier vegetation (Wang et al., 2016). Based on this, we expect an increase in above-ground short-term carbon fixation in the more productive graminoid vegetation at intensively-grazed sites compared to non-grazed tundra sites, but cannot make any quantitative statements on this for our sites due to a lack of data regarding this matter.

We summarize two points from our carbon data:

(1) In the upper part of the soil cores, TOC contents were higher in more intensively grazed sites, with a decrease with less grazing for both thermokarst-affected and Yedoma sites (Figure 9B). As the soil and sediment characteristics of the non-grazed U1 Yedoma site matched findings from other Yedoma studies (Strauss et al., 2012; Jongejans et al., 2018; Windirsch et al., 2020), we consider this as the pre-grazing Holocene state (i.e., not heavily grazed state). This suggests that intensive grazing for 23 years has already increased TOC amounts stored in the active layer by a factor of two (Figure 4). At thermokarst sites, hydrology has an impact on TOC storage as well, which is seen in the strong correlation between TOC content and water content in the PCA (Figure 7). However, a direct comparison to a pre-grazing state for our sites is not possible.

Even if we cannot quantify the landscape and its carbon storage with no artificially introduced herbivores, we can make some assumptions based on other studies’ findings. These assumptions are hypothetical for our study area, but based on previously published studies. In a partially drained thermokarst basin within a Yedoma landscape, relatively warm soil conditions and hence a deep active layer would most likely be found in the drained basin (Windirsch et al., 2020). Normally, the active layer would first become shallower after lake drainage, with less water content leading to reduced heat flux into the soil and making room for cold winter air, and stabilize once tundra vegetation, shrubs, or even small trees re-establish, working as ground insulators in summer (Blok et al., 2010) in the absence of grazing. Over time, peat formation is likely to take place if wetness is still high, increasing carbon concentration and insulation effects on the underlying ground (Jones et al., 2012). In regularly flooded thermokarst basin areas, a graminoid vegetation holding vast extents of Eriophorum spp. will likely develop, while dry basin areas will normally feature Salix spp. shrubs (Regmi et al., 2012). On the surrounding uplands we expect a generally homogeneous appearance of soil carbon and vegetation, where differences are controlled by local parameters. These controls include differences in cryostratigraphy, substrate, and water availability, among others. On uplands we expect to find tundra vegetation, most likely with increasing shrubification throughout the last decades. Moreover we expect forested areas on those uplands with differing thaw depth between vegetation types and a thick but stable active layer beneath forests (Stuenzi et al., 2021). A deepening active layer is likely to occur underneath tundra as a result of global climate warming, which also increases shrubification (Wilcox et al., 2019) and therefore snow trapping. While shrubification may shield the ground from further summer thawing, the contrasting seasonal lengths of summer and winter may come into effect, reducing summer thaw but preventing winter refreezing at the same time. These literature-based qualitative scenarios suggest that the grazing impacts observed in our study, especially as regards alterations in vegetation and thaw depth, might represent a major system change; this may also be true for the increase in carbon storage at the intensively grazed U3 upland site. Hence, larger scale grazing could potentially help to keep carbon in the ground and reduce natural carbon emissions from permafrost, although in no case would this be a substitute for decarbonizing the current industrialized society. Due to technical constraints, the efficacy of this landscape management tool should not be overestimated.

Due to the spatial and temporal limitations in the data we collected, we cannot unambiguously distinguish between previous conditions at our sites, and the overall herbivory impact on soil TOC. Soil carbon distribution as well as local environmental differences influencing thaw depth might have existed before herbivore introduction; this cannot be accounted for in our study. However, the presence of rather fresh material present in the seasonal thaw layer at sites B3 and U3 (Supplementary Table S1) suggested a link to modern disturbances (Supplementary Table S1), for which animal trampling, bioturbation, or cryoturbation all provide a valid explanation.

Despite results being consistent with the expected effects of grazing on ground temperature proposed by Zimov et al. (2012), the missing core replication of the present study combined with the expected spatial heterogeneity of active layer depths indicated that this interpretation needs further data for validation. Further, more replicated studies, potentially combined with spatially comprehensive modeling exercises, are advised along with monitoring of exclosure sites with gradual herbivore introduction. Also a 23-year time period of herbivore activity is likely too short to affect deeper and permanently frozen deposits but mainly affects the active layer both in terms of chemical properties and thaw depth.

The larch forest clear cut at B1 marked a clear difference between our sampling sites, as stable boreal forest is known to protect permafrost against summer heat (Stuenzi et al., 2021), maintaining a deep but stable active layer. However, we still encountered the greatest seasonal thaw depth here. This greater thaw depth here agreed with the findings of relatively warmer but stable permafrost temperatures below forests in the Arctic and boreal region (Kropp et al., 2020) and was likely a relic of the previous forest-covered state. We found no sign of active layer refreezing. The thaw depth extended down to the peat layer, which started at 85 cm bs; the peat layer might have been functioning as a heat flux insulator both for summer thawing (from the top) and winter refreezing (from the bottom). Also, a detailed assessment of effects of different vegetation compositions typical for pastures in the Arctic on heat fluxes is needed.

Ultimately, more studies and monitoring setups on the impact of large herbivore grazing on permafrost soil carbon storage are needed across the Arctic. Hence, we hesitate to draw universal conclusions on herbivory effects from our data, including the extension of these results to much larger regions.