Terrestrial organic matter (tOM) plays a key role in aquatic ecosystems, influencing carbon cycling and greenhouse gas emissions. Here, we investigate how tOM affects methane production in littoral and pelagic sediments from the Mississippi River headwaters using a controlled microcosm approach. Contrary to expectations, tOM additions consistently enhanced methane production across both sediment types, with no significant differences between littoral and pelagic zones. Methane generation was mediated by specific sediment microorganisms—primarily methanogens and lignocellulose-degrading bacteria—which exhibited similar responses across habitats. Our results highlight the contribution of cytochrome-containing methanogens and their syntrophic interactions with fermentative bacteria, underscoring the importance of microbial community composition in regulating sediment methane dynamics. These findings suggest that increasing tOM loads to freshwater systems could have broader implications for methane emissions, driven by specific microbial interactions.

Table of Contents

• Coverage
• Dataset Description
  • Methods & Sampling
  • Data Processing Description
• Parameters
• Project Information
• Funding

Sample Collection

In June 2020, duplicate gravity cores were collected from five sampling locations within the Mississippi River headwaters (Fig. 1A; Fig. S1 in Sauer et al., 2024: https://doi.org/10.1101/2024.12.31.630949). Of these sites, two were pelagic, two were littoral, and one was riverine. The surrounding watershed lies within Itasca State Park (Minnesota, USA) and is dominated by mature mixed deciduous–coniferous forest. We performed loss-on-ignition (LOI) analyses on homogenized 0–15 cm subsamples from one core per site to determine sediment organic matter content (Table S1 in Sauer et al., 2024). The remaining core from each site was used to establish microcosms representing three treatment groups: high, low, and no terrestrial organic matter (tOM) addition (i.e., leaf litter).

Microcosm Set-up

We used the second core from each of the five sites for microcosms – following a 3x3 design with three tOM (dried leaf litter) treatments (0% tOM addition, 5% tOM addition, 15% tOM addition), each with three replicates. We extruded and homogenized the top 15cm of each core under micro-oxic conditions (using an inflatable anaerobic chamber filled with N2) and added the homogenized mixture to 125mL serum vials. We filled each vial to a depth of 2cm. For each core, we filled the first three vials for the 0% tOM addition treatment and calculated the average sediment weight of the three replicates. We used this value and the starting sediment organic fraction percentage to determine the mass of dried leaf litter needed to increase the organic fraction by 5% and 15%. For tOM additions, we collected and homogenized a dried (105°C for 12h; sieved through No. 5 mesh) mixture of leaves and needles from Acer, Quercus, and Pinus species. We assumed the leaf litter was 100% organic. While there was a statistically significant difference in our tOM addition percentages between treatment (t-test; p<0.001) there was overlap across the percentages and sites, as a result instead of referring to these as firm percentage additions we will refer to them as “low tOM spike” (target 5% tOM) and “high tOM spike” (target 15% tOM). To maintain anoxic conditions, we flushed all the vials with N2 before adding sediment and prior to sealing with a gas-tight butyl-rubber stopper and aluminum crimp.  Finally, we covered the vials in tinfoil and stored them at 4°C in the dark for the duration of the experiment – 180 days. While we recognize that these conditions do not reflect the in situ conditions, our aim was to study methane production dynamics under consistent growth conditions, irrespective of the starting conditions of the initial sediment and abiotic forces. 

Methane Concentrations

We collected 10mL from the headspace of each vial using a gas-tight syringe and injected the sample into Labco Exetainers that were pre-flushed with helium. Exetainers were stored inverted at 4°C until they could be further processed. We replaced the microcosm headspace after every gas pull with 11mL N2 gas to maintain headspace pressure and avoid methanogenesis inhibition caused by the accumulation of CH4 or CO2 (Grasset et al., 2018). We measured CH4 concentrations using a GC-2014 gas chromatograph equipped with a flame ionization detector and an 80/100 Porapak N column (6ft x 1/8in x 2.1mm SS). The carrier gas was argon run at a 25mL min-1 flow rate, and the calibration standards ranged from 200 to 60,000ppm. We then converted these ppm concentrations to molar concentration using the ideal gas law and Henry’s law. All production rates reported in the text were normalized per gram of C in the microcosm and reported as µmol CH4 gC-1 d-1. We used the values obtained from our total organic carbon analysis to normalize production rates by gram carbon in the supplemental materials. To calculate production rates per gram organic C OUT, we assumed the litter was 80% organic matter with a C:N ratio of 20 (Hornbach et al., 2021). Due to product backorders and shipping delays during the pandemic, we did not have enough Exetainers and were unable to pull gas samples from the Mississippi River site for days 4 and 7.

DNA Isolation, Sequencing, and Post-processing 

During the initial microcosm set up, we subsampled ~2g of sediment from the 0-15cm slurry and ~6g of leaf litter and stored these at -20°C until extraction. At the end of the experiment, we combined and homogenized the sediment from the three replicates for each treatment. Again, we collected a ~2g from the pooled and homogenized sediment and stored it at -20°C until extraction. We extracted triplicate DNA for each sample using ~0.25g of sediment or leaf litter using a Qiagen Dneasy PowerSoil Pro kit following manufacturer protocols including 4°C incubation steps. We performed negative controls by carrying out extractions on blanks, using only reagents with no sample. We determined the final bulk DNA concentrations using a Qubit dsDNA HS Assay kit and Qubit Fluorometer. We did not detect any DNA in our blanks (detection limit for the assay kit is 10pg/µL). We then pooled 10µL of each DNA yielding replicate and recalculated bulk DNA for the pooled sample. We sent all DNA yielding samples and blanks to the University of Minnesota Genomic Center (UMGC) for sequencing. 

The UMGC prepared libraries for the samples for Illumina sequencing using a Nextera XT workflow and 2x300bp chemistry. They targeted the V3-V4 hypervariable region of the bacterial and archaeal 16S SSU rRNA gene using primers 341F (5’- CCTAYGGGRBGCASCAG-3’) and 806R (5’- GGACTACNNGGGTATCTAAT-3’) (Yu et al., 2005). The amplicon preparation performed at the UMGC have been shown to be quantitatively more accurate and qualitatively complete than existing methods (Gohl et al., 2016). We recovered a total of 352,403 raw reads from 24 samples, including blanks and leaf litter. We processed these reads using Mothur (v.1.48.0) following the MiSeq SOP (Kozich et al., 2013; Schloss et al., 2009). We aligned our reads using the SILVA database (v.138) and removed chimeras with vsearch (v2.17.1) (Edgar et al., 2011; Quast et al., 2013). Finally, we classified the sequences as operational taxonomic units (OTUs) using a 97% similarity threshold and assigned taxonomy using the SILVA database (Glassman & Martiny, 2018; Stackebrandt & Goebel, 1994). After processing we had 273,151 reads across the 24 samples. 

All further analyses were conducted in R (v4.3.2) using the following packages: tidyverse, phyloseq, vegan, DESeq2, pheatmap, MicrobiomeStat, topicmodels, and ldatuning (Grün & Hornik, 2011; Kolde, 2019; Love et al., 2020; McMurdie & Holmes, 2013; Nikita, 2020; Oksanen et al., 2009; R Core Team, 2018; C. Zhang, 2022). Prior to analyzing the community composition of these sites, we filtered the data by removing any OTU that did not have 2 or more counts in at least 5% of samples. We also removed the nine OTUs which had reads in blank samples, each OTU having a single read. After filtering, the average number of reads per sample was 13,762 and the minimum and maximum read depths were 8,858 and 20,302, respectively. OTU data have a strong positive skew due to many zero counts. To diminish these effects, we used a variance stabilizing transformation (Love et al., 2020). Log-like transformations such as this bring count data to near-normal distributions, produce larger eigengap values, and lead to more consistent correlation estimates all of which influence downstream analysis (Badri et al., n.d.). To compare microbial community structure before and after the microcosm experiment, we conducted a principal component analysis (PCA) of the entire microbial community and the methanogen populations. From both PCAs we pulled the scores of the top two PCs and used those as variables for difference in composition when determining if sediment community composition influences CH4 production. We determined methanogens based on taxonomy and compared the methanogen composition based on energy conservation strategies (i.e., those with or without cytochromes) (Buan, 2018; Ou et al., 2022; Thauer et al., 2008). Finally, we took two approaches to compare the change in populations from initial sediments to post-microcosm. First, we calculated the percent change in each OTU (after adding 10% the lowest observed pre-and-post microcosm abundance to any 0 values). We then selected the top 100 OTUs that had the greatest percent change from initial sediment to post-microcosm ). 

