Introduction

Climate change is having a profound impact on Arctic biogeochemical cycles, hydrology, and geomorphology (Kokelj et al. 2021; Littlefair and Tank 2018; Zolkos et al. 2018; Zolkos and Tank 2019). Since 1965, Arctic average surface air temperature has risen about 3 °C exemplifying so called “Arctic Amplification” conditions, wherein the Arctic has experienced warming at a rate twice that of the global average since the 1990’s (Biskaborn et al. 2019; Richter-Menge et al. 2019; Schuur et al. 2008). Between 2007 and 2016, permafrost soils have warmed on average about 0.3 °C (Biskaborn et al. 2019), and it is estimated that up to 90% of permafrost that lies close to the surface will be thawed by 2100 (Gruber 2012). There is an estimated 1700 Pg of organic carbon (OC) held in permafrost region soils in the northern high-latitudes, and mobilization and degradation of this well-preserved OC has the potential to further increase atmospheric CO2 concentrations through a positive feedback cycle (Richter-Menge et al. 2019). Thus, understanding the fate of mobilized permafrost OC is of great relevance to the global carbon cycle due to the size of the permafrost OC store and the rapid pace of its warming and mobilization (Schuur et al. 2008; Report of the Intergovernmental Panel on Climate Change 2021).

About 15% of landmass in the Northern Hemisphere is underlain by permafrost, and about 16% of these permafrost regions are also covered by water features (Obu 2021; Vonk et al. 2015a). In ice-rich terrains, permafrost thaw causes ground subsidence following a reduction in soil or sediment volume associated with melting ground ice, referred to as “thermokarst”. In hillslope regions, thermokarst can cause the mass movement of materials downslope, exposing deeper soil layers and enabling transport into recipient aquatic networks. As water interacts with previously frozen soils, dissolved organic matter (DOM) can be mobilized to enter aquatic environments where it may be degraded via photochemical and microbial processing, or adsorbed onto mineral surfaces (Spencer et al. 2015; Vonk et al. 2015b; Ward et al. 2017; Olefeldt & Roulet 2012). The age of the mobilized dissolved organic carbon (DOC; the carbon component of DOM) can be variable based on its origin, and can range from Pleistocene-aged, to Holocene-aged, through to mobilization from the modern active layer (Lacelle et al. 2013; Drake et al. 2018; Rogers et al. 2021). The introduction of this ancient DOC into contemporary aquatic environments has ramifications for carbon cycling as permafrost thaw DOC has been shown to often be highly biolabile and to be rapidly utilized by microbial communities (Mann et al. 2012; Spencer et al. 2015; Vonk et al. 2015b). However, this input of DOC into aquatic environments may have a wide range of impacts and fates. For example, if mineral soils are also released along with DOC, the DOC can adhere to mineral surfaces and settle out with the sediments, effectively facilitating carbon burial and an overall reduction in DOC concentration (Collins et al. 1995).

Previous work in the Peel Plateau region in northwestern Canada has investigated overall DOC fluxes and biolability of thermokarst slump derived organic matter, as well as particulate organic carbon and inorganic chemistry (Littlefair et al. 2017; Littlefair and Tank 2018; Shakil et al. 2020; Zolkos et al. 2018). Thaw slump size, and the subsequent age of material exposed, are positively correlated with connectivity with downstream locations, and thaw slump size has been increasing over time (Kokelj et al. 2021). The overall impact of thaw slump material additions is dependent on the watershed size, with smaller watersheds experiencing greater impacts proportionally (Kokelj et al. 2021). Previous studies in the Peel Plateau showed Pleistocene aged glacial till deposits that were released from the largest slump sizes were associated with lower DOC concentrations at within-slump locations (Littlefair et al. 2017). This phenomenon was hypothesized to be related to the potential of rapid sorption of the fine-grained inorganic sediments from Pleistocene aged material, to the DOC also released from the retrogressive thaw slump formations. The biolability of thaw slump derived material was shown to be greater than that of unimpacted upstream samples, with less distinction among smaller thermokarst features (Littlefair and Tank 2018). Furthermore, DOM composition was the primary control of biolability in the Peel Plateau, supporting past studies from Northeastern Siberia, with permafrost sourced material shown to be highly biolabile (Drake et al. 2018; Littlefair & Tank 2018; Spencer et al. 2015). Ultimately, this permafrost sourced DOM has been shown to lead to significant CO2 outgassing in headwater streams at permafrost thaw impacted sites in Siberia (Drake et al. 2018), and so tracing its inputs and fate is of great importance for our understanding of the changing Arctic carbon cycle.

