Storm-driven hydrological, seasonal, and land use/land cover impact on dissolved organic matter dynamics in a mid-Atlantic, USA coastal plain river system characterized by 21 T FT-ICR mass spectrometry

phytoplankton-derived compounds). Collectively, these results demonstrate how arti ﬁ cial reservoirs alter the characteristics of DOM exported from rivers with implications for understanding carbon export and fate at river-estuary interfaces.


Introduction
Dissolved organic matter (DOM) is the dominant energy source (Fisher et al., 1973), a transporter of pollutants (Sigleo and Means, 1990), and the largest organic carbon reservoir (Hansell and Carlson, 2015) in the aquatic environment.Large amounts of riverine DOM (approximately 2.5 × 10 14 g-DOC year −1 ; Hedges et al., 1997) are transported to estuaries and the ocean surface with impacts for global ocean productivity, carbon budgets, and carbon cycles.Riverine DOM originates from a range of allochthonous and autochthonous sources including terrestrial plants and soils (Hedges, 1992;Hedges et al., 1997;Barnes et al., 2018), primary and secondary production (Peterson et al., 1994;Raymond and Bauer, 2001;McCallister et al., 2004), groundwater (Brunke and Gonser, 2003;Inamdar et al., 2011;Lee et al., 2019), and atmospheric deposition (Wozniak et al., 2012;Mitra et al., 2017;Czarnecki et al., 2023).The molecular characteristics and quantities of DOM from these sources vary in space and time.DOM retrieved from plant litter and upper organic-rich soil layers have been found to be enriched in modern aromatic and lignin-like materials (Barnes et al., 2018;Qin et al., 2020).Bio-sourced DOM, resulting from primary production and heterotrophic degradation, contains higher abundances of protein-like, carbohydrate-like, and microbial humic-like compounds with relatively low molecular weights (Peterson et al., 1994;McCallister et al., 2004;Minor and Oyler, 2023).Consequently, DOM composition along a river-to-estuary transit depends on the relative contributions from autochthonous and allochthonous sources to rivers and their watersheds.Importantly, DOM composition influences its cycling and fate.Ahamed et al. (2023), for example, found that DOM bioavailability was influenced by the C/N ratio, C number, and thermodynamic favorability of a molecule.Further, the compound aromaticity has been associated with both higher photodegradability (Stubbins et al., 2010;Goranov et al., 2020) and resistance to microbial degradation (Bostick et al., 2021).Understanding controls on DOM composition exported from rivers is therefore important for the evaluation of carbon cycling in coastal environments especially in the context of continued human modification of the environment and Earth's changing climate.
Hydrological events (i.e., storm, rainfall, snowmelt, etc.) have emerged as a key factor impacting the quantity, characteristics, and cycling of riverine DOM via changing the dominant water flow paths (Barnes et al., 2018;Hosen et al., 2021a, b), enhancing the leaching of organic-rich soil layers and surface vegetation (Hernes et al., 2008;Singh et al., 2014), and shortening DOM residence time (Boyer et al., 2001;Hood et al., 2006).Raymond and Saiers (2010) estimated that 86% of annual DOC export occurs during hydrological events in New England (USA) forested watersheds, spotlighting the importance of event-driven effects on DOM export and cycles.The significant relationship between DOM export and storm-driven high discharge events results in a well-known theoretical framework, the pulse-shunt concept (Raymond et al., 2016).The concept highlights that high discharge triggered by storm events can (1) prompt a large flux of modern terrestrial DOM from surficial sources (i.e., throughfall, litter leachate, and soil water) into drainage networks (known as the pulse process), (2) decrease water residence time allowing the DOM to escape from microbial and photochemical transformations and reach estuaries and coastal oceans in a less degraded form (known as the shunt process), and (3) relocate the biogeochemical hotspots for DOM processing from the upper to the lower reaches of the river systems (Raymond et al., 2016;Wagner et al., 2019).The pulseshunt process has been found to be a primary control on DOM concentration and composition regardless of spatial variations in temperate upland landscape units (i.e., Raymond and Saiers, 2010;Raymond et al., 2016;Hosen et al., 2021a).Wagner et al. (2019) further hypothesized that the pulse-shunt concept alters patterns of molecular diversity found along the river continuum (the River Continuum Concept; Vannote et al., 1980) yielding more homogenized characteristics across river orders.
However, recent work suggests the pulse-shunt concept is not always the dominant driver of riverine DOM dynamics in coastal plain watersheds.For example, increased river discharge during storm events had no significant relationship with DOM quantities in the Delaware Estuary (Osterholz et al., 2016) and the Waccamaw River system (Majidzadeh et al., 2017).The hydrological effects on riverine DOM variations were strongly impacted by changes in seasons and natural LULC units (forested/agricultural lands and wetlands) (Majidzadeh et al., 2017;Hosen et al., 2018;Fasching et al., 2019).In rivers with primarily agricultural and urban land uses, Fasching et al. (2019) modeled that the pulse process becomes invalid if watersheds contain more than 30%-40% non-cropland (i.e., wetlands, forested, and grasslands) within the catchments.Spatial variability among LULC units and associated soil biogeochemistry variations have been further shown to impact DOM sources and mobility to rivers, subsequently affecting DOM composition (Stahl et al., 2021;Shousha et al., 2022;Ma et al., 2023).The magnitudes and characteristics of allochthonous and autochthonous DOM sources to rivers also vary temporally due to influences on light, temperature, and associated biological (and human) activity (Power et al., 2018;Bhattacharya and Osburn, 2020).The relative importance of storm-driven discharge, time of year, and LULC units on DOM dynamics may thus fluctuate depending on location, anthropogenic effects, and particular circumstances.With predicted increasing urbanization and farmland expansion in coastal plain landscapes (Griffith et al., 2003;Fisher et al., 2007), incorporating human impacted LULC units (developed/agricultural lands) and associated drainage networks as functional units along with forested lands, seasons and hydrological factors is essential to better model DOM variations in coastal plain riverine systems.
Artificial reservoirs (i.e., millponds, artificial lakes, stormwater retention water bodies), another class of major anthropogenic interventions influencing riverine processes, may also disrupt the pulse-shunt process.Dams in millponds and artificial lakes coupled with frequent irrigation can decrease water flow rates and retain stormwaters, leading to extended water residence times (Friedl and Wüest, 2002;Andres et al., 2019).While storm events can pulse terrestrial materials (i.e., DOM and nutrients) and slightly shorten water residence times in these artificial reservoirs, the storm-derived shunt effects may be muted relative to those observed in free-flowing streams.The prolonged water residence times in millponds and artificial lakes create a lentic habitat that fragments the river continuity (Friedl and Wüest, 2002) and allows for very high primary productivity with associated microbial activity (Mansour et al., 2018), particularly during warmer months [i.e., spring and summer; (Oliver et al., 2016;Andres et al., 2019)].These artificial reservoir systems, therefore, likely (1) contain more algal-influenced DOM than those of free-flowing river systems, (2) export large quantities of algal OM after storms that has accumulated in the interim, and (3) receive ample inputs of nutrients with storms that promote even higher primary production rates.As an overall result, in river systems influenced by millponds and artificial lakes, the seasonal pattern of autochthonous production becomes more important in determining the DOM characteristics, with storm events likely remaining as the dominant controller on DOM quantities.There are at least 2.6 million artificial waterbodies across the United States landscape (Smith et al., 2002), with rapid construction predicted to increase those numbers (Oswald, 1995).Considering damming effects of artificial reservoirs within the pulse-shunt concept can provide important implications for explaining DOM dynamics and regulating its sources in river systems across the United States and beyond.
This study investigates the integrated roles of LULCs, storm/ precipitation-driven changes in river discharge, and time of year on the concentration and molecular composition of DOM in the Murderkill River (MR), a typical mid-Atlantic coastal plain river system.Artificial lakes and remnants of historical millponds are also important components of the MR system that act as stormwater retention water bodies leading to very long water residence times (Andres et al., 2019).Absorption and fluorescence spectroscopy of DOM were used along with 21 T Fourier-transform ion cyclotron resonance mass spectrometry (FT-ICR MS) to evaluate DOM composition, sources, and dynamics in the MR.The 21 T FT-ICR MS is the highest resolving mass analyzer for DOM characterization, and routinely achieves resolving power that exceeds 2,000,000 at m/z 400, which enables identification of tens of thousands of DOM species at sub-ppm mass measurement accuracy (Hendrickson et al., 2015;Smith et al., 2018;Bahureksa et al., 2022).Our first working hypothesis was that storm/ precipitation events drive variations in DOM quantities and composition with increased river discharge after storms, yielding higher quantities of DOM characterized as more aromatic, oxygenated, and humic DOM relative to lower discharge periods.We also hypothesized that seasonal patterns of autochthonous production would be a very important factor changing the composition of DOM pools in the MR system due to the damming effect of artificial reservoirs, with a higher contribution of protein-like, aliphatic or highly unsaturated DOM compounds during the more biologically active spring through late summer (Stenson et al., 2003;Koch et al., 2007;Nebbioso and Piccolo, 2013).

