High rates of daytime river metabolism are an underestimated component of carbon cycling

Study sites and data collection

During 2017 and 2018, we carried out 14 experiments in rivers located in temperate, tropical, and subarctic biomes to capture a gradient of river productivity and climatic characteristics (Table 1, Fig. 1). Apart from the Mekong and Sekong rivers in Cambodia that were impacted by plantations, rice cultivation, grassland, and urban areas (56% impacted land cover in the Mekong and 38% in the Sekong), the selected rivers were predominantly in pristine areas (impacted land-use ≤ 8%), although two rivers in Mongolia were affected by livestock grazing (with 26% of land cover at the Khovd and 59% in the two Zavkhan rivers).

We conducted traditional O₂ concentration metabolic assessments, assessments of isotopic fractionation, and 24 h characterization of δ¹⁸O₂ at each site_. We measured changes in dissolved O₂ concentrations and temperature every 10 min over at least 24 h with at least one MiniDOT logger (PME, Vista, California, USA). We calibrated for drift using the average measurement values made in 100% saturated water for at least 30 min before and after each deployment to allow adjustment to temperature and placed sensors in the river for at least 30 min prior to using data to allow equilibration to temperature (following methods detailed in ref. ⁵²).

We collected δ¹⁸O₂ samples by hand every 2 h during the same 24-h period of the O₂ concentration measurements in pre-evacuated 100 mL vials loaded with 50 µl HgCl₂ as a preservative and sealed with septum stoppers (Bellco Glass Inc., Supelco, Vineland NJ). We analyzed samples for δ¹⁸O₂ at the Nevada Stable Isotope Lab of the University of Nevada, Reno with a Micromass Isoprime (Middlewich, UK) stable isotope ratio mass spectrometer. We followed the method described by ref. ¹⁷ and injected 1.0–2.5 mL of headspace gas taken from the serum bottles using a gastight syringe (SGE, Australia) into a Eurovector (Pavia, Italy) elemental analyzer equipped with a septum injector port, and a 1.5 m long molecular sieve gas chromatography column. Water-δ¹⁸O was also collected at each site every 2 h and analyses were performed using a Picarro L2130-i cavity ringdown spectrometer at the Nevada Stable Isotope Lab of the University of Nevada, Reno. δ¹⁸O₂ values are reported in the usual δ notation vs. VSMOW in units of ‰, with an analytical uncertainty of ±0.2‰ for δ¹⁸O₂, or an analytical uncertainty of ±0.1‰ for water-δ¹⁸O.

We characterized physical characteristics at each site to provide parameters to estimate whole-system metabolism. We measured conductivity, slope, and flow velocity and depth at ten transects using a flow meter when wadeable or with an Acoustic Doppler Velocimeter (Sontek, Xylem, San Diego, CA) when rivers were not wadeable. At each site, we measured light as photosynthetically active radiation (PAR) every 10 min, using Odyssey PAR loggers (Data Flow Systems, Christchurch, New Zealand) calibrated with a Li-Cor PAR sensor (Lincoln, Nebraska, USA).

At each site, we also directly measured biofilm ash-free dry mass (AFDM) from 8 to 12 rocks (⁵³). The material was scrubbed from the rocks, agitated, filtered (Whatman glass microfiber GF/F filters). Rock area was estimated with calibrated pictures processed with the ImageJ processing program (National Institutes of Health and the Laboratory for Optical and Computational Instrumentation LOCI, University of Wisconsin). For AFDM analyses, samples were dried, and weighed before and after combustion.

Additionally, we collected data on the percentage of impacted land use in the watershed above each sampling site: for the Mekong and the Sekong we used Landsat satellite imagery from ref. ⁵⁴, for the US and Mongolian sites land use characteristics were derived from the National Land Cover Database⁵⁵ and for Patagonia we used the Chilean national land use inventory maps from ref. ⁵⁶.

δ¹⁸O₂ stable isotope fractionation during respiration in sealed recirculating chambers

Models based on oxygen isotopes are sensitive to the oxygen isotope fractionation factor (α_R) during respiration used; α_R can vary widely among sites and is influenced by temperature and water velocity³⁰. We used in our models the range of α_R values measured by³⁰ using sealed Plexiglas recirculating chambers as in ref. ⁵⁷. These measurements were done at the same time as the 24 h δ¹⁸O₂ sample collections in the rivers of this study. We placed rocks, sediment, macrophytes (macrophytes dominated in the Zavkhan 1 site) inside the chambers, depending on the site’s dominant substrata (see ref. ³⁰ for more details on chamber measurements). We collected water samples in the chambers for δ¹⁸O₂ analyses before and after the incubations and the O₂ isotope fractionation factor was calculated using Eq. (2).

where δ is the O₂ isotopic composition of dissolved oxygen at the end of the dark incubation, δ_i is the O₂ isotopic composition of dissolved oxygen at the beginning of the dark incubation, F the fractional abundance of O₂ concentration remaining at the end of the dark incubation, and α is the isotopic fractionation factor during respiration.

Ecosystem metabolism O₂ single station modeling

We modeled metabolism as a function of GPP, ER, and reaeration with the atmosphere, using the single-station open-channel metabolism method⁴ using the same approach as¹⁵, given in Eq. (3).