In previous studies centered on Arctic permafrost samples, a suite of analytical techniques have been used to try to source permafrost DOM inputs including optical parameters at these sites on the Peel Plateau (Littlefair and Tank 2018). Such techniques provide an insight into the bulk optically active portion of the DOM pool but do not have the resolving power of ultrahigh resolution mass spectrometry techniques which have shown the ability to distinguish DOM from a wide range of endmembers including anthropogenic (e.g. agricultural and urban) as well as terrestrial and aquatic sources (Kellerman et al. 2018; Kurek et al. 2021; Vaughn et al. 2021). Fourier transform-ion cyclotron resonance mass spectrometry (FT-ICR MS) has shown the ability, due to its high mass accuracy and resolution, to trace the unique fingerprint of permafrost derived DOM in past studies of yedoma inputs (Drake et al. 2018; Rogers et al. 2021; Spencer et al. 2015). In this study we set out to assess if FT-ICR MS could fingerprint thermokarst DOM inputs on the Peel Plateau. First, we aimed to resolve if permafrost derived DOM exhibits exclusive characteristics, as hypothesized, by examining the DOM composition between upstream and thaw slump derived samples to assess differences in bulk level metrics as well as unique molecular formulae present. Second, we hypothesized these unique characteristics could be tracked downstream and could be used in future studies as a tracer of thermokarst inputs on the Peel Plateau. Such an approach seeks to derive the ability to track future thermokarst DOM inputs by tracking the presence of specific molecular formulae within the DOM pool.

Materials and methods

General site description

The Peel Plateau is located in the continuous permafrost zone of northwestern Canada (Fig. 1). The region is underlain by ice-rich terrain that is comprised of unamended Pleistocene-aged glacial till deposited following the recession of the Laurentide Ice Sheet, overlain by previously thawed Holocene aged materials of several meters thickness into which permafrost was reestablished, and a thin, organic active layer at the ground surface (Duk-Rodkin 1992; Kokelj et al. 2017; Lacelle et al. 2013). The Richardson Mountain range lies to the west and produces a height of land defining the eastward drainage networks that have incised the ice-rich glacial deposits. The westernmost-reaches of the Plateau have exposed bedrock and colluvial deposits and transition through shrub-tundra and eventually subarctic taiga forest in the easternmost reaches (O’Neill et al. 2015). Warming in the early Holocene thawed the upper permafrost layer resulting in the development of Holocene-era soils, which were subsequently frozen as climate cooled and permafrost aggraded (Burn 1997). Warming air temperatures and increased seasonal precipitation has increased contemporary active layer thickness and fluvial erosion, accelerating the growth and development of retrogressive thaw slump features (Kokelj et al. 2015). These rapidly enlarging, chronic landslides can thaw permafrost up to tens of meters thick, and transport previously frozen materials downslope, reconfiguring slopes and fluvial networks (Kokelj et al. 2021). This process mobilizes organic matter into local aquatic ecosystems (Environment Canada 2015; Littlefair et al. 2017), with the contribution of active layer, Holocene-origin, and Pleistocene-origin permafrost dependent in part on thaw slump feature size and depth of permafrost thawed (Kokelj et al. 2021).

Fig. 1
figure 1

Location of thermokarst thaw slump sample sites on the Peel Plateau, Northwest Territories, Canada

Full size image

Thermokarst slump site selection

Six retrogressive thaw slump (RTS) features were chosen to represent various sizes of RTS present on the Peel Plateau. These sites have been the focus of a number of recent studies examining a range of geochemical impacts (Bröder et al. 2021; Kokelj et al. 2021; Littlefair et al. 2017; Littlefair and Tank 2018; Shakil et al. 2020). Such RTS features have grown in frequency and areal extent in recent decades in ice-rich permafrost regions and are chronic landslides which grow due to ablation of an ice-rich headwall (van der Sluijs 2023). The chosen slumps are categorized by previously-documented headwall sizes which range from: large size slumps represented by SE (23 m) and FM3 (10.3 m) (Zolkos et al. 2019); medium size slumps represented by sites HA (7.3 M), HB (7.1 m), and SF (7.6 m) (Zolkos & Tank 2019); and the smaller RTS feature category represented by site CB (5.8 m) (Bröder et al. 2021; Littlefair & Tank 2018; Shakil et al. 2020; Zolkos et al. 2018). Notably, these RTS features also occur across the east–west landscape gradient described above (Fig. 1). The RTS site names utilize the established naming scheme for these sites from past studies (Bröder et al. 2021; Kokelj et al. 2021; Littlefair et al. 2017; Littlefair and Tank 2018; Shakil et al. 2020). At each RTS site, water samples were collected at three locations: upstream of the RTS feature (i.e. control, before the influence of any slump inputs), within-slump runoff water, and downstream from the RTS feature, immediately below the point of all slump inflow.

