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). 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.