Second, we used Latent Dirichlet Allocation (LDA) or topic modeling to examine the structural differences across microcosm treatments. LDA is a mixture model which correlates microbial communities with relevant environmental factors of interest. The advantages of using LDA over other mixed model approaches is that LDA allows fractional membership allowing samples to be composed of multiple sub-communities. The application of LDA models to microbiome datasets has been described in detail by Sankaran and Holmes 2019 (Sankaran & Holmes, 2019). Briefly, we determined the relevant number of topics or sub-communities (k=34) using the FindTopicsNumber function in ldatuning using a Gibbs sampling method and CaoJuan 2009, Arun 2010 metrics, and Deveaud2014. Then we conducted the LDA model again using Gibbs sampling and the topicsmodels package. We then converted our LDA model back to a phyloseq object to further assess the differential abundance in the 34 sub-communities across select parameters (e.g., treatment addition, treatment quantity, and sample site). For this, we use the linda function in MicrobiomeStat with an alpha of 0.01, winsorization of outliers at 3%, and zero count data handling set to imputation – where in zero counts are given values with respect to sequencing depth (Zhou et al., 2022). From both approaches we aggregated the list of 100 most changed OTUs and those with >1% OUT-sub-community probabilities in the significantly different sub-communities (n=3; 46 OTUs) and evaluated how the abundances of those OTUs explained methane production rates. 