Field sample collection

Slumps HA, HB, and SE were accessed via helicopter in late June (early season) and early August (late season). Slumps CB, FM3, and SF were accessed by foot, via a road access point from the Dempster Highway in early July (early season) and August (late season) 2021. These slumps were sampled, at three different locations (upstream, within-slump, and downstream), during two seasons (early and late season) except for HA and SF. Due to the inherent challenges of field sampling in Arctic conditions, slump HA was only sampled in the early season, and slump SF was only sampled in the late season. All water samples were collected in 1L acid washed polycarbonate (PC) bottles from just under the stream surface. Sample water was filtered with pre-rinsed 0.45 µm PES filters and the filtered water samples were immediately kept in the dark and subsequently frozen until analysis.

Dissolved organic carbon (DOC) and chromophoric dissolved organic matter analysis

A Shimadzu TOC-V analyzer was used to measure DOC concentrations of samples, and the mean of the best three out of five injections was used with a coefficient of variation < 2% (Littlefair et al. 2017; Murphy et al. 2013). Chromophoric dissolved organic matter (CDOM) absorbance was analyzed using a Biotek EPOCH2 spectrometer for UV–visible absorbance between 250 and 750 nm. Minor interferences were eliminated using a baseline correction procedure (Green & Blough 1994). The decadic absorbance at 254 nm (m−1) was divided by the DOC concentration (mg C L−1) to calculate the specific UV absorbance at 254 nm (SUVA254; L mg C−1 m−1). The SUVA254 value has been shown to be positively correlated with DOM aromaticity (Weishaar et al. 2003). a250:a365 ratios are inversely correlated with DOM molecular size (Peuravuori and Pihlaja 1997) and positively correlated with biolability (Olefeldt and Roulet 2012), and were calculated using established protocols comparing the ratio of the absorbance values between the 250 to 365 nm wavelengths (de Haan and de Boer 1987). The CDOM spectral slope coefficient between wavelengths 275 and 295 nm (S275−295; nm−1 × 10–3) was calculated using established procedures and has been related to molecular weight and aromaticity with higher values more indicative of smaller, less aromatic DOM (Helms et al. 2008).

Molecular level characterization of dissolved organic matter

DOC concentration values were used to determine the volume of sample required for solid phase extraction procedures as outlined in Dittmar et al. 2008. DOM extraction procedures were aimed to have a target concentration of 40 µg C mL−1 and were eluted with 1 mL of methanol. The extracted sample was frozen at − 20 °C before molecular level compositional analysis on a 21 Tesla FT-ICR MS at the National High Magnetic Field Laboratory in Tallahassee, Florida (USA) using negative electrospray ionization (Behnke et al. 2021; Hendrickson et al. 2015; Smith et al. 2018).

The molecular formulae determined by FT-ICR MS analysis were processed using a mass range of 130–1000 m/z and noise signals were excluded using > 6 σ root mean square baseline. Molecular formulae were then assigned in the PetroOrg Software (Corilo 2015). Elemental constrains of C4–45H4–92O1–25N0–4S0–2 (Kurek et al. 2022b) were used and calculated with no greater than 300 ppb mass accuracy (Behnke et al. 2021). Molecular formulae were divided into functional compound classes based on elemental compositions. The groups used included carbon, hydrogen, oxygen (CHO); nitrogen (CHON); sulfur (CHOS), and sulfur and nitrogen (CHONS). Condensed aromatics (0.67 < AImod), polyphenolics (0.50 < AImod < 0.67), highly unsaturated and phenolic low O/C (HUP Low O/C, AImod < 0.50, H/C < 1.5, O/C < 0.5), highly unsaturated and phenolic high O/C (HUP High O/C, AImod < 0.50, H/C < 1.5, 0.5 < O/C), and aliphatic (1.5 < H/C) were grouped based on established protocols (Kurek et al. 2022a, b; Šantl-Temkiv et al. 2013).

Statistical analysis

Base R was used with the ggplot2 package for data analysis and visualization (R Core Team 2019; Wickham 2016), in conjunction with Microsoft Excel, and plotted in MATLAB to create boxplots separated by sampling location to represent minimum, maximum, median, and mean values for DOM optical and compositional variables. A principal component analysis (PCA) was performed on the data using FT-ICR MS and optical derived parameters and visualized using the factoextra package in R (Kassambara and Mundt 2017).