where GPP is gross primary production in g O₂ m⁻² d⁻¹, ER is ecosystem respiration in g O₂ m⁻² d⁻¹, ({K}_{{O}_{2}}) is the reaeration coefficient (d⁻¹). PPFD is photosynthetic photon flux density (µmol m⁻² s⁻¹), z is mean stream depth (m), and ∆t is time increment between logging intervals (d). We used Bayesian inverse modeling approach to estimate the probability distribution of parameters GPP and ER that produce the best model fit between observed and modeled O₂ data. We fixed site-specific ({K}_{{O}_{2}}) estimates using K600 (d⁻¹) (normalized beyond gas-specific Schmidt number conversions among gases⁵⁸) based on prior work characterizing K using BASE⁵⁹, and converted these prior estimates of K600 to ({K}_{{O}_{2}})using appropriate temperature corrections. We estimated daily GPP and ER from diel O₂ data only (Eq. (3)) to be used as prior estimates of daily GPP_O2 and ER_O2 in the coupled O₂ and δ¹⁸O₂ model (Eqs. (4a) and (4b))¹⁵, where the mean and SD of GPP and ER from the O₂ ^_only method were used as prior estimates of GPP_O2 and ER_O2 in the dual O₂ and δ¹⁸O₂ model described below.

Ecosystem metabolism: Diel δ¹⁸O₂ modeling

We also modeled metabolism using an updated version of the model developed by ref. ¹⁵ coupling high-frequency O₂ concentration data with δ¹⁸O₂ collected every 2 h throughout the same 24 h period of the O₂ concentration measurements. With this model, daily rates of ecosystem metabolism are derived from diel changes in δ¹⁸O₂ and O₂, where values of δ¹⁸O₂ are converted to g ¹⁸O m⁻³ (¹⁸O₂ in Eq. 4b) and modeled as a function of water isotope values, isotope fractionation, reaeration with the atmosphere, ER, and GPP. As with Eq. 3, the ratio of light at the previous logging time (({{PPFD}}_{left(t-1right)})) relative to the sum of light over 24 h (({sum {PPFD}}_{24h})) is used to characterize times when GPP is zero and only ER is taking place (Eqs. (4a) and (4b)):

Where GPP_O2 and ER_O2 (g O₂ m⁻² d⁻¹) refer to the values obtained from diel O₂ only, dielMET (g O₂ m⁻² d⁻¹) is the diel metabolism term that allows for the estimation of diel ER and GPP from ¹⁸O₂, K_O2 is the O₂ gas exchange rate (d⁻¹), z is mean stream depth (m), PPFD is photosynthetic photon flux density (µmol m⁻² s⁻¹), Δt is time step between measurements (d), ¹⁸O₂ is the concentration of ¹⁸O in dissolved O₂ (g ¹⁸O m⁻³), AF_DO is atomic fraction of dissolved O₂ (mol¹⁸O:mol O_2, measured)_, AF_w is atomic fraction of H₂O (mol ¹⁸O:mol O_2, measured), AF_atm is atomic fraction of atmospheric air (mol¹⁸O:mol O_2, literature), α_g is the fractionation factor during air–water gas exchange (0.9972, from ref. ⁶⁰), α_R is the fractionation factor during respiration measured in the chambers (varied by site³⁰; Fig. 1), α_p is the fractionation factor during photosynthesis (1.0000 from ref. ⁶⁰).

The inverse modeling approach finds the best estimates of parameters to match measured and modeled dissolved O₂. The model assumes that the measured changes in O₂ concentration represent the actual net diel changes in O₂ concentration and uses an additional parameter, dielMET, that is a function of the isotopic enrichment occurring during respiration, derived from diel ¹⁸O₂. This parameter increases daily ER_O2 and GPP_O2 of the same amount, adding and subtracting dielMET, to obtain daily δ¹⁸O_2-ER and δ¹⁸O₂-GPP, respectively.

We estimated the posterior distributions of unknown parameters (ER_O2, GPP_O2, and dielMET) using a Bayesian inverse modeling approach¹⁵ and Markov chain Monte Carlo sampling with the R metrop function in the mcmc package^61,62. Each model was run for at least 200,000 iterations using nominally informative priors based on the range of ER_O2 and GPP_O2. For dielMET, we used a minimally informative uniform prior distribution (0–100 g O₂ m⁻² d⁻¹). We removed the first 10,000 iterations of model burn-in and assessed quality of model fit. Model runs using the minimum, average, and maximum α_R values measured in the field recirculating chambers were also compared, and we selected the α_R and report associated model metabolism estimates that generated the lowest sum of squared differences between the observed and modeled O₂ and ¹⁸O₂ diel values.

Temperature-normalized comparisons

To test the effect of temperature from the daily δ¹⁸O₂-ER and δ¹⁸O₂-GPP rates and account for daily variations in temperature, we normalized estimates from models to 20 °C (and report them as 20δ¹⁸O₂-ER and 20δ¹⁸O₂-GPP) for comparison with O₂-derived metabolism estimates following³³ with Eq. (5):

Where t₁ is site temperature and rate is the measured rate (i.e., GPP or ER) at t₁.

Statistical analyses

We used multiple linear regression to find the best predictor of the magnitude of diel 20δ¹⁸O₂-ER and differences between sites. To select the best model, we performed a stepwise variable selection and selected the best model based on the lowest AIC. Tested variables included percentage of impacted land use (%), 20δ¹⁸O₂-GPP (g O₂ m⁻² d⁻¹), conductivity (µS/cm), ash-free dry mass (AFDM, g), slope (%), water depth (m), and flow velocity (m/s) measured in the field. We used ANOVA to test the relative contribution of each variable selected with the AIC to total variance. Analyses were run with the R software⁶¹.

Reporting summary

Further information on research design is available in the Nature Portfolio Reporting Summary linked to this article.

Source: Ecology - nature.com