Carbon to Nitrogen Ratios 

We measured total organic carbon (C), total organic nitrogen (N), and isotopic C and N (δ13C, δ15N) on acidified (1M HCl for 24 hr) sediments pre and post microcosm using a Conflow IV open split interface connected to a Delta V Advantage IRMS (Table S2 Sauer et al: https://doi.org/10.1101/2024.12.31.630949). We corrected the isotopic values using a set of known external standards, including sucrose, GA40, and G41. We exhausted the remaining sample for each site and leaf litter and were unable to determine C and N in 2 samples: Pelagic B high tOM spike and leaf litter. 

Methane Concentrations

During the incubation, we sampled the headspace gas for CH4 concentrations at days 1, 4, 7, 19, 60, 90, 150, and 180 using a GC-2014 gas chromatograph. We converted these ppm concentrations to molar concentration using the ideal gas law and Henry’s law, and we normalized production rates per gram C in the microcosm and reported as µmol CH4 gC-1 d-1 (Table S3 in Sauer et al: https://doi.org/10.1101/2024.12.31.630949). 

DNA Isolation, Sequencing, and Post-processing 

We isolated DNA from ~2g of sediment from the initial 0-15cm slurries for each location, ~6g of leaf litter, and ~2g of homogenized sediments from the three replicates for each tOM treatment microcosm post experiment. We sent DNA to the University of Minnesota Genomic Center for sequencing and targeted the V3-V4 hypervariable region of the 16S SSU rRNA using primers 341F (5’- CCTAYGGGRBGCASCAG-3’) and 806R (5’- GGACTACNNGGGTATCTAAT-3’) (Yu et al., 2005). We processed the resulting reads using Mothur (v.1.48.0) following the MiSeq SOP (Kozich et al., 2013; Schloss et al., 2009) and aligned our reads using the SILVA database (v.138) (Quast et al., 2013). We performed all subsequent analyses in R, including Latent Dirichlet Allocation (LDA), which correlates microbial communities with relevant environmental factors. A detailed description of LDA is for microbiomes is provided in Sankaran and Holmes (2019) and also outlined in our SI methods (Sankaran & Holmes, 2019).

 

Coverage: Lake Superior and Lake Erie

NSF Award Abstract:
The Great Lakes hold about 20% of the freshwater on Earth and have been increasingly impacted by human activities in recent decades. Lake Erie suffers from large, annually recurring, toxic cyanobacterial blooms in summer, whereas Lake Superior experiences smaller, localized cyanobacterial blooms after storm events. Cyanobacterial blooms have harmful ecological, human health, and economic implications. These blooms are a global phenomenon, observed in lakes and oceans, and can lead to low oxygen conditions and the production of toxins, both of which can be harmful for ecosystems. Understanding how different types of cyanobacteria influence nutrient cycling remains a major knowledge gap. This project aims to provide a deeper understanding of the long-term state of the Great Lakes ecosystem. The research approach combines new and established methods. Project results and implications will be shared with local and regional water interests in partnership with the Pittsburgh Collaboratory for Water Research, Education, and Outreach, the Great Lakes Commission Harmful Algal Blooms Collaborative, and the Lake Erie Area Research Network. Education is a central part of this project and training opportunities target next generation of scientists, including postdoctoral, graduate, and undergraduate students. The students and postdoc will receive state-of-the-art training in the rapidly developing fields of biogeochemistry and geomicrobiology, while working with an interdisciplinary team of scientists.

This study will examine nitrogen cycling, phytoplankton community composition, and the nitrogen isotopic composition of chloropigments in order to evaluate cyanobacterial productivity in the modern Laurentian Great Lakes as well as the historical record of cyanobacterial blooms over the past several hundred years. The nitrogen isotope composition of chloropigments is expected to provide a powerful new proxy for understanding primary productivity and the relative importance of cyanobacteria to export production and nitrogen cycling. This proxy would be valuable not only for management of modern systems but has important implications for increasing our understanding of the role of cyanobacteria throughout Earth history. This project would test this molecular isotopic proxy in contemporary aquatic ecosystems to assess its efficacy for: (1) determining the relative contributions of cyanobacteria vs eukaryotic algae (e.g., diatoms) to primary production; (2) evaluating export production of cyanobacterial productivity (including blooms); and (3) constraining historical cyanobacteria productivity in the sedimentary record. Comparison of a system characterized by eutrophication and seasonal cyanobacterial blooms (Lake Erie) with one characterized by picocyanobacteria productivity, but the near-absence of large-scale cyanobacterial blooms (Lake Superior), will provide information about the range of impacts that cyanobacteria can have on carbon and nitrogen cycling. Further information regarding nitrogen cycling will be derived from analysis of solid and dissolved nitrogen species throughout the annual cycle, as well as seasonal studies of sediment processes to measure associated sediment nitrogen removal rates through different processes.

This award reflects NSF's statutory mission and has been deemed worthy of support through evaluation using the Foundation's intellectual merit and broader impacts review criteria.