Results and discussion

Permafrost slump impacts on DOC and CDOM parameters

Concentrations of DOC in upstream locations ranged from 0.85 to 12.65 mg L−1 (mean = 5.19 ± 4.21 mg L−1 std dev.; Table 1; Fig. 2a). In all of the box plots, each of the box and whisker shapes represent ten total analyzed samples. Of all upstream locations, the greatest DOC concentrations were found in the early season compared to the matched (i.e. same location) late season sample for all but the SE slump (Table 1). Within-slump DOC concentrations ranged from 1.45 to 22.85 mg L−1 (mean = 10.08 ± 6.89 mg L−1; Table 1; Fig. 2a). The within-slump locations exhibited both the greatest DOC range and mean values when compared to the other sampling locations. The greatest DOC concentrations were found in the late season in all of the within-slump matched pair locations besides slump SE. Downstream DOC concentrations ranged from 1.16 to 13.58 mg L−1 (mean = 5.52 ± 2.82 mg L−1; Table 1; Fig. 2a). Of all downstream locations, the highest DOC concentrations were observed in the early season matched pair locations (Table 1).

Table 1 Average and standard deviations for upstream, within-slump, and downstream permafrost slump locations for DOC and optical parameters for early season (ES), late season (LS), and combined averages (Avg)
Full size table
Fig. 2
figure 2

Box plots representing upstream (blue), within-slump (black), and downstream (red) slump location for (a) dissolved organic carbon concentration (DOC), (b) a250:a365, (c) SUVA254, and (d) S275-295. Solid circles represent means and horizontal lines represent medians

Full size image

Previous studies have found a broad range of results regarding DOC impacts from thermokarst features. In East Siberia (Spencer et al. 2015; Vonk et al. 2013) and Alaska (Abbott et al. 2014) slumps significantly increased DOC concentrations in waterways. However, within the Peel Plateau region, DOC concentrations have been previously shown to typically decline downstream of RTS impacted locations, even when within-slump values showed a modest increase in DOC concentrations (Littlefair and Tank 2018). This difference in DOC trends has been attributed to the sorption of RTS released DOC onto inorganic fine-grained, clay rich, glaciogenic sediments during the journey from the RTS complex into the riverine system (Littlefair et al. 2017). Overall, in previous studies in the Peel Plateau region, the DOC range recorded (5.4–26.1 mg L−1) was slightly higher than, but comparable to this study (0.85–22.85 mg L−1) (Littlefair et al. 2017). Variation between the studies is expected as some RTS sites were different than the previous study, which targeted more sites in the easterly, subarctic Taiga reaches of the plateau, and less in the westerly, mountainous reaches. Meteorological conditions on the sampling date also impact DOC concentration, with increased air temperature being positively correlated with DOC concentration (Littlefair et al. 2017).

Ratios of a250:a365 for upstream locations ranged from 4.27–6.66 (mean = 5.55 ± 0.83), and from 5.71–7.74 (mean = 6.50 ± 0.69) for the within-slump samples (Table 1; Fig. 2b). Downstream locations ranged from 5.22 to 8.28 (mean = 6.71 ± 1.05; Table 1; Fig. 2b), the largest range and average of all sample locations. At all sample locations, the lowest and highest a250:a365 measurements were in the early season (Table 1). CDOM characteristics like a250:a365 have been utilized to estimate DOM source and subsequent degradation, as well as to infer compositional properties of DOM pools (Helms et al. 2008; Peter et al. 2012; Spencer et al. 2007). As a250:a365 has been shown to be inversely correlated with molecular weight and directly correlated with biolability of DOM, this highlights the within-slump and downstream locations are expected to be more biolabile than upstream locations (Peuravuori & Pihlaja 1997; Olefeldt & Roulet 2012).

The range for SUVA254 in upstream locations was 0.94–6.75 L mg C−1 m−1 (Table 1; Fig. 2c). The within-slump locations had a smaller range from 0.50 to 3.06 L mg C−1 m−1, and the downstream range was even narrower at 0.84–2.48 L mg C−1 m−1 (Table 1; Fig. 2c). SUVA254 values were not substantially impacted by seasonality except in the upstream sample location, where the greatest SUVA254 values were found in the early season (Table 1). SUVA254 average values for upstream, within-slump, and downstream got progressively smaller (means = 2.63 ± 1.79, 1.88 ± 0.84, 1.57 ± 0.63 L mg C−1 m−1 respectively; Table 1; Fig. 2c). Previous studies of SUVA254 on the Peel Plateau showed the same relationship, wherein upstream locations had greater SUVA254 values indicating that they were more aromatic than within-slump locations (Littlefair and Tank 2018). This aligns with the expected biolability of downstream samples as evidenced by their a250:a365 values based on biolability of thawing permafrost trends seen in Eastern Siberian soils with lower MW being associated with greater biolability (Spencer et al. 2015; Vonk et al. 2013). Notably, the upstream SUVA254 values from the previous study were generally more than 1 L mg C−1 m−1 greater than in the current study (Littlefair and Tank 2018). Similar to for DOC concentration, this could be related to the less vegetated nature of the more westerly sample locations in this study (i.e. less vascular plant and litter inputs with elevated aromaticity available for leaching) (O’Neill et al. 2015).