Study sites
Sampling sites are located at pond spillways and streams within the MR system, Delaware (DE), USA (Figure 1).The MR watershed is a mid-Atlantic coastal plain environment that covers approximately 280 km 2 within the greater DE Bay watershed and is predominantly composed of agricultural (52.4% of land cover), developed (15.5% of land cover), and forested (10.8% of land cover) lands (Rogerson et al., 2011).Soil biogeochemistry and riverine inorganic nutrient composition were found significantly different across different LULC units in the system (Fischer, 2014) and could be important factors observed in this study in context of LULC effects.Sampling sites were selected to sample waterways draining sub-watersheds disproportionately impacted by these agricultural (n = 3), developed (n = 1), and forested (n = 2) LULC units according to the most recent Department of Natural Resources and Environmental Control (DNREC) MR watershed report (Rogerson et al., 2011; Figure 1).The sites are selected from within the mainstem of the MR, which contributes approximately 64.2% of the total freshwater inflow into the downstream MR estuary, and from Spring Creek contributing an additional 23.4% of the discharge (Dewitt and Daiber, 1974).The freshwater discharge in the MR system was found to strongly respond to storm events within a day (Voynova et al., 2015), thus yielding ideal conditions to study the storm-driven hydrological controls on DOM compositions.Sampling sites also capture various waterbodies within the MR system including (1) free-flowing streams (Spring Branch, Spring Creek, and Market Street), (2) small millponds (Coursey Pond and McCauley's Pond;Andres et al., 2019), and (3) artificial lakes (Andrews Lake).The water flow rate in both millpond and artificial lake sites (fluctuating around 1.0 m 3 /s; Andres et al., 2019) are significantly slower than that of freeflowing streams (average = 8.1 m 3 /s; Dewitt and Daiber, 1974), particularly during summer at millponds due to increased irrigation (<0.5 m 3 /s; Andres et al., 2019).Notably, the millponds and lake sampling sites are upstream of the free-flowing streams, and their influence on DOM composition may be observed at the stream sites.

Sampling, LULC, and hydrology
Surface water samples (~500 mL) were collected seasonally and before/after storm events (n = 13 for Andrew's Lake, n = 14 for Coursey Pond, McCauley's Pond & Spring Creek, n = 15 for Market Street & Spring Branch) using acid-cleaned (1 M HCl) high density polyethylene (HDPE) bottles (Thermo Fischer Scientific Nalgene ® ) at 6 selected sites draining representative LULC in the MR system from 21 November 2019 to 23 August 2021 (Figure 1).Seasons were defined as-spring (March 20 to June 19), summer (June 20 to September 22), fall (September 23 to December 20), and winter (December 21 to March 19) according to the astronomical calendar in the Northern hemisphere.Hourly tidally filtered river discharge (m 3 s −1 ), the discharge rate calculated after removing short-term variations in instantaneous discharge due to tidal effects from the Delaware Bay (Ruhl and Simpson, 2005;Dzwonkowski et al., 2014), for the MR was obtained from the United States Geological Survey database at Frederica, Kent County, DE (waterdata.usgs.gov).The mean tidally filtered river discharge during the collection time was calculated and used to represent the river discharge on collection dates.The sampling events right after/during storm events under high discharge (relative to non-storm conditions in the same season) were classified as post-storm conditions (2.4-17.5 m 3 s −1 ).Sampling events showing lower river discharge (relative to post-storm in the same season) occurred between two storm events representing nonstorm conditions (0.3-6.1 m 3 s −1 ).Samples were vacuum filtered within 24 h via Whatman GF/F 0.7 µm glass fiber filters that were pre-combusted at 450 °C.Approximately 50 mL filtrate was refrigerated (at ~4 °C) until quantification of DOC and spectroscopic analyses (<1 week).The remainder was stored in aluminum-foil-covered, acid-washed HDPE bottles (Thermo Fischer Scientific Nalgene ® ) and frozen at −20 °C for solid phase extraction and molecular formula level characterization via FT-ICR MS.

Dissolved organic matter (DOM)
Filtrate (~20 mL) for DOC quantification were acidified to pH 2 using 2 M HCl.DOC concentrations ([DOC]) were measured by high temperature catalytic oxidation (Emery et al., 1971;Sharp et al., 1993) using a Shimadzu TOC-Vcph analyzer.Prior to analysis, the acidified samples were purged for 1 min with purified air to remove dissolved inorganic carbon.[DOC] was quantified as non-purgeable organic carbon by non-dispersive infrared detection.Instrument responses were calibrated for [DOC] using a potassium hydrogen phthalate solution.Aliquots of surface seawater (SSR; Batch 21, Lot# 04-21) and mid-depth seawater reference (MSR; Batch 18, Lot# 04-21) samples from the Consensus Reference Material Project of Hansell Organic Biogeochemistry Laboratory (University of Miami, Florida, USA) were analyzed during each run to monitor instrumental precision and accuracy.Analyses of the SSR and MSR deviated by less than 5% from the reported values for these standards (73-75 μM and 61-63 µM DOC, respectively).
Absorption and excitation-emission matrix spectra (EEMs) were simultaneously obtained using a spectrofluorometer (Aqualog, Horiba) to optically characterize the chromophoric and fluorophoric fractions of DOM (CDOM and FDOM, respectively).Low DOC (<2 µM) water (Millipore) was used as a blank to correct spectral intensities for CDOM and FDOM analyses.CDOM absorption spectra from 230 to 700 nm in 2 nm intervals were measured and auto-corrected to real signals by subtracting the dark-current signals via the Aqualog software.Specific ultraviolet absorbance (SUVA 254 ) was evaluated by normalizing the absorbance (a m −1 ) at 254 nm using [DOC] (mg-C L −1 ) to estimate DOM aromaticity following previous work (Weishaar et al., 2003;Helms et al., 2008).The ratio of absorbance at 250-365 nm (called E 2 :E 3 ) was computed and used as a proxy for molecular weight and sizes (De Haan and De Boer, 1987;Helms et al., 2008).EEMs were measured with a 150-W Xenon lamp at 2 nm excitation wavelength intervals between 230 and 700 nm, with emission data collected at 4.65-nm intervals from 245 to 822 nm.EEMs were normalized to Raman Units, and corrected for inner filter effects, 1st and 2nd order Rayleigh scattering using processing tools in Aqualog software.The corrected EEMs were interpolated (Lee et al., 1997;Elcoroaristizabal et al., 2015) and smoothed using the "staRdom" (Pucher et al., 2019) and "eemR" packages.Coble peaks (Coble et al., 2014), humification index (HIX; Ohno, 2002), and biological index (BIX; Huguet et al., 2009) of processed EEMs were analyzed using the "staRdom" package.HIX is a proxy for DOM humic contents and was computed as the fluorescence intensity in emission 300-345 nm region divided by the sum of intensity in emission 300-450 nm and 435-480 nm regions (Ohno, 2002).BIX as an indicator of recent autochthonous DOM contribution was calculated as the ratio of fluorescence at excitation 310 nm and emission at 380 nm to that at excitation 310 nm and at emission 430 nm (Huguet et al., 2009).
Molecular level characterization of DOM was performed with a custom-built hybrid linear ion trap FT-ICR mass spectrometer equipped with a 21 T superconducting solenoid magnet (Hendrickson et al., 2015;Smith et al., 2018) at the National High Magnetic Field Laboratory (NHMFL), Florida, USA.A pair of non-storm and post-storm sample filtrates at selected sampling sites (n = 5, excluding Spring Creek) in spring, summer, fall, and winter (n = 40; Supplementary Figure S1) were thawed under 4 °C in the refrigerator.Thawed sample filtrates were acidified to pH 2 using 1 M HCl and solid phase extracted following Dittmar et al. (2008) using a 100 mg pre-packed PPL column (Agilent Technologies).The resulting methanol extracts (or PPL-isolated DOM samples) were diluted to a [DOC] of ~50 mg-C L −1 in methanol (HPLC-grade, Fischer Scientific, USA), stored in acid-washed, pre-combusted (at 450 °C), aluminum-foil-covered glass vials, and shipped on ice to the NHMFL for FT-ICR MS analysis.Sample extracts were infused at 500 nL min −1 into a micro-electrospray source (Emmett et al., 1998) with the external multipole ion guide for 1-5 ms and released m/ z-dependently (via decrease of an auxiliary radio frequency potential the multipole rods and the end-cap electrodes) to the harmonized ICR cell with 6 V trapping potential (Boldin and Nikolaev, 2011;Kaiser et al., 2013).Data was conditionally coadded in broadband mode (m/z 150-1800) using the Predator data station (Blakney et al., 2011) with 100 time-domain acquisitions averaged for all experiments, with a detection period.Mass spectra were phasecorrected (Xian et al., 2010) and internally calibrated with 10-15 highly abundant homologous series that span the molecular weight distribution based on the "walking" calibration method (Savory et al., 2011).Experimentally measured masses were converted from the International Union of Pure and Applied Chemistry mass scale to the Kendrick mass scale (Kendrick, 1963) for rapid identification of homologous series for each heteroatom (i.e., species with the same C c H h N n O o S s content and degree of aromaticity, differing only by degree of alkylation; Hughey et al., 2001).The heteroatom classes, types (referring double bond equivalents, DBE, calculated from neutral elemental compositions; McLafferty et al., 1994), and carbon numbers were tabulated for subsequent generation of relative heteroatom class relative abundance distributions, graphical relative-abundance weighted images, and van Krevelen diagrams (van Krevelen, 1950;Kim et al., 2003).Peaks with signal magnitude greater than 6 times the baseline root-mean-square (rms) noise at m/z 500 were exported to peak lists.Molecular formula assignments and data visualization were performed with PetroOrg software (version 18.0.5;Corilo, 2014), and assignments with an error <0.5 parts-per-million (ppm) were considered in this study.For all samples, the achieved resolving power is at m/z 400 > 2,000,000.All spectra presented result in 16,272-28,810 assigned elemental compositions with 0.034-0.046ppm RMS error.The sample name, number of assigned species, achieved resolving power, and mass measurement accuracy for all samples are listed in Supplementary Table S3.All FT-ICR MS raw files, calibrated peak lists, and assigned elemental compositions are publicly available via the Open Science Framework (osf.io/kwpr6) through DOI: 10.17605/OSF.IO/KWPR6.

Statistical analyses 2.4.1 Linear regression analysis
Linear regressions were performed in R (version 4.2.2) using "lm" (Wilkinson and Rogers, 1973;Chambers, 1992) to illustrate the relationships between river discharge and [DOC] as well as optical parameters evaluated based on EEMs measurements.Statistical hypothesis tests were conducted simultaneously with linear regression model fitting in R (version 4.2.2).

Parallel factor analysis (PARAFAC)
A PARAFAC model was developed for all processed EEMs (n = 97 including replicates) in R (version 4.2.2) using packages "staRdom" and "eemR" (Pucher et al., 2019).The number of components was determined by statistical check of correlation between each component and visual examination of EEMs and residuals from models with increasing or decreasing numbers of components (Stedmon and Bro, 2008).The model with determined components was validated via random split-half analysis (Murphy et al., 2013).Intensities of each component were normalized to the maximum fluorescent intensity of each component (Stedmon and Bro, 2008;Pucher et al., 2019).The component EEMs were tested against the OpenFluor database (www.openfluor.org;Murphy et al., 2014).

FT-ICR MS data analysis
The modified aromaticity index (AI mod ; Koch and Dittmar, 2006), DBE (McLafferty et al., 1994;Bae et al., 2011), O/C, and H/C ratios were calculated to approximate the structure of each DOM formula.The study assigned sample compounds into 4 main aromaticity-based groups [adapted from (Koch and Dittmar, 2006;Osterholz et al., 2016)] with each being further divided into 3 O/C ratio-based subcategories including (1) Due to the selectivity in the electrospray ionization process and sample processing, peak signal intensity/abundance of individual samples is not equivalent to its concentrations.The peak relative abundance, however, has been found useful for comparing DOM compositions measured under the same electrospray conditions, instrumental parameters, and solvents (Fievre et al., 1997;Sleighter and Hatcher, 2007;D'Andrilli et al., 2010;Pan et al., 2020), and thus it is used in this study to evaluate the percentage of individual categorized compound classes to the total assigned peak abundance.Tukey's honestly significant difference test (Tukey's HSD) coupled with ANOVA was performed on R (version 4.2.2) using "TukeyHSD" (Miller, 1981;Yandell, 1997) in this study to examine pairwise differences in H/C, O/C, and compound classes among river discharge, season, and LULC conditions.The percentage data were arcsine transformed (or the angular transformed) prior to use in Tukey's HSD test.

Principal component analysis (PCA)
PCA was performed in R (version 4.2.2) using "prcomp" (Mardia et al., 1979;Venables and Ripley, 2002) and visualized using package "factoextra".PARAFAC loading scores of each identified component from all samples were input as variables in PCA.The Pearson's product-moment correlations (r) between PARAFAC component loadings and PCA component scores were tested using "cor.test"function on R (Best and Roberts, 1975;Hollander et al., 2015).

DOC and optical measurements
[DOC] ranged from 45.2 to 1,525.8 μM and positively correlates with river discharge (Figure 2A).SUVA 254 ranged from 2.7 to 5.8 L mg-C −1 m −1 , and E 2 :E 3 ratio values ranged from 4.0 to 6.4.HIX and BIX indexes evaluated from EEMs data ranged from 5.3 to 24.4 and 0.49 to 0.75, respectively.A complete list of [DOC] and described optical parameters for all samples is detailed in Supplementary Table S1.The relationship between optical parameters and river discharge during collection times is displayed in Figures 2B-E with SUVA 254 and HIX showing statistically significant positive correlations, and E 2 :E 3 and BIX showing statistically significant negative correlations (linear regression, p < 0.001 in all cases).
Data inputs in PCA were further categorized according to discharge, season, and LULC conditions to determine the major drivers of EEM-PARAFAC component variances (Figure 4).PC1, as indicative of FDOM abundances, explains the variance between samples collected under non-storm and post-storm conditions (Figure 4A), with post-storm EEM-PARAFAC sample scores trending toward the positive axis of PC1 and those of nonstorm toward its negative axis.PC2 explains the seasonal variances of EEM-PARAFAC loadings (Figure 4B).Fall and winter EEM-PARAFAC sample PC2 scores are primarily positive and driven by C4, C3, and C1 vectors, while those of summer samples trend negative and align with the C5 and C6 vectors (Figure 4B).The LULC unit mean sample scores overlap (Figure 4C) implying a minor role for LULC in explaining FDOM compositional variations.

Ultrahigh resolution mass spectral data
Between 16272 and 28810 molecular formulae (Supplementary Figure S4A) were assigned to FT-ICR MS peaks identified in each PPL-isolated DOM sample spectrum, with a mean error per Spectral excitation and emission properties of six fluorescent components identified by PARAFAC analysis.The excitation and emission intensity loadings of each PARAFAC components were shown in Supplementary Figure S2., D).A van Krevelen diagram of all assigned molecular formulae is shown in Figure 5, coupled with histogram distribution of H/C and O/C ratios grouped by discharge conditions.The distributions of O/C ratios shift toward the higher values from non-storm to poststorm PPL-isolated DOM samples, with mean post-storm O/C significantly higher than that of non-storm (ANOVA-Tukey's p < 0.001; Figure 5).Non-storm PPL-isolated DOM sample H/C ratios skew higher than post-storm samples, and the mean poststorm H/C value is significantly lower than that of non-storm (ANOVA-Tukey's p < 0.001; Figure 5).The H/C and O/C ratios vary across different season and LULC conditions (Supplementary Figures S5, S6).Under post-storm conditions, PPL-isolated DOM H/C distributes toward lower values in the order: spring > summer > winter > fall (Supplementary Figure S6A) and developed > agricultural > forested (Supplementary Figure S5A; ANOVA-Tukey's p < 0.001).Mean post-storm PPL-isolated DOM O/C ratios decrease in the order fall > winter > summer > spring (Supplementary Figure S6B; ANOVA-Tukey's p < 0.01) and forested > agricultural > developed (Supplementary Figure S5B; ANOVA-Tukey's p < 0.001).For non-storm conditions, the mean H/C ratios of PPL-isolated DOM collected at forested sites is significantly higher than those collected at agricultural and developed sites (Supplementary Figure S5A; ANOVA-Tukey's p < 0.001), with no significant difference of mean O/C values found between different LULC conditions (Supplementary Figure S5B).The mean H/C of non-storm PPL-isolated DOM samples increases in the order: fall < spring < summer < winter (Supplementary Figure S6A; ANOVA-Tukey's p < 0.001), and their mean O/C decreases in the order: winter > spring > summer > fall (Supplementary Figure S6B: ANOVA-Tukey's p < 0.01).
All assigned molecular formulae of the individual PPL-isolated DOM sample contained one or more O, N, and S atoms, with 54.7-79.6CHO%, 10.0 to 21.3 CHON%, 4.0 to 23.8 CHOS%, and 0.5 to 4.1 CHONS% of total relative abundances of assigned molecular formulae across all PPL-isolated DOM samples (Supplementary Table S3; Supplementary Figure S7).There is no significant difference of CHO%, CHON%, CHOS%, and CHONS% found between non-storm and post-storm conditions (ANOVA-Tukey's p > 0.05; Supplementary Table S8).CHO% significantly changes seasonally coupled with LULC effects (ANOVA-Tukey's p < 0.05), with highest mean CHO% in fall PPL-isolated DOM samples collected at forested sites (77.6%), lowest mean CHO% in spring PPL-isolated DOM samples at developed sites (59.9%;Principal component analysis based on a correlation matrix for all samples using PARAFAC loadings (Supplementary Figure S3) as input data.Shown are principal component 1 (PC1) and principal component 2 (PC2), accounting for 83.6% and 13.0% of the total variance, respectively.Sample PARAFAC scores are grouped by (A) discharge conditions, (B) season, and (C) LULC.Ellipses represent mean scores of classified groups.The arrows represent loading vectors for each of the components with length of the vectors corresponding to the loading intensity.6B).The highest mean sum of CHON%, CHOS%, and CHONS% presented in spring PPL-isolated DOM samples under agricultural conditions (38.2%; Figure 6B), with lowest mean values in fall forested PPL-isolated DOM samples (21.6%; Figure 6B).
Differences in arcsine transformed percentage abundances of compound categories across different discharge, season, and LULC  S9, S10.
conditions were examined via Tukey's HSD test (Figures 7, 8; Supplementary Tables S11-S15).The percentage abundances of most compound categories showed no significant difference across post-storm and non-storm conditions, except HO CondAromatic formulae groups (Supplementary Table S13).The percentage abundances of HO CondAromatic formulae in post-storm PPLisolated DOM samples is significantly higher than those of nonstorm conditions (Supplementary Table S13; ANOVA-Tukey's p =  S11.0.047).No significant difference in percentage abundances of all compound categories were found across different LULC units (Supplementary Table S15; ANOVA-Tukey's p > 0.05), however, significantly more HO aliphatic compounds were observed in spring and winter samples compared to fall (Supplementary Table S14; ANOVA-Tukey's p < 0.05).The combined influences of season and discharge conditions significantly affect the contributions of LO CondAromatic, MO Aromatic, MO CondAromatic, HO Aliphatic, and HO Aromatic compounds in all PPL-isolated DOM samples (Figure 7; Supplementary Table S11; ANOVA-Tukey's p < 0.05).Significantly lower percentage abundances of LO CondAromatic (mean = 0.23%), MO Aromatic (mean = 8.72%), MO CondAromatic (mean = 2.50%), and HO Aromatic (mean = 4.36%) formulae were found in winter non-storm PPL-isolated DOM samples relative to those of other season-discharge pairings (Figure 7; Supplementary Table S11; ANOVA-Tukey's p < 0.05).The lowest percentage abundance of HO Aliphatic formulae is in fall post-storm PPL-isolated DOM samples (mean = 0.63%), with highest values under winter non-storm conditions (mean = 1.12%; Figure 7).Combining discharge, season, and LULC conditions for comparison, the highest HO Olefinic % was found in winter forested PPL-isolated DOM samples under post-storm conditions (mean = 24.5%),with lowest values (mean = 9.7%) in winter developed PPL-isolated DOM samples under non-storm conditions (Figure 8; Supplementary Table S12; ANOVA-Tukey's p < 0.05).

Storm events and seasons primarily control [DOC] and FDOM compositions
The optical data described herein provide evidence for both storm and seasonal influences on [DOC] and FDOM composition.Rainstorm events can trigger rising river discharge in the mid-Atlantic US coastal plains (Sheridan, 1997;Supplementary Figure S1) and modify [DOC] and DOM composition exported to rivers, including the MR system.Positive correlations between river discharge and [DOC] (Figure 2A) as well as the contributions of high molecular weight (as measured by E 2 : E 3 ), aromatic (as measured by SUVA 254 ) CDOM (Figures 2B, C) and humic (as measured by HIX) FDOM (Figure 2E) are consistent with the storm influence (Inamdar et al., 2011;Raymond et al., 2016;Hosen et al., 2021a).Similarly, the inverse correlation between BIX (relative contributions from recent autochthonous production) and rising river discharge (Figure 2D) results presumably due to storm-mobilized contributions from lands (Dixon et al., 2014;Raymond et al., 2016;Pang et al., 2021).However, the PARAFAC PCA post-storm samples distribute among both positive and negative PC2 values with a mean value near 0, suggesting no dominant influence from the plant-litter derived humic-like C3 and C4 which have very positive PC2 loadings (Figure 4A).This lack of trend between post-storm samples and plant/ soil-derived humic-like FDOM is in apparent conflict with the positive  S12.
correlation between storm-induced high discharge and HIX, highlighting potential ambiguities in FDOM indices and the need to evaluate EEMs datasets thoroughly using multiple measures.In this case, the PCA results may be skewed toward summer post-storm events which show high C6 contributions revealing a temporal influence on FDOM regardless of storm status (Figures 4A, B).Here, the combined seasonal-storm influence likely results from the long residence times in the MR system artificial reservoirs and the storm-induced delivery of inorganic nutrients needed for primary production (Voynova et al., 2015).The higher C6 contributions during summer likely represent an algal response to the increased nutrient loadings (Chen et al., 2018) to the study streams and ponds (Andres et al., 2019;Musser, 2020).Notably, the summer post-storm discharge was the lowest across all seasons (Supplementary Figure S1) potentially due to growing season agricultural irrigation practices that further reduce discharge from millponds (Andres et al., 2019), suggesting that storm pulses may not have been strong enough to shunt the nutrients and DOM downstream.Below, we elaborate on these storm and season influences on FDOM.
Previous work suggests storm events to increase DOC quantities via changing the dominant river input flow paths to shallower soil layers (through micropores; Friesen, 2020) and overland flows (Barnes et al., 2018).With increasing depths, soil OM profiles commonly become highly decomposed reflected by low C densities, increasing δ 13 C, decreasing C: N values (Baldock et al., 1997;Baisden et al., 2002;Wynn et al., 2005) and lesser proportions of lignin (Baldock et al., 1992;Kiem and Kogel-Knabner, 2003;Klotzbücher et al., 2016).The storm-induced shallow soil and overland flow paths thus enhance the leaching of the organicrich layers and surface vegetation (Kaiser and Guggenberger, 2005;Hernes et al., 2008), resulting in the observed higher inputs of aromatic, humic, and high-molecular-weight DOM.By contrast then, deeper water flow paths during non-storm conditions leach soils with lesser organic content from mineral rich, deeper portions of the subsurface (Jeanneau et al., 2015;Barnes et al., 2018) and yield DOM that is aged, highly processed soil-OM, and enriched in compounds with low molecular weights, aromaticity, and C:N ratios (Neff et al., 2006;Sanderman et al., 2008Sanderman et al., , 2009)).The loss of lignin with increasing degree of decomposition (and depth) may explain the decreasing aromaticity (SUVA 254 ) and molecular weights (E 2 : E 3 ) in this study during non-storm conditions when deeper river flow paths prevail.The overland and shallow subsurface flows during storm events are more likely to encounter less degraded lignin-and tannin-derived aromatic components due to enhanced leaching of shallow soil layers and surface vegetation (Hernes et al., 2008).
The positive relationship between storm associated river discharge and allochthonous DOM proxies has been observed in freshwater upland systems (i.e., Raymond and Saiers, 2010;Raymond et al., 2016) and estuarine systems in mid-Atlantic USA coastal plains (i.e., Osterholz et al., 2016;Hosen et al., 2018).Here, we find a similar relationship within freshwaters from a small mid-Atlantic coastal plain watershed providing evidence for the ubiquity of the relationship.However, the linear correlations of river discharge with DOM optical parameters in our study show considerable variances (r 2 ≤ 0.44), with particularly low r 2 values in river discharge's linear fitting with aromaticity (SUVA 254 ; r 2 = 0.23), molecular weight (E 2 : E 3 ; r 2 = 0.21), and humic contents (HIX; r 2 = 0.22).The low r 2 values are indicative that river discharge is not the sole factor shaping CDOM and FDOM characteristics in the system.In the MR system, the river discharge obtained from free-flowing streams may not accurately capture the flow rate and retention time within upstream artificial reservoir sites (i.e., Coursey Pond, McCauley's Pond, and Andrews Lake).The relatively long residence time of artificial lake (Stedmon et al., 2006;Kellerman et al., 2014;Kothawala et al., 2014) and millpond (Ortega et al., 2018;Andres et al., 2019) waters can enhance autochthonous production during warmer months (i.e., summer and spring) and subsequently amplify the seasonal patterns of DOM in downstream rivers.Therefore, other factors, such as the season of collection and associated biogeochemical transformations, must play a role along with discharge conditions in determining DOM quantities and composition to elucidate the considerable variances in correlation fitting, particularly for a river system.
While the PCA shows storm events to be an important driver determining total FDOM abundance (PC1, Figure 4A), season of collection was more important for distinguishing the relative abundance of different PARAFAC-FDOM components (PC2, Figure 4B).The humic-like winter and fall FDOM characteristics likely result from elevated inputs of the microbial breakdown products of terrestrial higher plants and soils.This observation is evidenced by the prevalence of the soil derived and lignin-related humic-like components (C1 and C3) and the occasional appearance of the microbial humic-like FDOM (C4) known to be produced during microbial degradation of allochthonous OM in rivers (Yamashita et al., 2013;Vines and Terry, 2020;Smith et al., 2021).In summer, the increasing amounts of tryptophan/tyrosine-like C6 component likely derive from the rapid growth of phytoplankton responding to increased light and temperature (Seitzinger et al., 2002;Seidel et al., 2016;Wen et al., 2022) well known to increase chlorophyll a in summer ubiquitously, and specifically in the MR system artificial reservoir waters (Fischer, 2014;Andres et al., 2019).The influence of the C6 component may be dwarfed by allochthonous inputs in free-flowing systems but is likely to be enhanced in other systems with artificial reservoirs.PARAFAC-FDOM loadings of spring samples are intermediate to those observed for fall/ winter and summer in line with it being a transition period between the three seasons (Figure 4B).Therefore, our results indicate that spring DOM pools in the MR system comprise a mixture of legacy DOM components from fall/winter and less of the presumably more labile, algal-influenced DOM characteristics of summer.The C2 component associated with potentially labile, microbial humic DOM (Murphy et al., 2014;Osburn et al., 2015;Chen et al., 2021;Smith et al., 2021) has a PC2 value near zero indicative of its consistent relative abundance regardless of season, suggesting it to be representative of the microbial community that responds to DOM sources indiscriminately.

Discharge, season, and LULC effects on PPL-isolated DOM molecular characteristics
The significant shifts of O/C (increases) and H/C (decreases) ratios from non-storm to post-storm conditions found in FT-ICR MS data (Figure 5) can similarly be explained in the context of storm-induced flow path changes and differences in soil OM along those paths.Increased contributions from MO and HO aromatic compounds, particularly in fall and winter (Figure 7), may be responsible for the O/C and H/C shifts with storm conditions.The highest average O/C and lowest average H/C are observed during fall, forested, post-storm conditions, and an inverse result (lowest average O/C, highest average H/C) is found for fall, forested, non-storm events (Supplementary Figures S5, S6).Sources of these storm-triggered MO/HO aromatic and condensed aromatic compounds include oxygenated polyphenolic plant, soil-derived lignin/tannins, and their degradation products (Siqueira et al., 1991).Lignin and tannin in leaf litter becomes degraded on timescales of weeks to months (Schofield et al., 1998;Lorenz et al., 2000;Klotzbücher et al., 2011), and their degradation products can thus be mobilized with storms throughout the year but especially during fall and winter after peak leaf fall, a time when algal production and associated export will be diminished.Diagenetic processing of OM reduces the relative contributions from lignin/tannins and oxygenated compounds in deeper soil horizons (Baldock et al., 1992(Baldock et al., , 1997;;Klotzbücher et al., 2016;Bao et al., 2017).Deeper river flow paths during non-storm conditions release DOM from organic-poor and deeper soil horizons depleted in tannin-like, lignin-like, and oxygenated compounds (Sanderman et al., 2009).Selective adsorption of aromatic hydrophobic OM with depth in soils (Kaiser andZech, 1998, 2000;Kothawala et al., 2012) may further explain the poor export of aromatic compounds from deeper soil horizons despite deep-soil OM showing higher aromaticity (e.g., Qualls and Haines, 1992;Yano et al., 2004;Hernes et al., 2007;Kothawala et al., 2012).
Time of year primarily drives the elemental composition of DOM formulae assigned by FT-ICR MS with LULC units presenting as a secondary factor.The weighted relative abundance of CHO molecular groups peaks in fall, particularly in forested lands, and drops to the lowest values during spring at the most developed site (Figure 6A).Since aromatic (mid and high oxygen levels) and condensed aromatic (low and mid oxygen levels) compounds are the most abundant in fall post-storm events (Figures 7A-C, E), the high CHO% is likely due to increasing lignin-like and tannin-like aromatic compounds derived from elevated leaf fall, plant litter, and enhanced soil release triggered by storms, particularly for forested lands (Figure 6B; ANOVA-Tukey's p < 0.05).The significant low CHO% in spring PPL-isolated DOM samples at Andrews Lake (developed LULC) is countered by the elevated contributions from N-and S-containing formulae (Figure 6; ANOVA-Tukey's p < 0.05).The higher contributions of N-and S-containing formulae in spring may be attributed to DOM interaction with inorganic nitrogen inputs from agricultural/suburban fertilizer use and/or enhanced inputs from N-and S-containing biological materials associated with new growth (i.e., DNA, proteins in phytoplankton, microbes, plants) (Mannino and Harvey, 2000;Meng et al., 2013;D'Andrilli et al., 2019).The lack of differences in elemental composition with storm status suggests that the differing water flow paths do not impact DOM composition at the coarse elemental category scale.
Seasonal factors play an important role in variations of PPL-isolated DOM FT-ICR MS characteristics in the MR system, with algal DOM muting the effects of increasing river discharge during storm events.
The winter samples show the biggest changes in molecular formula group abundances due to storm events (Figure 7) despite the winter river discharge difference between post-storm and non-storm events (5.63 m 3 /s) being less than that of the fall (15.77m 3 /s; Supplementary Figure S1).Reduced algal production in winter may render allochthonous sources more important quantitatively in the manmade pond and lake reservoirs.This reduction in algal DOM accentuates the differences in allochthonous DOM exported from deep (non-storm) and shallow (storm) flow paths with the shallow flow paths exporting compounds with higher aromaticity and oxygenation (Figures 7A-C, E) as has been observed in other work (Lambert et al., 2016;Osterholz et al., 2016;Da Silva et al., 2021).By contrast, we suggested that algal DOM produced in artificial reservoirs is quantitively important in fall, spring, and summer).
Consequently, millponds and artificial lakes can amplify seasonal effects of autochthonous production and complicate predictions of DOM composition based on an allochthonous OM-centric interpretation of the pulse-shunt concept (Figure 9).Millponds (Coursey Pond, McCauley Pond) and artificial lakes (Andrews Lake) have longer residence time compared to freeflowing streams (Spring Branch, Market Street sites) due to anthropogenic dams (Dewitt and Daiber, 1974;Andres et al., 2019), leading to the production and accumulation of algal DOM at high rates (Andres et al., 2019).Further, large amounts of nutrients (i.e., > 50 μM NO 3 ) are transported to the MR system during storms (Voynova et al., 2015) including these artificial reservoirs.Because the damming system does not operate as a true shunt process, residence times remain long enough after storms for the retained nutrients to fuel the primary productions further within these artificial reservoirs during warmer months.In fact, post-storm algal blooms are commonly observed at millponds during summer when increased irrigation dramatically reduced their water flux (Andres et al., 2019).That algal-produced DOM is characterized by relatively higher H/C and lower O/C compounds and higher heteroatom (N-, S-containing) content (Mangal et al., 2016;Wen et al., 2022), mixing with oxygenated, aromatic DOM from overland and shallow flow paths during storms.The resultant DOM mixture makes for a diversity of formula groups being represented in spring, summer, and fall, leading to the lack of statistical differences in aromatic and oxygenated formula groups between post-storm and non-storm events during those periods (Figure 7).This FT-ICR MS observation is in line with the high contribution of tryptophan/tyrosine-like FDOM in this study regardless of storm conditions (Figure 4B), providing further evidence that artificial reservoirs must be considered in the context of the pulse-shunt process.In systems without high productivity artificial reservoirs, DOM composition may be less affected by seasonal influences.
It should be noted that the seasonal differences in relative abundances of molecular formulae groups observed from FT-ICR MS data can be further affected by the seasonal patterns of soil-OM properties.The warmer temperature and enhanced microbial activity (and associated higher decomposition rates) have been suggested to increase soil-OM solubility (Tipping et al., 1999;Kalbitz et al., 2000;O'Donnell et al., 2016) and decrease its rate of sorption to soils (Kaiser et al., 2001).Increasing aromatic and condensed aromatic DOM has been observed in DOM leached from deeper soil horizons at 20 °C compared to at 5 °C or below (O'Donnell et al., 2016).The observed high contributions of aromatic and condensed aromatic compounds during spring and summer compared to those of winter during non-storm events can thus be attributed to temperature-induced changes in the mobility of soilderived DOM transported from deeper soil horizons by deep flow paths.Furthermore, the abundances of aromatic and condensed aromatic DOM from soil shallow organic-rich layers vary depending on seasonal leaf/plant-litter fall cycles.The higher contributions of high oxygenated aromatic and condensed aromatic compounds observed in fall and winter during post-storm conditions (Figure 7) are likely linked to the peak leaf fall and associated degraded products leached from shallow soil layers by storm triggered shallow river flow paths (Singh et al., 2014).Our observations imply the important influence of soil seasonal patterns on DOM molecular characteristics in combination with storm impacts.
LULC impacts can only be observed for HO Olefinic groups in combination with seasonal and storm influences (Figure 8), indicating their minor effect on DOM molecular characteristics in the MR system.Differences in the surface vegetation (including potential anthropogenic modification and maintenance; Rogerson et al., 2011) and soil chemistry across different LULC units in the MR system (Fischer, 2014) likely explains the observed changes in percentage abundances of HO Olefinic compounds (i.e., lignin-like, peptide-like compounds; Koch and Dittmar, 2006;Rivas-Ubach et al., 2018).Despite LULC units not being a primary driver of DOM molecular variations in the MR system, their significant impacts have been found in other coastal plain river system and thus should not be overlooked.Hosen et al. (2018), for example, found that forested lands release more aromatic and recalcitrant DOM while agricultural lands exports more proteinlike and bioavailable DOM during high discharge events in the Tuckahoe Creek watershed, USA.Bhattacharya and Osburn (2020) The summary schema of DOM dynamics in the Murderkill River system: (A) allochthonous inputs controlled by season, storm events, and LULC units; (B) conceptual model of DOM variations in rivers containing artificial reservoirs with season and storm events as dominant factors.The arrows direct to the highest magnitudes/values, and their sizes represents the importance of effects.Allochthonous input arrows in (B) mimic the trends described in (A).
further indicated that LULC controls the availability and mobilization of terrestrial DOM exported to the Neuse River Basin, USA, and thus primarily drives the variability in DOM composition and concentration.Therefore, LULC cannot be neglected as a factor impacting DOM characteristics in coastal plain river systems, and variations in its relative influence deserve further study.

Summary and implications for understanding riverine DOM biogeochemistry
The pulse-shunt concept proposed by Raymond et al. (2016) hypothesized that increase river discharge triggered by storm events transforms rivers from active to passive pipes and elevates the export of modern terrestrial DOM originating from surficial sources (i.e., throughfall, litter leachate, and soil water) (Inamdar et al., 2011;Barnes et al., 2018).Elaborating on the pulse-shunt concept, recent studies (i.e., Osterholz et al., 2016;Lambert et al., 2016;Da Silva et al., 2021) show that the DOM characteristics associated with high discharge events tend to be oxygen-rich, aromatic, humic-like compounds.In support of the pulse-shunt concept, our study shows that storm-triggered increases in river discharge yield increased [DOC] (Figure 2A) and export DOM that can be characterized as terrestrially derived including highly oxygenated, aromatic compounds (Figures 2,5,7).We suggest tannin, lignin, and their early diagenetic degradation products mobilized with storminduced shallow and overland flow paths to be logical sources matching those molecular characteristics.Our work elaborates on this concept by showing a seasonal impact on the composition of exported DOM that we attribute to: (1) the effects of artificial reservoirs that lengthen water residence times throughout the year including after storm-pulse events, and (2) temperature influences on soil-OM processing.These storm-induced processes and seasonal impacts are summarized in a conceptual model (Figure 9) based on our observations.This study highlights the damming effects of artificial reservoirs on the composition to DOM exported to rivers (Figure 9).We propose that the long water residence times of artificial reservoirs dampen the shunt process of pulsed materials (i.e., nutrients, terrestrial DOM) and enhance the role of seasonal patterns of autochthonous production (Figure 9) in determining DOM molecular composition.Warm weather pulses, (1) increase discharge from the reservoirs releasing bio-produced algal materials retained in artificial reservoirs, and (2) deliver nutrients to the reservoirs that further stimulate primary production and associated microbial activities (Figures 4B, 7).Our study thus has important implications for predicting how the quantities and characteristics of DOM exported from rivers to estuaries may vary with the anthropogenic installation or removal of artificial reservoirs.

FIGURE 1
FIGURE 1 Map of sampling sites in the Murderkill River watershed, Delaware (DE), United States.Artificial reservoirs are annotated with black star signs.Black dots represent sampling sites classified as draining developed lands; green squares represent sites classified as draining agricultural lands; and orange triangles represent sites classified as draining forested lands.

FIGURE 2
FIGURE 2Relationship between mean tidally filtered river discharge for the date of sample collection as measured at Frederica, DE, USGS station, and (A) DOC concentration, (B) SUVA 254 , (C) E 2 : E 3 , (D) BIX, and (E) HIX for all samples collected in the Murderkill River system over the course of this study.

FIGURE 5
FIGURE 5 Van Krevelen Diagram showing the H/C and O/C ratios of each assigned molecular formula in PPL-isolated DOM samples (n = 40).The plots above and to the right of the van Krevelen diagram are marginal histogram plots of the O/C (bandwidth = 0.015) and H/ C (bandwidth = 0.030) ratios, respectively, of assigned molecular formulae grouped by discharge conditions.The dashed lines represent the mean value of H/C and O/C ratios under non-storm (red dash line) and post-storm (blue dash line) conditions.

Figure
Figure6A).The combined percentage abundances of CHON, CHOS, and CHONS groups are only significantly different across different season and LULC conditions (ANOVA-Tukey's p < 0.05; Figure6B).The highest mean sum of CHON%, CHOS%, and CHONS% presented in spring PPL-isolated DOM samples under agricultural conditions (38.2%; Figure6B), with lowest mean values in fall forested PPL-isolated DOM samples (21.6%;Figure6B).The percentage abundances of compound categories are detailed in Supplementary TableS3(Supplementary FiguresS8-S12).Mid oxygen olefinic (MO Olefinic) formulae are the dominant class of compounds across all PPL-isolated DOM samples, ranging from 37.0% to 70.7% of total molecular formulae abundances.High oxygen olefinic (HO Olefinic) formulae are the second most abundant class of compounds, varying between 9.7% and 26.0% of total molecular formula abundances across different PPL-isolated DOM samples.Mid oxygen aliphatic (MO Aliphatic), mid oxygen aromatic (MO Aromatic), and low oxygen olefinic (LO Olefinic) formulae are also well represented in PPL-isolated DOM samples

FIGURE 6
FIGURE 6Boxplots with the solid horizontal line indicating the median percentage abundances for all PPL-isolated DOM samples (n = 40), with lower 25th and upper 75th percentile, of (A) CHO and (B) sum of CHON, CHOS, and CHONS groups in spring (white), summer (blue), fall (orange), and winter (gray) (n = 10) for each season) classified by LULC conditions.Outlier values are > 1.5 times the interquartile range beyond the upper end of the box.The mean percentage abundances are labeled in the black numbers below the boxes.Letters indicate results of Tukey's honestly significant differences post hoc multiple mean comparisons (ANOVA-Tukey's p < 0.05).Fall statistics are presented in Supplementary TablesS9, S10.

