The evolution of stream dissolved organic matter composition following glacier retreat in coastal watersheds of southeast Alaska

Climate change is melting glaciers and altering watershed biogeochemistry across the globe, particularly in regions dominated by mountain glaciers, such as southeast Alaska. Glacier dominated watersheds exhibit distinct dissolved organic matter (DOM) characteristics compared to forested and vegetated watersheds. However, there is a paucity of information on how stream DOM composition changes as glaciers retreat and terrestrial ecosystem succession ensues. Importantly, it is unclear over what timescales these transformations occur. Here, we used bulk, isotopic and ultrahigh resolution molecular-level techniques to assess how streamwater DOM composition evolves in response to glacier retreat and subsequent terrestrial ecosystem succession. For this, water samples were collected from eleven streams across a chronosequence spanning a temporal gradient 0 to ~ 1400 years since glacier retreat in coastal, southeast Alaska. During the first ~ 200 years since glacier retreat, stream DOM showed marked and consistent changes in bulk, isotopic, and molecular-level composition. In particular, there was a decreased relative abundance (RA) of ancient, energy-rich (e.g., elevated aliphatic contribution), low aromaticity (e.g., low SUVA254 and AImod) DOM and an increased RA of soil and vegetation derived aromatic DOM (e.g., more depleted δ13C, elevated condensed aromatic and polyphenolic contribution) that had a modern radiocarbon age. After ~ 200 years of ecosystem development, DOM composition was comparable to that observed for other temperate and arctic forested watersheds without permafrost influence. These results underscore the timelines on which glacier retreat may have substantial impacts on watershed biogeochemistry and coastal ecosystems that receive DOM from these rapidly changing landscapes.

Abstract Climate change is melting glaciers and altering watershed biogeochemistry across the globe, particularly in regions dominated by mountain glaciers, such as southeast Alaska. Glacier dominated watersheds exhibit distinct dissolved organic matter (DOM) characteristics compared to forested and vegetated watersheds. However, there is a paucity of information on how stream DOM composition changes as glaciers retreat and terrestrial ecosystem succession ensues. Importantly, it is unclear over what timescales these transformations occur. Here, we used bulk, isotopic and ultrahigh resolution molecular-level techniques to assess how streamwater DOM composition evolves in response to glacier retreat and subsequent terrestrial ecosystem succession. For this, water samples were collected from eleven streams across a chronosequence spanning a temporal gradient 0 to * 1400 years since glacier retreat in coastal, southeast Alaska. During the first * 200 years since glacier retreat, stream DOM showed marked and consistent changes in bulk, isotopic, and molecularlevel composition. In particular, there was a decreased