S275-295 values ranged from 12.37 to 17.15 × 10–3 nm−1 in upstream locations, 16.10–19.96 × 10–3 nm−1 in within-slump locations, and 15.23–20.60 × 10–3 nm−1 in downstream locations (Table 1; Fig. 2d). Mean S275-295 peaked in the within-slump locations (mean = 18.17 ± 1.52 × 10–3 nm−1), which were steeper than the upstream locations (mean = 15.51 ± 1.42 × 10–3 nm−1), and was intermediate for the downstream locations (mean = 17.54 ± 1.92 × 10–3 nm−1). In both the upstream and the downstream locations, the shallowest values were measured in the early season, and the steepest values were measured in the late season (Table 1). There was limited variability with respect to seasonality in the within-slump locations (Table 1). Steeper S275-295 values have been linked with increasing internally produced and/or photodegraded DOM (Helms et al. 2008; Spencer et al. 2012), thus suggesting a less aromatic nature of RTS inputs at the within-slump locations versus upstream locations.

Molecular signatures of permafrost slump dissolved organic matter

Upstream locations had the lowest range of assigned molecular formulae (11,172–14,238) and lower total average number of formulae (mean = 12,620 ± 864; Table 2; Fig. 3a) compared to the within-slump and downstream locations. The within-slump locations had both the greatest range (12,879–18,219) and highest average number of assigned molecular formulae (mean = 14,650 ± 15,861; Table 2; Fig. 3a). Downstream locations exhibited an intermediate range and average (11,314–16,983 and mean = 13,431 ± 1701; Table 2; Fig. 3a). In both the upstream and the downstream locations, the lowest number of assigned formulae were in the early season and the highest number of assigned formulae were in the late season (Table 2). Recent work has highlighted that DOM from more diverse sources typically has a greater number of assigned molecular formulae (Kurek et al. 2023; Textor et al. 2019).

Table 2 Average and standard deviations for upstream, within-slump, and downstream permafrost slump locations for FT-ICR MS parameters, for early season (ES), late season (LS), and combined averages (Avg)
Full size table
Fig. 3
figure 3

Box plots representing upstream (blue), within-slump (black), and downstream (red) slump location for (a) number of formulae, (b) mean sample mass, (c) modified aromaticity index (AImod), and (d) nominal oxidation state of carbon (NOSC). Solid circles represent means and horizontal lines represent medians

Full size image

The greatest range in average mass was measured in the downstream locations (467.4236–567.2569 Da, mean = 520.8625 ± 28; Table 2; Fig. 3b), while upstream locations showed the narrowest range (516.5799–573.3658), but the greatest average mass (mean = 547.4840 ± 20) (Table 2; Fig. 3b). The within-slump mass average was the lowest at 510.8770 ± 25 (Table 2; Fig. 3b). These findings are in agreement with previous studies in the Kolyma River Basin, which found permafrost sources had a lower mass average than riverine samples (Rogers et al. 2021). In all sampling locations, the lowest mass averages were found in the early season and the greatest average mass values were found in the late season samples (Table 2).

The modified aromaticity index (AImod), is a parameter derived to highlight aromatic and condensed aromatic formulae in FT-ICR MS spectra (Dittmar et al. 2008), and has previously been shown to be inversely correlated with biolability of molecular formula in permafrost samples (Spencer et al. 2015; Vonk et al. 2013; Textor et al. 2019). The AImod was greatest for the upstream locations (0.31 ± 0.03) and the lowest in the within-slump locations (0.26 ± 0.03), with the downstream locations falling between the two (0.28 ± 0.03; Table 2; Fig. 3c). These values are similar to results from the Kolyma River Basin which found lower AImod values associated with permafrost samples when compared to recent vegetation and soil derived riverine samples (Rogers et al. 2021). Seasonality had a modest impact on AImod values, with early season locations showing a greater AImod in both upstream and downstream locations, and no effect of seasonality on within-slump average values (Table 2).