FIGURE 7
FIGURE 7 Boxplots indicating the median percentage abundances of (A) low oxygen condensed aromatic (LO CondAromatic), (B) mid oxygen aromatic (MO Aromatic), (C) mid oxygen condensed aromatic (MO CondAromatic), (D) high oxygen aliphatic (HO Aliphatic), and (E) high oxygen aromatic (HO Aromatic) compound classes in spring (white), summer (blue), fall (orange), and winter (gray) across all PPL-isolated DOM samples (n = 40) grouped by different discharge classifications (n = 5 for each season-discharge status pair).The lower 25th and upper 75th percentile are represented by the upper and lower borders of the boxes, and outlier values are > 1.5 times and/or < 3 times the interquartile range beyond either end of the box.The mean values within each identified condition are labeled in black numbers below the boxes.Statistically significant differences of percentage abundances for the individual compound class under different conditions are indicated by letters at ANOVA-Tukey's p < 0.05.Full statistics are presented in Supplementary TableS11.

FIGURE 8
FIGURE 8 Mean percentage abundances (labeled in black numbers) of the high oxygen olefinic (HO Olefinic) compound class in four seasons (n = 10 for each season) across all PPL-isolated DOM samples (n = 40) categorized by different discharge and LULC conditions.Error bars represent the highest and lowest values within each identified category.Letters show the results from the Tukey's honestly significant different test with ANOVA-Tukey's < 0.05.Full statistics are present in Supplementary TableS12.

TABLE 1
Excitation and emission maxima, characterizations of modeled PARAFAC components, and related references matched from the OpenFluor spectra database with similarity scores specified.