Introduction
Climate change induced glacier retreat is likely to have major yet poorly constrained impacts on stream dissolved organic matter (DOM) composition, impacting the biogeochemistry of downstream and coastal ecosystems (Canadell et al. 2019;Hood et al. 2015;Lawson et al. 2014a;Milner et al. 2017;O'Neel et al. 2015;Pain et al. 2020). Currently, mountain glaciers provide a large seasonal flux of DOM downstream (0.58 ± 0.07 Tg C yr -1 ; Hood et al. 2015). Yet, many mountain glaciers, have, or are now beginning to retreat beyond the point of peak discharge (i.e., further loss in ice mass results in a declining glacier-derived contribution to streamflow; Milner et al. 2009;Vaughan et al. 2013). With continued mass loss projected to occur in the coming decades, the glacierderived DOM flux will diminish (Hood et al. 2015). These watershed-scale changes will impact the quantity and quality of streamwater DOM (Canadell et al. 2019;Fellman et al. 2014;Lafrenière and Sharp 2004;Pain et al. 2020), since the low elevation deglaciated portion of these catchments will become re-vegetated (Matthews 1999;Milner et al. 2007). Despite unabated glacier loss, the timescales on which contemporary retreat and subsequent terrestrial ecosystem succession alters streamwater DOM composition is not sufficiently understood to predict centennial changes to biogeochemical dynamics across the range of deglaciating watersheds.
Efforts to predict the impact of glacier loss on stream DOM composition have often used a space for time approach, comparing endmember watersheds with varying glacier cover to forested watersheds that have been deglaciated for millennia (Behnke et al. 2020;Fellman et al. 2010Fellman et al. , 2014Hood and Berner 2009;Lafrenière and Sharp 2004;Pain et al. 2020). Research highlights that glacier-derived DOM is chemically distinct from most aquatic ecosystems (Kellerman et al. 2021). Glacier-derived DOM is D 14 C deplete (i.e., aged; * 2000-10,000 ybp) yet highly bioavailable, associated with a high relative abundance (RA) of aliphatic compounds Arimitsu et al. 2018;Bhatia et al. 2010;Fellman et al. 2015a;Lawson et al. 2014b;Singer et al. 2012;Spencer et al. 2014a;Stubbins et al. 2012;Zhou et al. 2019a). Studies suggest that aged glacier DOM maybe derived from atmospheric deposition including anthropogenic aerosols, as well as potentially subglacial palaeoorganic material and microbial processes in supraglacial or subglacial environments (Bhatia et al. 2010(Bhatia et al. , 2013Fellman et al. 2015a;Lawson et al. 2014b;Musilova et al. 2017;Spencer et al. 2014b;Stubbins et al. 2012). Glacier DOM composition is the antipode to DOM in deglaciated, forested watersheds. This terrestrialderived DOM is typically less bioavailable, modern, has a high RA of aromatic and polyphenolic moieties, together with fluorescence properties indicative of fresh plant biomass and soil organic matter Behnke et al. 2020;Fellman et al. 2010;Kellerman et al. 2021;Lafrenière and Sharp 2004;Pain et al. 2020). Despite a comprehensive understanding of these endmembers and theoretical models proposed for the evolution of DOM composition with glacier retreat (Behnke et al. 2020), assessments of the nuanced changes in DOM composition between these end-states and quantification of the timescales associated with these changes are currently lacking. This is due to an absence of DOM compositional studies assessing watersheds that cover temporal gradients in glacier retreat. Therefore, it is unclear how DOM composition changes in recently deglaciated catchments through the processes of terrestrial ecosystem succession.
The Coastal Mountain Glaciers of southeast Alaska have experienced some of the highest rates of ice mass loss globally (Hugonnet et al. 2021;Zemp et al. 2019), with glacier runoff projected to diminish this century (Bliss et al. 2014). Glacier Bay, Alaska has acted as a model chronosequence study area for assessing how and over what timescales glacier retreat influences watershed ecology and biogeochemistry (Milner et al. 2007). The bay provides a unique opportunity as it has experienced rapid and well documented glacier recession since the end of the Little Ice age (ca. 1760; Motyka and Beget 1996), and thus has streams covering a well constrained and broad gradient in glacier retreat and terrestrial ecosystem succession. To date, studies have assessed changes to ecology, including above ground biomass, as well as soil type and quality, together with instream geochemistry including macronutrients, dissolved organic carbon (DOC) concentrations and optical proxies for DOM aromaticity (e.g., Bormann and Sidle 1990;Buma et al. 2017;Chapin et al. 1994;Engstrom et al. 2000;Milner et al. 2000Milner et al. , 2007Nagorski et al. 2014). However, it remains unclear how the chemical properties of DOM evolve in Glacier Bay in response to glacier retreat and subsequent terrestrial ecosystem succession, despite this assessment being central to understanding DOM reactivity and thus carbon dynamics in these rapidly changing watersheds.
The goal of this study was to investigate the centennial evolution of DOM composition following the onset of glacier retreat. We compare the DOM composition of streams in Glacier Bay, and near Juneau, Alaska, covering a stream age (time since glacier retreated from the stream mouth) gradient of * 0-1400 years (Milner et al. 2000). Water isotopes (d 18 O) were used to identify the sources of streamwater, particularly distinguishing between glacier-derived runoff, including snow/ice melt contributions, and rainfall. We analyzed DOM composition through bulk (DOC concentration and optical properties), isotopic (d 13 C-DOC and D 14 C-DOC), and molecular-level techniques (ultrahigh resolution Fourier transform ion cyclotron resonance mass spectrometry; FT-ICR MS). We draw on the full spectrum of data across the gradient in stream age to evaluate how DOM composition changes in response to glacier retreat and subsequent terrestrial ecosystem succession, and assess over what timescales these transformations occur. We hypothesis that over the centennial timeframe glacial contributions to streamflow will decline and streamwater DOM composition will shift from a glacier dominated, aged (D 14 C-DOC deplete) relatively aliphatic enriched DOM composition towards a modern, terrestrial-derived DOM signature (e.g., d 13 C-DOC deplete, relatively enhanced condensed aromatic and polyphenolic contribution).