Average nominal oxidation state of carbon (NOSC) values were more negative in the within-slump locations (− 0.15 ± 0.06), indicating that slump-derived DOM was more reduced than DOM from upstream locations (− 0.06 ± 0.05; Table 2; Fig. 3d). A more negative NOSC value is also correlated with less potential energy derived from carbon oxidation per electron (LaRowe and van Cappellen 2011). Downstream locations had NOSC values similar to within-slump locations (− 0.14 ± 0.04), but once again fell in between within-slump and upstream location averages (Table 2, Fig. 3d). Late season samples had more negative NOSC values at all sample site locations (Table 2). This suggests that DOM from upstream locations was dominated by oxygenated terrestrial DOM, while DOM from within-slump locations represented a mixture of microbially-derived reduced DOM and degradation products from permafrost thaw inputs (Leewis et al. 2020).

Previous research has shown the relative abundance of CHO assignments comprise the majority of molecular formulae in DOM derived from permafrost soils (Rogers et al. 2021; Spencer et al. 2015). However, elevated relative abundances of CHOS, CHON, and CHONS formulae have also been found in permafrost samples compared to stream samples dominated by modern terrestrial inputs (Drake et al. 2018; Spencer et al. 2015). On the Peel Plateau, upstream sample locations showed enriched relative abundances of CHO molecular formulae compared to the within-slump locations, and the downstream locations were between the upstream and within-slump ranges (Table 2; Fig. 4a). Within-slump locations showed the greatest relative abundance of CHOS, CHON, and CHONS molecular formulae, while the upstream locations were relatively more depleted with respect to heteroatom containing formulae (Table 2; Fig. 4b–d). The downstream locations once again fell in between the upstream and within-slump values with respect to heteroatom containing molecular formulae (Table 2; Fig. 4b–d). Limited impact of seasonality was apparent with respect to heteroatom containing formulae throughout the sampling locations (Table 2).

Fig. 4
figure 4

Box plots representing upstream (blue), within-slump (black), and downstream (red) slump location for relative abundance of (a) CHO, (b) CHOS, (c) CHON, and (d) CHONS compositions. Solid circles represent means and horizontal lines represent medians

Full size image

Elemental ratios of H/C, S/C, and N/C were highest in the within-slump locations and lowest in the upstream locations while having the inverse trend for elemental O/C ratios (Table 2). H/C, S/C, and N/C in downstream locations fell in between upstream and within-slump locations, while O/C downstream average was the same as the within-slump average (Table 2). H/C, S/C, N/C and O/C did not show strong or consistent seasonality at any of the sampling locations (Table 2). Elemental N/C ratios are important in the Peel Plateau due to the apparent nitrogen limited condition of Arctic microbial DOM processing (Littlefair & Tank 2018; Sistla et al. 2012; Wild et al. 2015). Furthermore, greater nitrogen availability is correlated with increased microbial metabolic efficiency which corresponds to more energy used for growth and reduced CO2 released due to respiration (Mooshammer et al. 2014; Sinsabaugh et al. 2016). Therefore, nitrogen availability can directly impact carbon cycling in Arctic microbial ecosystems. Increased N/C ratios of thermokarst formations in the Peel Plateau is consistent with previous research in the region which found increased N/C ratios to be associated with greater biolability (Littlefair & Tank 2018). Elevated H/C and S/C ratios at the within-slump locations highlight the relative enrichment of hydrogen and sulfur-containing formulae commonly observed in aquatic environments draining from continuous permafrost (Drake et al. 2018; Kurek et al. 2023).

The compound class of polyphenolics has been associated with DOM from aromatic dominated terrestrial sources (Kellerman et al. 2018; Wagner et al. 2015). The greatest range and average for polyphenolic relative abundance was measured in the upstream locations (6.8–17.2%; 11.0 ± 3.2%), with declining range and averages from within-slump (7.5–13.3%; 9.7 ± 1.7%) to downstream (7.6–10.7%; 9.1 ± 1.0%; Table 2; Fig. 5a). Early season averages for polyphenolic relative abundances were greater than late season averages at all sampling locations (Table 2). The highly unsaturated and phenolic grouping with a high O/C ratio (HUP high O/C) showed the lowest averages in the relative abundance in the within-slump (31.4 ± 5.0%) locations and the greatest average in the upstream locations (39.7 ± 2.1%), and the downstream average fell in between (32.9 ± 4.1%; Table 2; Fig. 5b). The greatest relative abundance of HUP high O/C averages values were in the early season for upstream and downstream samples, while seasonality did not greatly impact within-slump averages (Table 2). The highly unsaturated and phenolic grouping with a low O/C ratio (HUP low O/C) showed the greatest relative abundance average in the downstream locations (44.1 ± 3.3%) with descending averages from within-slump to upstream locations (42.2 ± 3.2; 40.1 ± 4.8; Table 2; Fig. 5c). HUP low O/C averages were the greatest in the late season at all sampling locations (Table 2).

