Implications of sample treatment on characterization of riverine dissolved organic matter †

High-resolution mass spectrometry techniques are widely used in the environmental sciences to characterize natural organic matter and, when utilizing these instruments, researchers must make multiple decisions regarding sample pre-treatment and the instrument ionization mode. To identify how these choices alter organic matter characterization and resulting conclusions, we analyzed a collection of 17 riverine samples from East River, CO (USA) under four PPL-based Solid Phase Extraction (SPE) treatment and electrospray ionization polarity ( e.g. , positive and negative) combinations: SPE (+), SPE ( (cid:1) ), non-SPE ( (cid:1) ), and non-SPE (+). The greatest number of formula assignments were achieved with SPE-treated samples due to the removal of compounds that could interfere with ionization. Furthermore, the SPE ( (cid:1) ) treatment captured the most formulas across the widest chemical compound diversity. In addition to a reduced number of assigned formulas, the non-SPE datasets resulted in altered thermodynamic interpretations that could cascade into incomplete assumptions about the availability of organic matter pools for heterotrophic microbial respiration. Thus, we infer that the SPE ( (cid:1) ) treatment is the best single method for characterizing environmental organic matter pools unless the focus is on lipid-like compounds, in which case we recommend a combination of SPE ( (cid:1) ) and SPE (+) to adequately characterize these molecules.


Introduction
The development and increasing utilization of high-resolution mass spectrometry (HRMS) techniques (e.g., FTICR-MS, Orbitrap-MS, IMS-MS) has allowed scientists to characterize the molecular compounds that constitute natural organic matter 1 (NOM). These approaches have been applied to samples across a diverse range of environments, [2][3][4][5] and have yielded new insights into NOM processing. More recently, attempts have been made to couple these high-resolution chemical analyses with microbiome data to directly link measurements of microbial function with changes in chemical species. [6][7][8] One of these approaches involves the analysis of dissolved organic matter (DOM) in a thermodynamic framework through which the nominal oxidation state of carbon (NOSC) can be calculated. Specically, an empirical relationship between the NOSC and the standard molar Gibbs free energies of the oxidation half reactions of organic compounds 9 allows researchers to quantify the thermodynamic favorability of DOM pools as electron donors for microbial respiration from HRMS data. One study using this approach demonstrated that aerobic respiration increased with increasing DOM thermodynamic favorability under carbon limiting conditions. 6 Similarly, a study of carbon dynamics within an anoxic aquifer revealed the preservation of DOM pools that yielded insufficient energy when coupled to the reduction of sulfate. 10 Other recent studies taking a similar approach have further revealed relationships between DOM thermodynamics and various ecological measurements, including biogeochemical function and microbial community assembly. 7,[11][12][13] There are many pre-treatment and instrument options to choose when utilizing HRMS techniques that could bias the resulting dataset. Environmental samples can be analyzed without enrichment or treatment but are commonly pre-treated using PPL-based solid phase extraction (SPE) to concentrate DOM or remove species that may interfere with the ionization of organic compounds 14,15 (e.g., salts). During SPE, DOM is retained on a sorbent then eluted using a solvent such as methanol. Polymer-based sorbents (e.g., PPL) are commonly used during SPE due to their high extraction efficiency and representative character of retained DOM. 16,17 However, DOM recovery by SPE is incomplete and studies have shown that single-sorbent SPE methods can result in the preferential loss of specic compound classes such as organo-sulfur compounds, 18 aliphatic amines/amides, and tannin-like compounds with a high oxygen content. 19,20 Electrospray ionization (ESI) is a common ionization technique in HRMS because it can ionize a wide range of polar, hydrophilic molecules with diverse functional groups common in DOM, results in minimal-to-no fragmentation of DOM molecules (i.e., is a 'so' ionization technique), and can be used in either positive or negative mode. 21 DOM can contain acidic functional groups that can be readily deprotonated in ESI (À) or basic functional groups (e.g., amines) that can be readily protonated and ionized efficiently in ESI (+). Thus, the chosen ESI mode can further govern the nal HRMS DOM spectra and the chosen ESI mode has been found to result in different mass spectra for the same sample. 22,23 Ohno et al. (2016) 24 found that ESI (+) was better at detecting aliphatic and carbohydrate-like molecules and that the ionization of these molecules was generally suppressed in samples rich in aromatics ionized in ESI (À) mode. The selected ESI mode can further impact downstream data through the generation of unintended adducts which can confound molecular formula assignment. For example, ESI (À) can result in the formation of chloride-containing ions due to salinity interference. 16 Finally, Hawkes et al. (2020) 25 performed a cross-lab comparability study of HRMS data from different instruments to provide valuable metrics for future data benchmarking and further revealed patterns based upon ionization modes.
Although other studies have analyzed impacts of instrumentation, sample preparation, or ionization modes, no study has systematically examined the combined impacts of SPE and ionization mode choices or discussed changes to thermodynamic or microbial metabolic interpretations. Here, we analyze a set of 17 surface and pore water samples from the East River, CO, using four different Fourier transform ion cyclotron resonance mass spectrometry (FTICR-MS) methods (SPE or non-SPE, positive-ion or negative-ion ESI mode) with the goal of providing guidance to researchers on what method is most appropriate for their research question. While it is difficult to completely standardize these treatments (particularly regarding differences arising due to SPE), we have selected these because they represent common routes of DOM analyses in complex mixtures. We show that both sample pre-treatment and ESI mode choice inuences the detection of different DOM molecules and can therefore bias the interpretation of microbial DOM processing within a given sample.