Sample sites and catchment characteristics
Eleven streams in southeast Alaska were sampled from the 14th to 16th July 2013 during which time there was no rainfall. This period represents the peak glacier melt season where glacial contributions to streamflow are at a maximum. Ten streams are situated in Glacier Bay (Fig. 1). There is limited access to land terminating glacier outflows in Glacier Bay. Therefore, the eleventh stream, Herbert Glacier outflow (GO 0 ), which is a major outlet glacier of the Juneau icefield, is situated * 60 km southeast of Glacier Bay in the Juneau area ( Fig. 1). Both areas have a cool maritime climate, with average highs for May-September of 14.8 and 15.6°C for Glacier Bay and Juneau, respectively (https://www.ncdc.noaa.gov/).
Besides Herbert Glacier, stream ages are taken directly from Milner et al. (2000). To each stream age estimate, 16 years was added to account for the period between our sampling campaign and the Milner et al. (2000) sampling campaign (2013 and 1997, respectively). At Herbert Glacier, sampling was conducted within 200 m of the glacier terminus, hence its streamwater composition is considered a glacier endmember (i.e., stream age * 0 years). Watersheds in the northern part of Glacier Bay are in the early stages of terrestrial ecosystem succession (\ 100 years). Northern Glacier Bay watersheds have some portion of glacier cover or have been recently deglaciated, and have substantial barren ground and low vegetation cover (with potential for remnant ice and permanent snowfields; Nagorski et al. 2014). Older streams, to the south and near the mouth of Glacier Bay, are in the later stages of succession ([ 100 years) and have established above ground biomass and developed biotic systems, including wetlands (Milner et al. 2007;Nagorski et al. 2014). The oldest and most southern Glacier Bay catchment, Carolus River (* 1393 years), is an old growth forest watershed (of spruce, hemlock and cottonwood), which also has * 26% of its land area covered by wetlands (Nagorski et al. 2014). This stream is [ 1000 years older than the second oldest stream (Rush Point Creek), and its slope and elevation fall within the range of the other watersheds (Nagorski et al. 2014), thus it is used as a comparison to a mature forest watershed. Together the streams cover a gradient in stream age from 0 to 1393 years (GO 0 to Carolus River, respectively). At all sites, factors such as, presence of lakes, elevation, aspect, and slope may affect water quality (i.e., d 18 O), along with the speed of succession and vegetation community composition. All watersheds were grouped according to stream age: glacier outflow (Herbert Glacier) stream age 0 years (GO 0 ); early succession 50-100 years (ES 50-100 ); late succession [ 100-210 years (LS 100-210 ); and old growth forest 1393 years (Carolus River; OGF 1393 ). Since terrestrial ecosystem succession is a continuous process, statistical analysis was conducted on data without use of discrete categories (see Sect. ''Statistical analysis'').

Field sampling
Streamwater samples were collected using precleaned (10% HCl v/v for 48 h) 2 L NalgeneÒ bottles, rinsed three times with streamwater before collection. Collected water was immediately filtered through 0.7 lm pre-combusted (450°C, [ 5 h) GF/F filters and stored frozen in a series of pre-rinsed, acid washed (10% HCl v/v for 48 h) polycarbonate bottles (125 mL, 2 9 1 L), in the dark at -20°C. Water isotope samples were filtered using the same method, but stored in the dark at room temperature (20°C) in pre-rinsed 25 mL glass bottles with no head space (Fellman et al. 2014).

Water isotopes
Water isotope analysis (d 18 O) was conducted in order to provide an indication of the relative contributions of streamflow from glacier runoff (d 18 O deplete) and modern precipitation (snow and rain d 18 O enrichment compared to glacier runoff; Fellman et al. 2014). Samples were analyzed on a Picarro L2120-i analyzer and results normalised Vienna Standard Mean Ocean Water (VSMOW).

Dissolved organic carbon and carbon isotope analysis
Concentrations of DOC were measured using standard protocols on a Shimadzu TOC-V CHS analyzer (Stubbins and Dittmar 2012). DOC concentrations were calculated based on the average of 3-5 injections with a \ 2% coefficient of variance (Mann et al. 2012). For each stream, stable and radiogenic carbon isotopes (d 13 C-DOC and D 14 C-DOC, respectively) were analyzed. d 13 C-DOC was measured using a O.I. Analytical Model 1010 TOC interfaced to a PDZ Europa 20-20 isotope ratio mass spectrometer (Sercon Ltd.), as described previously ). D 14 C-DOC was measured using standard methods at the National Ocean Sciences Accelerator Mass Spectrometry Facility at the Woods Hole Oceanographic Institution, USA. Briefly, DOC was UV-oxidised and the resultant CO 2 was cryogenically purified and analyzed for D 14 C (Raymond et al. 2004).

Optical characterization
The fluorescence index (the ratio of emission intensity at 470-520 nm obtained at with an excitation wavelength of 370 nm) was calculated for each sample. This optical metric provides an indication of bulk DOM composition by its linkages to watershed sources of DOM (Cory et al. 2010;McKnight et al. 2001). Streamwater absorbance at 254 nm was measured on a Genesis spectrophotometer. Decadic absorbance was DOC-normalized to convert to specific ultraviolet absorbance at 254 nm (SUVA 254 ), an indicator of the aromaticity of bulk DOM (Weishaar et al. 2003).

Ultrahigh-resolution mass spectrometry
For each stream, sample water was acidified to pH 2 using reagent grade HCl, then extracted via solidphase extraction onto pre-cleaned 100 mg PPL cartridges (Agilent Technologies), using the procedure described in Dittmar et al. (2008). The volume of sample extracted was adjusted dependent on the DOC concentration of the sample, with a target loading of 40 lg C considering an assumed extraction efficiency. Cartridges were then rinsed twice with pH 2 ultrapure water and dried with ultra-high purity nitrogen gas. The dried cartridges were eluted with HPLC grade methanol into 2 mL pre-cleaned (10% HCl v/v 48 h) and combusted (550°C, 5 h) amber glass vials. The vials were stored at 4°C, in the dark, until analysis.
Eluted methanol was analyzed using electrospray ionization in negative mode on the 15 T FT-ICR MS (Bruker Solarix) at the University of Oldenburg, Germany. 500 scans were accumulated for each sample's mass spectra. Spectra were internally calibrated and molecular formulae were assigned to peaks with greater than five times the baseline signal to noise ratio using MolecularFormulaCalc-Ó, which assigns formulae within the bounds C 5-50 H 5-100 N 1-30-O 0-70 S 0-2 P 0-2 . The RA of each formula was calculated as the percent contribution of that peak to the total signal of all peaks within the sample. The nominal oxidation state of carbon (NOSC) was calculated for each formula (Riedel et al. 2012). Compound classes were assigned to molecular formulae based on stoichiometry and the modified aromaticity index (AI mod ; Koch and Dittmar 2006;Koch and Dittmar 2016). AI mod values of 0.5-0.67 and [ 0.67 were classed as polyphenolic and condensed aromatics, respectively; AI mod of \ 0.5 and H/C \ 1.5 were highly unsaturated and phenolic compounds (HUPs); H/C C 1.5-2.0, O/C B 0.9 and N = 0 were aliphatic; H/C C 1.5-2, O/C B 0.9 and N [ 0 were peptidelike; and H/C C 1.5-2 and O/C [ 0.9 were sugar-like (Koch and Dittmar 2006). Note each molecular formula may represent numerous structural isomers, since compound classes are based purely on stochiometric similarity (Hertkorn et al. 2007). The percent RA of a compound class was calculated as the sum total of the RA of all peaks belonging to that class.

Statistical analysis
Principal component (PC) analysis was used to connect the major DOM compositional changes between samples to stream age and thus watershed deglaciation. PC analysis was conducted using the FactoMineR package (Lê et al. 2008) and considered d 18 O, DOC, d 13 C-DOC, D 14 C-DOC and all optical and FT-ICR MS variables. FT-ICR MS metrics (NOSC, AI mod , and mass), together with compound and single heteroatom classes (CHO-only, CHON, CHOS and CHOP) were weighted according to RA. Metrics used in the PC analysis were scaled to one standard deviation, so that variables on different measurements scales were comparable. Asymptotic regression of PC1 against stream age was conducted using Base R ('SSasymp' function) and confidence intervals were calculated using nlstools (Baty et al. 2015). Spearman rank correlations, between the RA of each molecular formula and stream age, were calculated using the Base R 'cor.test' function. This was performed to isolate how an individual compounds RA changes with time since glacier retreat. Only significant correlations (p value [ 0.05) and formulae that were present in C 5 samples were included in the examination of spearman rank output. All statistical analyses were carried out using R (Version 1.1.463; R Core Team 2020) and figures were made using ggplot2 (Wickham 2016).

Results
Geochemical setting and bulk dissolved organic matter properties Streamwater d 18 O values ranged from -16.05 to -12.74% (Herbert Glacier outflow: GO 0 , to old growth forest: OGF 1393 , respectively; Table 1), reflecting a declining contribution of glacier or remnant ice and/or permanent snowfields to streamflow with increased stream age (Fig. 2a). Concentrations of DOC increased with stream age (range 0.35-1.19 mg C L -1 GO 0 to OGF 1393 ; Table 1, Fig. 2b). Our glacier and old growth forest stream (GO 0 and OGF 1393 , respectively) had d 18 O and DOC values that were within the ranges reported for their respective watershed types in southeast Alaska during the glacier melt season (i.e., June-September). For example, d 18 O values in southeast Alaskan glacier outflows and old growth forest streams ranged from -15.6 to -17.0% (n = 8) and -12.6 to -13.8% (n = 7), respectively (Fellman et al. 2014). In past studies, DOC concentrations ranged from 0.13 to 0.37 (n = 14) and 1.10-14.01 (mean 7.3 ± 4.6, n = 15) mg C L -1 for glacier outflow and old growth forest watersheds, respectively (Behnke et al. 2020;Fellman et al. 2014;Nagorski et al. 2014;Spencer et al. 2014b). The DOC concentration for OGF 1393 was at the lower end of previously reported values for old growth forest watersheds; this could be a function of soil properties, soil organic carbon stocks and/or the low streamflow during the sampling period, since transport of DOC to the stream would have been limited (Canadell et al. 2019;Hood et al. 2020). The marked similarity in d 18 O and DOC values of our endmember samples (GO 0 and OGF 1393 , where n = 1) with other similar watersheds in southeast Alaska, highlights that our watersheds and sampling period were representative.
Mean d 13 C-DOC values were inversely associated with stream age, with the most deplete d 13 C-DOC value observed at OGF 1393 (Table 1, Fig. 3a).  Fig. 3). The fluorescence index (FI) of DOM ranged from 1.39 (OGF 1393 ) to 1.64 (GO 0 ) and showed a progressive increase in the contribution of plant and soil derived DOM to the streamwater pool with increased stream age. Concomitantly, bulk aromaticity, as indicated by SUVA 254 , increased (1.69-4.62 L mg C -1 m -1 ) from GO 0 to OGF 1393 (Table 1, Fig. 3d). SUVA 254 values for OGF 1393 were comparable with those reported for other old growth forest watersheds in southeast Alaska Evolution of molecular-level dissolved organic matter composition with glacier loss GO 0 had the lowest number of assigned molecular formulae of all the streams (4838 compared to mean of 5830 ± 511 for the remaining streams; Table 1). The mean number of formulae increased from GO 0 to early succession (ES 50-100 ; mean # 5506 ± 365) and LS 100-210 (mean # 6053 ± 592) streams (Table 1). Typically, OGF 1393 had fewer formulae than LS 100-210 streams (# formulae 5802; Table 1). Across all streams, N, S or P containing formulae accounted for 32.8 ± 8.0% of the RA of all formulae. Heteroatom content declined from GO 0 to OGF 1393 , whilst CHOonly content increased (Table 1). The RA weighted average mass, AI mod , NOSC and the RA of aromatic and HUPs compound classes all increased from GO 0 to OGF 1393 (Table 1, Fig. 4). Concomitantly, the aliphatic and peptide-like contributions decreased (Table 1, Fig. 4). This was in broad agreement with previously reported values in southeast Alaska, where glacier outflows have a far greater RA of aliphatics and fewer condensed aromatic and polyphenolic compounds than old growth forest watersheds (Behnke et al. 2020;Stubbins et al. 2012). Sugar-like compounds contributed very little DOM (0.02-0.31% RA) to all streams and thus are not discussed further (Table 1). GO 0 had substantially higher aliphatic and peptide-like content (combined 49.7% RA) than all the other streams (mean 9.77 ± 5.21% RA), resulting in a marked decline in RA between GO 0 and ES 50-100 streams ( Table 1, Fig. 4).
PC analysis showed strong trends in molecularlevel DOM composition with stream age (Table 2, Fig. 5). There was clear division of DOM compositional variables in PC space (Fig. 5a). PC1 explains 70.3% of the variance in the dataset and was strongly negatively associated with the percent RA of aliphatic and peptide-like compounds, as well as FI and d 13 C-DOC (Table 2). These metrics clustered in PC space highlighting that bulk (optical and carbon isotopes) and molecular-level (FT-ICR MS variables) metrics were highly correlated (Fig. 5a). PC1 was positively associated with SUVA 254 and several FT-ICR MS variables including number of assigned formulae, NOSC, mean mass, AI mod , condensed aromatics (% RA), polyphenolics (% RA) and HUPs (% RA; Table 2). Furthermore, PC1 was positively associated with d 18 O, DOC and D 14 C-DOC (Table 2, Fig. 5a). Therefore, enriched d 18 O (i.e., precipitation rather than glacier-derived streamflow) was associated with more modern (D 14 C-DOC enriched) and higher concentrations of DOC, where DOM was higher in mass, oxidized (i.e., higher NOSC) and contained larger relative contributions of HUPs, aromatic and polyphenolic compounds (Fig. 5a).
Unlike the clear division of parameters in PC space (Fig. 5a), samples were distributed across PC1 reflecting the evolution in streamwater and DOM composition across the stream age gradient (Table 1; Fig. 5b). GO 0 and ES 50-100 streams were found at negative values of PC1, with LS 100-210 streams and OGF 1393 at positive values (Fig. 5b). The position of samples along PC1 was correlated with stream age and was represented by an asymptotic regression curve ( Fig. 5c). This highlighted that the rate of change in streamwater and DOM composition was initially rapid (i.e., within the first * 200 years) and that the rate of change decreased through time until a climax composition was reached (Fig. 5c). The change in DOM composition was characterized by a decline in RA of aliphatic and peptide-like compounds and an increase in RA of HUPs, polyphenolic and aromatic compounds (Table 1, Figs. 4,5). This trend was mirrored in bulk metrics (  (Figs. 3, 5). The 95% confidence interval of the regression curve indicated uncertainty in possible climax compositions and timing (Fig. 5c). It is unlikely that there is a single end-state composition across Glacier Bay, given the variability around the regression curve. However, the strong association of PC1 with stream age and clear trends with d 18 O mean that it is likely that watersheds will trend towards a comparable terrestrial vegetation dominated composition after * 200 years, as demonstrated by the minimal change in PC1 position after this point (Fig. 5c). Therefore, it is apparent that the major landscape control on DOM composition was time since glacier retreat, through its influence on catchment vegetation and soil organic carbon stocks.
Many of the individual molecular formulae present between the streams were highly correlated with stream age (# formulae correlated 1512; Fig. 6). Spearman rank correlations between the RA of molecular formulae (present in C 5 samples) and stream age showed both a change in RA and the detection of new molecular formulae with increased stream age (Table 1, Fig. 6). Higher mass, more oxidized (higher NOSC) condensed aromatic and polyphenolic formulae were enriched in the DOM pool and increased in RA in older streams (Table 1, Fig. 6b). Whereas, low mass, aliphatic formulae declined in RA and were eventually undetected with increased stream age (Fig. 6c). The correlation of molecular formulae across the streams showed a consistency in composition and that DOM Fig. 4 The RA weighted average mass, AI mod and NOSC (a-c), together with compound class scores (d-f) seperated by stream age category. GO 0 and OGF 1393 streams are a n of one. Early and late succession streams are abrievated to ES 50-100 and LS 100-210 , respectively. Dots represent individual datapoints composition evolves as stream age increased, trending towards an enrichment in aromaticity.
Despite marked trends in DOM composition with stream age, there was variability in the DOM composition within this trend (Fig. 5). This is reflected by the departure of streams along the PC2 axis, especially in LS 100-210 streams (Fig. 5b). PC2 explains 13.2% of the variability in the dataset and was positively associated with the RA of CHON formulae (   Fig. 1 for full name): GO (GO 0 , 0 years), ES 1-4 (stream age [ 50-100 years, early succession), LS 1-5 (stream age [ 100-210 years, late succession) and OGF Carolus River (1393 years). c Asymptotic regression for stream position on PC1 axis with stream age, grey shading represents the 95% confidence interval of the asymptotic regression parameters  variable positioning of streams along PC2 likely reflected minor compositional differences between these streams that are not driven by glacier retreat.

Discussion
Glacier retreat driven changes in dissolved organic matter composition Bulk optical, isotopic and molecular-level data showed that time since glacier retreat was a major control on streamwater DOM composition (Fig. 5).
The distinct ancient aliphatic composition of glacierderived DOM was lost from streams within the first * 200 years since glacier retreat, suggesting DOM contributions to streamwater from glacial sources (e.g., glacier ice, remnant ice and snowfields) declined with increased stream age (  Figs. 3, 4). The observed trends in the streamwater DOM composition with glacier retreat were consistent with previous studies that assessed watersheds with gradients in glacier cover, where low glacier cover was associated with more enriched bulk D 14 C-DOC values (Fellman et al. 2014;Zhou et al. 2019a) and DOM that contained a large proportion of organic matter (OM) from fresh plant biomass and soil organic material (i.e., d 13 C-DOC depletion and increased % RA of aromatics and HUPs). Despite the ancient aliphatic signature of DOM being a ubiquitous feature of glaciers and their environments, the sources of OM that explain glacier outflow DOM composition remain uncertain (Behnke et al. 2020;Bhatia et al. 2013;Kellerman et al. 2021;Lawson et al. 2014b;Singer et al. 2012;Spencer et al. 2014a;Stubbins et al. 2012;Zhou et al. 2019a). GO 0 had low aromaticity (i.e., low SUVA 254 , AI mod , and % RA of condensed aromatic and polyphenolics; Table 1, Figs. 3, 4), suggesting minimal DOM contribution from terrestrial plant and soil sources (Behnke et al. 2020;Koch and Dittmar 2006). Thus, significant contributions of OM from overridden, subglacial palaeoorganic material are unlikely. It is possible that microbial processing of overridden soils may contribute to the highly aliphatic, ancient signature of glacier outflows and thus GO 0 (Bhatia et al. 2010(Bhatia et al. , 2013Lawson et al. 2014b;Musilova et al. 2017). Similarly, microbial autotrophy in supraglacial or subglacial environments could contribute to the glacier DOM-pool (Bhatia et al. 2010;Musilova et al. 2017;Smith et al. 2017). Nonetheless, there was similarity between the DOC concentrations, as well as the DOM age and composition of GO 0 (together with DOM in supraglacial environments and glacier outflows globally) with water soluble organic carbon from atmospheric aerosols and DOM in precipitation (e.g., OM from fossil fuels, secondary aerosols and biomass burning; Fellman et al. 2015a;Gurganus et al. 2015;Spencer et al. 2014a;Stubbins et al. 2012;Wozniak et al. 2014). Therefore, our findings support the hypothesis that atmospheric deposition on glacier surfaces is a major source of organic material exported from the glacier outflow, particularly in southeast Alaska (Fellman et al. 2015a;Spencer et al. 2014b;Stubbins et al. 2012).
Early succession streams (50-100 years) had higher aliphatic contributions and FI values than LS 100-210 streams and OGF 1393 (Figs. 3, 4), as well as typical non-glacial southeast Alaskan riverine DOM (Behnke et al. 2020;Fellman et al. 2014;Stubbins et al. 2012). Relatively depleted d 18 O values (-15.6 to -14.7%) suggested that glacier-derived contributions of DOM (e.g., from remnant ice or permanent snowfields) remained somewhat important in ES 50-100 streams (Fig. 2). Similarly, atmospheric deposition is not restricted to solely glacier surfaces (Spencer et al. 2014a;Stubbins et al. 2012), thus may have supplied DOM, which was still readily apparent (i.e., not overprinted by terrestrial sources), in ES 50-100 streams with no glacial contributions to streamflow. More enriched D 14 C-DOC values in ES 50-100 streams than GO 0 could mean that autochthonous production of DOM was a relatively important source of OM a short time after glacier retreat, when riparian vegetation OM inputs were minimal (Milner et al. 2007). Microbial processes have been shown to be important for setting the initial biogeochemical conditions required for colonization by vegetation (Ciccazzo et al. 2016;Frey et al. 2013;Milner et al. 2007;Nicol et al. 2005). Furthermore, biofilms develop in glacially-fed streams and following retreat in Glacier Bay, possibly impacting stream biogeochemistry (Milner et al. 2007;Wilhelm et al. 2013). In light of declining glacial contributions to streamflow with increased stream age, autochthonous production could explain the marked rise in number of assigned formulae between GO 0 and ES 50-100 streams (Noriega-Ortega et al. 2019) and maintained high abundance of aliphatic and peptidelike compounds (combined mean of ES 50-100 streams 15.54 ± 1.95% RA; Table 1).
Declining d 13 C-DOC together with increased aromaticity and RA of HUPs compounds with increased stream age, suggested that terrestrial plants and soils became increasingly important OM sources with glacier retreat and subsequent terrestrial ecosystem succession. Terrestrial OM is visible in all ES 50-100 streams but dominates the molecular-level DOM composition in streams [ 100 years old (i.e., LS 100-210 streams; Table 1 ,Figs. 4,5). Plants colonize and soil organic carbon content is shown to increase within the first 50 years after glacier retreat (Bardgett et al. 2007;Milner et al. 2007). In Glacier Bay, * 150 years after glacier retreat soil organic carbon stocks and above-ground biomass are typically well established (Milner et al. 2007). Thus, terrestrial signatures dominated the molecular-level DOM composition of LS 100-210 streams and hence the RA of aliphatic formulae was far less than observed at previous stages in succession (Table 1,Figs. 4,5). It is likely that aliphatic formulae associated with either glacial-derived DOM, atmospheric deposition, or in situ production may have been undetected by FT-ICR MS, due to the increased prevalence of terrestrial aromatic compounds. Furthermore, aliphatic compounds may not be present or may have been removed from the streamwater DOM pool, since with glacier retreat channel conditions are likely to become more favorable for stream microbes (i.e., increased water temperatures and decreased turbidity), as shown by increased microbial diversity, abundance and biomass in non-glacial compared to glacial streams (Battin et al. 2004;Brighenti et al. 2019;Milner 1987;Milner et al. 2009Milner et al. , 2017. Thus, with glacier retreat stream microbial communities may be better able to metabolize biolabile DOM inputs. Together microbial metabolism and overprinting by terrestrial sources likely explained the stark trends of individual molecular formulae with stream age (Fig. 6), particularly that some aliphatic compounds were undetected in older streams.
By * 200 years since deglaciation, molecularlevel DOM composition was remarkably similar to that observed for OGF 1393 (despite over 1000 years difference in stream age), as well as temperate and arctic riverine DOM without permafrost influence (Behnke et al. 2020;Fellman et al. 2014;Spencer et al. 2008;Stubbins et al. 2012). Thus, following a similar timeseries of terrestrial ecosystem succession in Glacier Bay, stream DOM composition has shifted from glacier dominated to terrestrial vegetation and soil dominated. Interestingly, the successional changes in DOM composition highlight the majority of the shift has occurred after only 200 years (Figs. 2,3,4,5c). Thus, changes may still occur after 200 years, but our study suggests a consistent and rapid shift in DOM compositional metrics within the initial 200 years since deglaciation. Furthermore, it is unlikely that all watersheds across Glacier Bay will result in a single climax DOM composition, equivalent to OGF 1393 , as exemplified by the considerable variability in Fig. 5c. Successional differences are evident across Glacier Bay, where environments of a similar age have a variety of vegetation community compositions, landcover, and soil characteristics (Buma et al. 2017;Fastie 1995;Nagorski et al. 2014). It is expected that the quantity and quality of DOM will covary with these differences, along with the timeframes required to reach stasis (Fig. 5c). For example, the extent to which succession leads to upland forests or alpine tundra formation compared to peatlands or forested wetlands will likely determine the rate and magnitude of DOC increase and DOM compositional change (Aitkenhead et al. 1999;D'Amore et al. 2016;Fellman et al. 2008a;Hood and Scott 2008). Notwithstanding the natural variability in climax DOM composition and timings across Glacier Bay, late succession streamwater DOM composition will be dominated by terrestrial DOM signatures and thus will be distinct compared to glacier outflows and early succession streams (Fig. 5;Behnke et al. 2020;Fellman et al. 2010;Kellerman et al. 2021). It is likely that the much of this change in DOM composition will occur within the first two centuries since glacier retreat, regardless of the climax DOM composition, as shown by the large and consistent changes between the streams studied here (Fig. 5). Therefore, the evolution of DOM composition in these watersheds, places a molecular-level clock on DOM biogeochemistry in downstream ecosystems following watershed deglaciation.
Landscape controls on dissolved organic matter composition Heteroatom content declined with increased stream age (Table 1, Fig. 5), which is consistent with previous research showing that glacier-derived DOM contains large percent contributions from heteroatoms (Behnke et al. 2020;Bhatia et al. 2010;Kellerman et al. 2021;Spencer et al. 2014a;Stubbins et al. 2012). Despite this, CHON content was higher than the mean in LS1 138 , LS2 153 and LS3 168 streams, helping to explain the gradient in PC2 (Fig. 5). CHON content was not driven by stream age, but rather it was likely a function of other watershed and landscape characteristics. For example, vegetation with nitrogen-fixing microbes (e.g., alder species) colonize the landscape during early to late (50-150 years) succession (Bormann and Sidle 1990). This has been shown to result in an increase in soil and lake water total nitrogen concentrations in Glacier Bay (Engstrom et al. 2000;Milner et al. 2007) and thus, conceivably, a rise in abundance of streamwater CHON compounds. Additionally, salmon colonize streams in Glacier Bay (Nagorski et al. 2014) and when they occur in high densities (particularly in July-August) have been shown to seasonally influence DOM concentrations and composition by delivering protein-rich DOM to surface waters in the region (Behnke et al. 2020;Fellman et al. 2008bFellman et al. , 2014Hood et al. 2007). Thus, the enrichment in CHON observed in LS1 138 , LS2 153 and LS3 168 streams could be due to salmon-derived inputs of DOM that are rich in organic nitrogen (Hood et al. 2007). This salmon-derived DOM has also been associated with an enrichment in streamwater d 13 C-DOC and increased FI values (Behnke et al. 2020;Chaloner et al. 2002;Fellman et al. 2014;Hood et al. 2007). However, the salmon influence on DOM composition was not qualitatively observed in LS1 138 , LS2 153 and LS3 168 streams during the 3-day study period in early July, since there was no consistent shift in d 13 C-DOC or FI values between the streams (Table 1). Nitrogen-fixing vegetation are eventually outcompeted by later succession conifers and shrubs, such as spruce, hemlock and blueberry species, resulting in a decline in soil and lake total nitrogen concentrations (Engstrom and Fritz 2006;Milner et al. 2007). Although vegetation species with nitrogen-fixing microbes (e.g., alder) have been strongly linked with stream total nitrogen concentrations (e.g., Compton et al. 2003;Shaftel et al. 2012), more research is required to explore linkages with CHON content of DOM. Nonetheless, nitrogen-fixing vegetation provides a plausible explanation for the transient increase in CHON content in LS 100-210 streams (Table 1, Fig. 5). The timing of this replacement may be linked to landscape variables, such as distance from the competing species' propagules, rather than stream age (Fastie 1995).

Implications for downstream carbon cycling with continued glacier loss
Marked changes to DOM composition have been shown within the first * 200 years since glacier retreat, where beyond this point in landscape development DOM composition underwent minimal change. The timeframe of the evolution in DOM composition was consistent with terrestrial ecosystem succession in Glacier Bay, where within the first 50 years there is increased biotic control on the landscape, particularly due to gradual re-vegetation and the formation of soil OM (Chapin et al. 1994;Milner et al. 2007). The observed trends in DOM composition across the streams studied here, suggest that there are reasonably predictable shifts to DOM composition as glaciers retreat and terrestrial ecosystem system succession ensues. We expect these trends are broadly applicable to other deglaciating catchments where low elevation forest communities typically establish following glacier retreat, but especially the watersheds that ring the Gulf of Alaska.
Southeast Alaska is experiencing rapid glacier loss, with hydrological regimes shifts predicted to occur this century (i.e., shift from glacier dominated towards typical snow, rain or groundwater fed hydrograph; Bliss et al. 2014;Edwards et al. 2013;Zemp et al. 2019). Our findings suggest that there will be predictable changes to stream DOM composition on this centennial timeframe. Notably, there is likely to be a shift from ancient, highly bioavailable aliphatic dominated DOM exported to downstream ecosystems towards terrestrial-derived, modern, relatively aromatic enriched DOM. Ancient glacier-derived DOM has been shown to be assimilated into proglacial foodwebs (Fellman et al. 2015b;Hågvar and Ohlson 2013). However, the extent to which this seasonal supply of glacier DOM currently subsidizes the downstream microbial food web is unclear. Furthermore, it is uncertain whether declining DOM bioavailability, associated with increased terrestrial DOM contributions concurrent with glacier retreat, will be in part compensated by an enhanced annual DOM supply downstream (Edwards et al. 2021). The shift towards aromatic-rich DOM derived from modern plant and soil material has unknown consequences for coastal carbon dynamics, heterotrophic productivity and thus the highly productive marine foodwebs of the Gulf of Alaska.
Funding This study was supported by NSF DEB (Grant Nos.