Fig. 5
figure 5

Box plots representing upstream (blue), within-slump (black), and downstream (red) slump location for (a) polyphenolics, (b) Highly unsaturated and phenolic with a high O/C (HUP, High O/C), (c) Highly unsaturated and phenolic with a low O/C (HUP, Low O/C), (d) condensed aromatics, and (e) aliphatics

Full size image

The relative abundance of condensed aromatics, indicators of highly aromatic terrestrial DOM (Kellerman et al. 2018; Wagner et al. 2015), were less defined by sample location, but the upstream locations had the greatest average (1.8 ± 0.8%), the downstream locations had the lowest average (1.5 ± 0.2%), and the within-slump locations fell in between (1.6 ± 0.5%; Table 2; Fig. 5d). The early season had the greatest condensed aromatic relative abundance averages at all sampling locations (Table 2; Schmidt et al. 2009; Sleighter and Hatcher 2008). The greatest average values of polyphenolic, HUP high O/C, and condensed aromatics in the unimpacted upstream locations thus highlights relatively greater terrestrial, more aromatic inputs at these sites. A similar trend has been shown in streams in Siberia between permafrost thaw dominated and modern terrestrial DOM dominated sites (Drake et al. 2018), as well as past research at these upstream locations has shown them to have greater modern terrestrial inputs (Littlefair et al. 2017).

The compound class of aliphatics has been linked to DOM from in situ production and microbial sources and is related to elevated biolability (Kurek et al. 2023; Spencer et al. 2015; Textor et al. 2019). Within-slump sample locations showed the most enriched average aliphatic content (15.2 ± 5.7), and upstream sample locations showed the most depleted average aliphatic content (7.4 ± 1.9), with the downstream values falling in between upstream and within-slump averages (12.4 ± 4.8; Table 2; Fig. 5e). There were limited impacts with respect to seasonality on aliphatic content for all sample locations (Table 2).

Tracing dissolved organic matter permafrost slump signals downstream

A PCA of FT-ICR MS and optically derived parameters resulted in two significant principal components, with PC1 explaining 58.4% of variance, and PC2 explaining 17.9% of variance within the dataset (Table 3 and Fig. 6). The main drivers that were positively associated with PC1 were AImod, HUP high O/C, NOSC, and CHO, and negatively aliphatics, CHONS, CHOS and S275-295; while HUP low O/C was the main positive driver for PC2 and condensed aromatics and polyphenolics the main negative drivers. Based on the PCA with the overlain sample points (Fig. 6b), it is readily apparent that PC1 is most strongly associated with sample location, with upstream locations clearly congregating positively on PC1, and within-slump locations plotting mostly negatively on this axis. Generally, downstream locations fall between upstream and within-slump locations from the same slump site and date. In contrast, PC2 appears to highlight the seasonality of the samples as a secondary driver, most notably for the upstream locations where early season samples plot more negatively on PC2. In general, westerly slumps with deeper headwalls (e.g., SE) plotted more negatively on PC1, while the sites located further east, or with shallower headwalls plotted more positively.

Table 3 PCA correlation matrix
Full size table
Fig. 6
figure 6

a Principal component analysis functions showing trends within the data, where the vectors convey the magnitude of their impacts. b The specific sample locations within the PCA space are mapped with circles representing early season and triangles representing late season samples. The blue corresponds to upstream samples, black with within-slump samples, and red with downstream sites

Full size image

The assigned molecular formulae across all sample locations were separated into formulae that were common to all upstream locations (upstream common), common to all within-slump locations (within-slump common), formulae found in all within-slump locations but not upstream locations (within-slump unique: WSU), and the within-slump unique formula that were also found in the downstream sites (within-slump unique in downstream: WSUD) (Fig. 7; Table 4; Supplemental Table 1). The WSUD formulae group represents the unique molecular fingerprint of permafrost RTS features that persists to the downstream sampling location after any instream in-situ processing including: microbial degradation, photodegradation, and abiotic processes such as sorption (Fig. 7d) (Collins et al. 1995; Kawahigashi et al. 2006; Littlefair and Tank 2018; Ward et al. 2017). Most notably, the WSU group had substantially greater N/C and S/C ratios, while the upstream common group had the lowest N/C and S/C ratios (Table 4). Transitioning into the WSUD, the CHONS formulae disappear from van Krevelen space when compared to the WSU group (Fig. 7c and d). This may coincide with the nitrogen-limited nature of microbial remineralization in the Peel Plateau, and thus within-slump nitrogen inputs into the Peel River could increase microbial respiration rates upon thaw (Littlefair and Tank 2018). Alternatively, these CHONS formulae may be preferentially sorbing during transport and they represent regions of van Krevelen space that readily adsorb (Subdiaga et al. 2020).