Sample collection
Surface and pore water samples (2 surface, 15 pore) were collected from the East River (CO, USA) on a 200 m reach encompassing one meander that lies within the U.S. Department of Energy-supported Lawrence Berkeley National Laboratory's Watershed Science Focus Area (Fig. S1 †). This combination of 17 previously-analyzed samples were selected because they had representatives across all treatment types; see  27 for more detailed information on sampling efforts and information on the remainder of the samples. Briey, pore water samples were collected from their respective depths (Table S1 †) using a 0.6 cm diameter, stainless steel pore water sipper with a screen length of 5 cm attached to a syringe (MHE Products, MI, USA). One tubing volume of approximately 30 mL was discarded before sampling and sample was immediately ltered through 0.22 mm Sterivex lters (housings made of Eastar co-polyester; Massachusetts, USA). The 20 mL aliquot for the data presented here was ltered into a pre-combusted amber glass vial and immediately placed in a cooler on ice for transport back to the lab.  Table S1 †.

Sample preparation and ESI-FTICR-MS data collection
Fourier Transform Ion Cyclotron Resonance Mass Spectrometry (FTICR-MS) was used to provide ultra-high resolution organic matter characterization. For the samples undergoing solid phase extraction, aqueous samples (NPOC 0.33-0.99 mg C per L) were acidied to pH 2 with 85% phosphoric acid and extracted with Bond Elut PPL cartridges (Agilent, Santa Clara, CA), following Dittmar et al. (2008). 16 For those samples which were not extracted, 250 mL of sample were mixed with 500 mL of LC-MS grade MeOH. A 12 Tesla (12T) Bruker SolariX Fourier transform ion cyclotron resonance mass spectrometer (Bruker, Billerica, MA) located at the Environmental Molecular Sciences Laboratory in Richland, WA was used to collect high-resolution mass spectra of the organic matter found in each sample. Samples were directly injected into the instrument using a custom automated direct infusion cart that performed two offline blanks between each sample. 28 A Bruker SolariX electrospray ionization (ESI) source was used in positive and negative modes with applied voltages of +4.4 kV and À4.2 kV, respectively. Ion accumulation time was optimized between 50 and 80 ms. One hundred and forty-four transients were coadded into a 4MWord time domain (transient length of 1.1 s) with a spectral mass window of m/z 100-900, yielding an average resolution of 265 K at m/z 400. Spectra were internally recalibrated in the mass domain using homologous series separated by 14 Da (CH 2 groups) (File S1 †). The mass measurement accuracy was typically within 1 ppm for singly charged ions across the mass window. Bruker Daltonics DataAnalysis (version 4.2) was used to convert mass spectra to a list of m/z values by applying the FTMS peak picking module with a signal-to-noise ratio (S/N) threshold set to 7 and absolute intensity threshold to the default value of 100%. Formularity 29 was used to align peaks using a 0.5 ppm threshold, and then assign chemical formulas based on exact mass with a mass measurement error <0.5 ppm while allowing for CHONS (Formularity formula settings: O > 0 AND (N + S) < 6 AND S < 3 AND P ¼ 0).
The R package msRanalysis 30,31 was used to process the report generated by Formularity and remove peaks that either were outside the desired m/z range (150-900 m/z) or had an isotopic signature, calculate nominal oxidation state of carbon (NOSC), assign putative compound classes, 32 and organize the data. We have included a table that describes the number of assigned peaks included aer each ltering step (  25 was downloaded in R from https://github.com/ BarrowResearchGroup/InterLabStudy using the provided R package. The Hawkes et al. dataset consists of molecular formulas obtained by analyzing four different organic matter standards (Elliot Soil Fulvic Acid, Pony Lake Fulvic Acid, Suwannee River Fulvic Acid, Suwannee River Natural Organic Matter) with either FTICR-MS or Orbitrap-MS under varied ionization modes. NOSC values were calculated and compound classes were assigned to each molecular formula in order to evaluate whether patterns observed using data from this manuscript were broadly transferrable.

Thermodynamics calculations
We calculated the average thermodynamic potential factor (F T ) for the oxidation of average DOM pools coupled to the reduction of O 2 and SO 4 2À at standard conditions. 9,33 The F T was calculated as follows: where DG r is the Gibbs energy of the reaction, m is the number of moles of ATP synthesized per formula reaction (0.15 and 2 used for aerobic respiration and SO 4 2À reduction, respectively 33 ), DG ATP denotes the Gibbs energy to synthesize ATP (50 kJ mol À1 used here 33 ), c is the average stoichiometric number for the reaction of interest, and R and T are the universal gas constant and temperature, respectively. To determine c for our complex DOM pools, we averaged the values used for different compounds in the analysis presented in LaRowe and Van Cappellan (2011) 9 (i.e., averaging the stoichiometric numbers used for amino acids, mononucleotides, saccharides, complex organics, etc.).

KEGG mapping
In order to evaluate differences in potential biogeochemical interpretations, we assigned KEGG compound identiers (CPD numbers) to observed molecular formulas by mapping them to the KEGG database using the provided REST API. 34 In order to avoid potential duplication of CPD numbers, only exact formula matches were considered. Using these CPD numbers, we iden-tied corresponding pathways and qualitatively compared the most represented pathways across treatments. Web scraping was performed using the R package rvest 35 and the KEGG mapping scripts are available on GitHub at https://github.com/ danczakre/FTICR-Methods-Comparison.

Statistics and plot generation
The statistics program R was used to perform all statistical analyses with the R package ggplot2 used to generate all plots. 36

Molecular formula detection signicantly varies based on sample preparation and ionization type
We evaluated 17 samples that were processed with or without SPE and analyzed in ESI (+) or ESI (À) mode to investigate the impact of method selection on FTICR-MS data. Initial multivariate analyses revealed clear differences in the collected data between methods (Fig. S2 †). Considering the molecular formula count for each treatment, the SPE (À) treatment yielded the most assigned formulas (6218) while non-SPE (À) yielded the fewest (1257) ( Table S2 †). More broadly, the SPE-treated samples contained greater molecular formula counts than non-SPE samples (Fig. 1) and shared more common formulas (Table S1 †). These results are expected because SPE is a common method to concentrate DOM in a sample and reduce the impact of salt during ionization. 18,39 . By ensuring that carbon concentrations are higher, and that salts have a smaller impact on ionization efficiency, more molecular formula will be assigned. Differences in formula count between the two ionization modes are likely impacted by differences in the types of compounds ionized; specically, river corridor organic matter is typically acidic and oen has functional groups rich in oxygen which renders it more likely to be detected using negative mode. 25,40 The noted increase in molecular formula observed when comparing the non-SPE (À) vs. non-SPE (+), however, is likely the result of the SPE selecting for compounds more readily observed in negative mode. In other words, the molecular formulas observed under non-SPE conditions are more readily detectable in positive-mode or experience less ionization competition due to salinity. Additionally, the acidication of samples prior to SPE would result in more compounds readily detectable in ESI (À) relative to ESI (+) and basic modication, rather than acidication, prior to SPE may enhance the number of identiable compounds in ESI (+). 41 In these samples, conductivity (used as a proxy for salinity) ranged from 267.8-387.4 mS (Table S1 †). PPL-based solid-phase extraction can also lead to compositional shis in the detected DOM given that some compound types have a higher affinity for the column (e.g., sulfur-containing, hydrophobic, low O/C compounds). 18,39 The affinity of sulfur-containing compounds for the PPLcolumn is apparent in our dataset; SPE retention resulted in more assigned CHONS formulas in non-SPE in both negative and positive ionization mode (7.8% and 18.1%, respectively).
We observed the largest variation the proportion of formulas assigned to be CHO-only, with values ranging from 10.2% in non-SPE (+) to 49.3% in non-SPE (À) and signicant differences in proportion between each condition (Mann-Whitney U pvalue ( 0.001) (Fig. 1a). The SPE (À) samples had a similar proportion of CHO-only molecular formula as the non-SPE (À) samples (40.6% and 49.3%, respectively; Mann-Whitney U pvalue: 0.392), whereas the ESI (+) samples diverged in CHO-only formula proportions (10.2% and 29.7% in non-SPE and SPE, respectively). The large proportion of CHO-only formula in ESI (À) are consistent with observations that ESI (À) is primarily used to detect O-rich molecular formula. 25,40 CHON-and CHOScontaining formulas also demonstrated variability across treatment conditions overall (ranging from 14.4% to 32.2% and 14.9% to 27.7%, respectively; Fig. 1a), while sample-resolved analyses reveal more complex patterns. ESI (+) samples, regardless of extraction method, consistently contained a higher proportion of CHOS-containing formula than ESI (À) samples (Mann-Whitney U test p-value ( 0.001). These patterns reect the potential for ESI (+) to capture more compositionally complex spectra than ESI (À), 23 whereas ESI (À) is more ideally suited for capturing CHO-only formula. Because environmental DOM is enriched by acidic, O-rich compounds, 24,40,42 the ESI (À) approach likely results in better characterization of the riverine DOM pool studied here.
Van Krevelen-based compound classes exhibited high variability across treatments as well (Fig. 1b). While the non-SPE (+) treatment yielded more molecular formulas than non-SPE (À), both non-SPE ionization modes resulted in similar proportions of compound classes when all samples are considered together. The SPE (À) treatment generated the highest proportion of lignin-and tannin-like formulas (43.3% and 9%, respectively) while the SPE (+) treatment yielded a high proportion of lipidlike compounds (31.9%; Fig. 1b). Sample-resolved analyses further showed these trends (Table 1). Non-SPE (À) samples had a higher proportion of protein-like formulas than the other treatments (Mann-Whitney U test p-value < 0.001) and non-SPE (+) had signicantly greater representation of concentrated hydrocarbon-like formulas (Mann-Whitney U test p-value ( 0.001). Lignin-like and tannin-like formulas dominated SPE (À) samples (Mann-Whitney p-value ( 0.001) and both non-SPE (À) and SPE (+) had increased lipid-like efficiency (Mann-Whitney p-value < 0.001). Given that both the lignin-and tannin-like compound classes are characterized by higher O : C ratios (>0.28 and >0.65 respectively), we argue that these patterns reect the enhanced capability the PPL cartridge to retain O-rich compounds. 32,39 This suggests that solid phase extraction has some combinatorial effect when used in conjunction with ESI (À) due to the elevated potential for ESI (À) to enrich for O-rich compounds, as we observed above. It is also important to note that the acidication of sample prior to SPE also likely selects for compounds more readily identiable in SPE (À).

Thermodynamics vary substantially based upon extraction methodology and ionization mode
The NOSC for a given sample was used to evaluate the potential thermodynamic implications of the methodological differences. The NOSC metric can reveal the potential thermodynamic favorability of a carbon substrate (or bulk DOM pool favorability when averaged together), with higher NOSC values theoretically yielding a lower overall DG C ox (i.e., more favorable) when coupled to the reduction of an electron acceptor. 9 We observed that each method yielded NOSC values signicantly different from each other method (p < 0.001) and that compounds detected in the non-SPE (+) treatment had higher NOSC values overall while the average NOSC for detected compounds in the SPE (+) was signicantly lower (Fig. 2). Interestingly, the impacts of SPE on NOSC within each ionization mode are ipped; in ESI (+) mode, the SPE samples have signicantly lower NOSC values whereas, in ESI (À) mode, the non-SPE samples have signicantly lower NOSC values. We hypothesize that this could be due to potentially driven by the higher average proportion of lipid-like molecular formulas in both of these treatments, which have historically had lower NOSC values. 10 Thus, non-SPE (+) processing of these samples would have yielded DOM compounds considered less thermodynamically favorable in many biogeochemical analyses. 6,10,43 To evaluate the transferability of observations derived from our dataset, we performed similar analyses on data previously utilized in Hawkes et al. (2020). 25 In brief, these data consisted of four separate organic matter standards (ESFA: Elliot Soil Fulvic Acid, PLFA: Pony Lake Fulvic Acid, SRFA: Suwannee River Fulvic Acid, SRNOM: Suwannee River Natural Organic Matter) collected across international high-resolution mass spectrometry instruments (e.g., FTICR-MS, Orbitrap-MS) in both ESI (+) and ESI (À) modes. The fulvic acid samples were prepared by isolating DOM onto XAD-8 resin and the SRDOM sample was  Fig. S3 †). The results help assert that the impact of ionization mode is similar across sample pre-treatment methods and sample types (e.g., ESFA is a soil standard vs. our freshwater samples).
To further explore the thermodynamic implications of the FTICR-MS methods, we calculated the average thermodynamic potential factor (F T ) for the oxidation of the average DOM pools coupled to reduction of O 2 and SO 4 2À at standard conditions ( Fig. 3). Briey, F T is a dimensionless value derived from transition state theory that couples the rate of reaction to the Gibbs free energy of reaction. A value of 1 indicates that bioenergetic limitation is ignored, the reaction is more kinetically controlled, and, if there is substrate available, there will likely be a positive reaction rate until substrate is consumed. 9,33 Thus, F T values are helpful for identifying differences in potential organic matter degradation rates between different environments or samples. Data resulting from the non-SPE (+) treatment consistently yielded the highest F T values for reduction of both O 2 and SO 4 2À , while the F T values calculated from SPE (+) analyses were generally the lowest of the four treatment conditions (Fig. 3). This indicates that the suite of compounds detected by non-SPE (+) would have the fewest thermodynamic constraints if utilized for sulfate reduction and that oxidation of the compounds detected by SPE (+) have the highest thermodynamic constraint when coupled to aerobic respiration. The compounds detected in SPE (À) consistently had high F T values when coupled to aerobic oxidation and had a large range for sulfate reduction (0.02-0.51; Fig. 3). As we observed with the NOSC values, we suggest that this is a result of non-SPE (À) and SPE (+) selecting for larger proportions of comparably thermodynamically unfavorable high H : C, low O : C formulas (e.g., lipid-like). As this approach is commonly coupled to FTICR-MS data, 44 it is important to understand how the pre-treatment and instrument method inuences nal F T calculations and further biogeochemical interpretations.

Divergent ecosystem interpretations arise due to extraction methodology and ionization mode
We observed variability in elemental and lignin composition across ESI AE and can leverage this to further evaluate the organic matter detected by each ionization mode. For example, lignin : N and C : N ratios have been used as an index of organic matter stability and as a means to estimate decomposition rate. 45,46 We calculated these widely used ratios for each of the datasets and found that the methods study here greatly impact these values (Fig. S4 †). The non-SPE (+) dataset had  a signicantly lower C : N (here, organic C to organic N) than each other dataset and the SPE (À) dataset had a signicantly higher lignin : N ratio due to the high detection of lignin-like compounds and no acidication prior to analysis (Fig. 1). This would lead to an interpretation that the DOM pool detected in non-SPE (+) has a higher residence time than those detected in the other methods. The higher lignin : N ratio in the SPE (À) dataset cascades into an interpretation where we would assume lower decomposition rates in this dataset relative to the others.
To identify whether the trends seen in the bulk dataset held for sample-resolved analyses, we chose a representative sample for direct comparison of selected measurements between the four methodologies (Table 1). A single representative sample was selected to focus specically on the thermodynamic relationships across the treatment methodologies in one physical location, rather than broad scale biogeochemical processes explored elsewhere. 26 Overall, there were clear differences in the formula count, with SPE (À) still resulting in the most assigned formulas (1971). Additionally, despite variability in NOSC values across the sample treatments, measured DOM pools for this sample have few to no thermodynamic restraints for the reduction of O 2 (i.e., F T values close to 1 - Fig. 3a). Conversely, these same treatment-derived differences in NOSC result in greater variability in the F T parameter when oxidation of DOM is coupled to reduction of SO 4 2À . All the DOM pools detected using the SPE (+) and non-SPE (À) treatment yielded an F T value of 0, indicating that sulfate reduction would be thermodynamically inhibited. In contrast, the non-SPE (+) produced the highest F T when coupled to SO 4 2À reduction, indicating fewer thermodynamic constraints. Lastly, similar to the bulk dataset analysis (Fig. S4 †), the ratio of lignin-like compounds to N (lignin : N) was signicantly higher in the SPE (À) dataset due to the large number of lignin-like compounds detected (Fig. 1B) and the C : N ratio was still the lowest in the non-SPE (+). Thus, this example further illustrates how DOM treatment can inuence our interpretation of likely redox reactions and potential decomposition rates occurring in each sample. Given the observed patterns of unique and shared molecular formula across methodologies, extraction method and ionization mode have signicant implications for biogeochemical modeling. In the case of recent multi-omics analysis, environmental metabolite data collected via high resolution MS can be incorporated into pathway models to understand carbon ux or organismal metabolism. 47 Differences in the types of formula detected during data collection could have profound impacts on inferred metabolic pathways. Given the low percentage of shared molecular formula across methods (Table S3 †), detected pathways may be divergent resulting in different predicted paths for carbon ux. Using the data presented within this manuscript, for example, we observe variation in the non-specic "Metabolic pathways" and "Biosynthesis of secondary metabolites" categories across the methods with SPE (À) having the highest absolute abundance and non-SPE (+) having the lowest (Fig.S6 †). We also observe pathways unique to methods; for example, "Degradation of aromatic compounds" and "Flavonoid biosynthesis" were only present in SPE (À) data while we only saw "Limonene and pinene degradation" in SPE (+) (Fig. S6 †). Other types of biogeochemical analyses, such as those relying on substrate-explicit modeling, will also be signicantly impacted due to the shi in understood thermodynamic availability of each compound. 12 Specically, the analysis proposed by Song et al. (2020) 12 relies heavily on the predicted Gibb's free energy of various carbon reactions to estimate reaction kinetics and stoichiometric coefficients of catabolic and anabolic reactions. Shis in the average NOSC values will cascade into divergent predictions of biogeochemical rates.

Conclusions
We achieved the highest number of assigned formulas (11 072) and representation from all molecular and van Krevelen-based compound classes when combining the four datasets (Fig. S5 †). Although combining methods would most holistically characterize DOM pools, experimental design clearly depends on (1) instrument availability, (2) available sample volume, and (3) funding resources. We therefore provide some of the following recommendations to help researchers prioritize their desired outputs. Our analyses rst conrm that SPE treatment yields the largest number of molecular formulas from these riverine DOM samples and that ESI (À) captures the most formulas and is relatively sufficient for capturing broad chemical diversity (Fig. 1). However, if lipid-like compounds (e.g., formulas with low O : C and high H : C) are the overall target, extracted samples run with ESI (+) would be a better recommendation. If the aim is a thermodynamic analysis of the DOM pools, SPE (+) appears to detect compounds with signicantly lower NOSC values (Fig. 2) than the other methods, potentially altering downstream thermodynamic predictions. We therefore recommend that SPE (À) or SPE (À) combined with SPE (+) is used for most applications. Non-SPE yields the lowest number of formulas and has the most limited applications due to salinity and other interferences. While it may be used if the laboratory has either time or funding constraints, we highlight the potential for saline interference in data capture. If samples are not treated with SPE, we recommend the use of ESI (+) instead of ESI (À) because it not only yields more formulas but also results in more formulas from each Van-Krevelen-based compound class, likely resulting in a more holistic characterization of the DOM pool (Fig. 1). However, if the study is focusing on thermodynamics, ESI (À) detects compounds more representative of the DOM pool (Fig. 2). These recommendations should be taken into consideration when performing experimental design. Furthermore, given that we observe signicant differences in methodological decisions, we believe that more systematic comparisons across methods (e.g., different ionization techniques, extractions, storages, etc.) are needed to aid scientists in developing well-designed experiments utilizing FTICR-MS.

Conflicts of interest
There are no conicts of interest to declare.