Fig. 7
figure 7

van Krevelen diagrams of a formulae present in all upstream samples (upstream common), b formulae present in all within-slump samples (within-slump common), c formulae found in all within-slump samples, but not in the upstream samples (within-slump unique: WSU), d the within-slump unique formulae that persist downstream (within-slump unique in downstream: WSUD), e Venn diagram of unique and common formulae for upstream, within-slump, and downstream locations; the number of formulae in each group is included

Full size image
Table 4 Unique and common formulae for upstream, within-slump, and downstream locations
Full size table

The WSU and the WSUD van Krevelen diagrams are more aliphatic relative to the upstream common group, exhibiting greater H/C average values of 1.23 and 1.22 respectively in comparison to 1.11 (Fig. 7; Table 4). The upstream common group has a greater relative proportion of high O/C aromatic formulae in comparison to the WSU and WSUD as evidenced by higher O/C values, higher AImod, and lower NOSC values (Fig. 7; Table 4; Supplemental Table 1). The WSU fingerprint is considered the permafrost thaw signature for the sampled sites as it encapsulates the molecular formulae that are solely found in RTS derived inputs. The difference between the WSU and the WSUD groups may be caused by sorption to minerals as previously hypothesized in the region (Littlefair and Tank 2018), or utilization by microorganisms (Spencer et al. 2015; Vonk et al. 2013) or a combination of these processes. The potential for photochemical degradation at these study sites to be a driver of changes in DOM composition between WSU and WSUD would appear to be limited due to high suspended sediment loads (Shakil et al. 2020) and no selective loss of aromatic moieties as evidenced by no overall change in AImod (Table 4).

The patterns of molecular level composition observed between locations in this study complement past research on the Peel Plateau that found greater biolability in within-slump locations (Littlefair and Tank 2018), as past studies have found greater biolability in relatively more aliphatic samples (Spencer et al. 2015; Textor et al. 2019). Future projections postulate a reduction in the extent of permafrost cover and a concurrent expansion of the active layer in soils in the Arctic with warming conditions (Lewkowicz and Way 2019; Schuur et al. 2008). As these conditions present themselves, an increase in thermokarst features including RTS can be expected to develop, and more within-slump conditions will be expected in the Peel Plateau (Lewkowicz and Way 2019; Kokelj et al. 2021). As thaw-slump derived DOM has been shown to have a greater prevalence of biolabile aliphatic formulae in this study, an increased relative contribution of aliphatic formulae to aquatic environments can be expected during conditions that favor the acceleration of thaw-driven mass wasting. Heterotrophic microbial degradation of aliphatic formulae has been shown to release greenhouse gases, however the balance between photosynthetic microbes and microbial respiration needs to be better quantified to assess the overall contribution to the global carbon cycle (Dutta and Dutta 2016). Furthermore, how these formulae will react when released into the downstream and coastal systems through riverine transport needs further study. Assessing the fate of this RTS DOM, it appears that some of this mobilized DOM may be sequestered through sorption to mineral surfaces and potentially buried in sedimentary deposits such as debris tongue or alluvial deposits (Littlefair et al. 2017). This compositionally distinct and energetically favorable DOM may also be rapidly utilized by microbial communities. However, to make these climatically relevant assessments both possibilities need to be examined further, and mechanistic sorption and bioincubation experiments undertaken incorporating ultrahigh resolution mass spectrometry to assess how DOM composition changes with these two mechanisms. Future studies linking WSUD formulae to biodegradable DOC would improve our understanding of how organic carbon from permafrost thaw is cycled through aquatic ecosystems and ultimately impacts the global carbon cycle. What is readily apparent from this study is that through utilizing FT-ICR MS it is possible to distinguish RTS DOM inputs into aquatic ecosystems in permafrost thaw-affected environments due to its unique molecular-level fingerprint, and that a portion of that unique fingerprint persists downstream. Utilizing this tracer will aid future studies seeking not just to track RTS permafrost inputs, but also examine how these inputs might be changing in a warming Arctic.