Lakes were located in northern Wisconsin, United States. We sampled 13 lakes between June and August 2019 and an additional 9 lakes between July and August 2020 (Table 1). During each sampling season, the sampled lakes represented the full range of the pH gradient. To construct sediment slurry incubations, we sampled sediment and hypolimnetic water from the deepest point of every lake. Sediment was collected from the top 15 cm of the sediment surface with an Eckman dredge and hypolimnetic water was collected from 0.25 m above the sediment surface with a Van Dorn water sampler. There was some degree of atmospheric oxygen exposure during sampling for both hypolimnion water samples and sediment samples, however we attempted to mitigate the length of exposure. Hypolimnion water samples were sealed immediately after collection the field. Sediment samples were collected in 5-gallon open-air buckets. Upon return to the laboratory, sediment was sampled from the bottom of the bucket for construction of sediment incubations and construction occurred within 5 h of collection. Prior to construction, all samples were stored in the dark at 4°C.

We constructed experimental sediment incubations in 300 ml glass serum bottles containing 50 ml of sediment and 50 ml of anoxic hypolimnion water. To quantify the difference in CH₄ production under increased OM supply, we conducted triplicate incubations with and without additions of OM. For incubations treated with the addition, we added 20 mg of dried algal biomass (Scenedesmus obliquus). S. obliquus was grown in Bold’s Basal Medium under constant aeration with atmospheric carbon dioxide concentrations. To generate a homogeneous dried material for additions, we centrifuged the wet material in 40 ml batches at 4,000 rpm for 10 minutes, extracted the pellet, and dried the pellet for 7 days. All dried material was added to a single vial and homogenized before addition to incubations. After additions, we capped and sealed the serum bottles with rubber septa and aluminum crimp seals and homogenized the sediment slurries by manual shaking for 1 min. We then purged the headspace of each incubation with N₂ gas for 5 min to maintain anoxic conditions. We incubated the bottles in the dark at 4°C to simulate conditions at the bottom of a lake and we incubated the bottles for approximately 28 days.

During the length of the incubation period, we estimated net CH₄ production rates from each incubation by iteratively sampling from the headspace for CH₄ concentration. The first sampling occurred 24 h after purging the headspace and sealing the incubations. We then sampled every 3–5 days, as well as at end of incubation (approximately 28 days), for a minimum of 4 sampling points for each incubation (Supplementary Figure S1). To sample for CH₄ concentration, we extracted 10 ml of headspace gas using a 25-gauge needle pierced through the rubber septa. Samples were immediately injected into airtight glass vials for measurement of CH₄ concentration on an Agilent 6890 Gas Chromatograph with parameters previously described in West et al. (2016). We then added 10 ml of N₂ back to the headspace to maintain pressure. We inferred net rates of CH₄ production, after accounting for the headspace dilution, from the slope of a linear regression fit to the time-series (Supplementary Figure S1). For incubations in which CH₄ did not accumulate over time (N = 4), we assumed a net CH₄ production of zero.

During the initial collection of lake sediments for construction of incubations, we sampled approximately 1 g of sediment for analysis of microbial community composition and abundance. Sediments were stored in sterile plastic 1.5 ml cryovial tubes at −80°C until DNA extraction. We conducted a single DNA extraction for each lake using 0.25 g of sediment and the DNeasy PowerSoil Pro kit (Qiagen, Hilden, DE) according to manufacturer instructions. The extracted DNA was used as template for 16S rRNA paired-end gene sequencing to determine differences in community composition. Extracted DNA was also used as template in digital droplet polymerase chain reactions (ddPCR) of the 16S rRNA gene and the alpha subunit of methyl coenzyme reductase (mcrA) gene to determine total bacterial + archaeal abundance and methanogen abundance, respectively.

For 16S rRNA paired-end gene sequencing, the V4 region of the 16S rRNA gene was amplified in a 15 μl PCR reaction with dual indexed primers: 515F (5′-GTGCCAGCMGCCGCGGTAA-3′) and 806R (5′-GGACTACHVGGGTWTCTAAT- 3′) (Kozich et al., 2013). Amplified DNA served as template for high-throughput paired-end (2 × 250 bp) sequencing on an Illumina MiSeq v2 Standard flow cell at the Michigan State University Research Technology Support Facility. Sequencing resulted in 1,279,135 high-quality reads across all lake samples. Raw sequence reads are available at the NCBI Sequence Read Archive (Accession # PRJNA813558). To determine differences in microbial community composition, we used the mothur (version 1.44.3) bioinformatics pipeline for merging of paired-end reads, quality filtering, clustering, and taxonomic classification (Kozich et al., 2013). To determine community membership, operational taxonomic units (OTUs) were defined at 97% similarity using the dist.seqs command and a 0.03 cutoff. Representative sequences for each OTU were aligned and classified against the RDP classifier database (version 18) (Wang et al., 2007).

To quantify abundance of the 16S rRNA gene and mcrA gene, we conducted ddPCR using the BioRad QX200 Droplet Digital PCR system. The indicator genes were amplified in separate 20 μl ddPCR reactions using a BioRad C1000 Touch thermocycler with EvaGreen as the reporting dye. For amplification of the 16S rRNA gene, each reaction contained 2 μl of 1/100 diluted sediment DNA template, 1X Q200 ddPCR EvaGreen Supermix, and 0.1 μM of each primer targeting 16S rRNA: 338F (5′-ACT​CCT​ACG​GGA​GGC​AG-3′) and 805R (5′-GAC​TAC​CAG​GGT​ATC​TAA​TC- 3′) (Yu et al., 2005). For amplification of the mcrA gene, each reaction contained 1 μL of 1/10 diluted sediment DNA template, 1X Q200 ddPCR EvaGreen Supermix, and 0.25 μM of each primer targeting mcrA: mcrAqF (5′-AYGGTATGGARCAGTACGA- 3′) and mcrAqR (5′-TGVAGRTCGTABCCGWAGAA-3′) (West et al., 2012). Thermocycling conditions for both indicator genes are reported in Supplementary Table S1. Additionally, for both assays, we performed three technical replicates on a subset of the samples to verify precision. Technical replicates indicated high precision of the assays (Supplementary Figure S2), and we normalized copy numbers to copies per gram of sediment to be comparable across samples.

During sampling for construction of incubations, we measured dissolved oxygen concentration and pH at the sediment-water interface using a YSI Professional Plus Multiparameter meter (Yellow Springs Instruments, Yellow Springs, OH, United States). We also collected additional sediment to quantify sediment OM content and epilimnion water for quantification of water column chlorophyll a (chl a) concentrations, as metrics of differences in OM availability and lake trophic status. We quantified percent OM of sediment using loss on ignition measurements of dried sediment samples (Heiri et al., 2001). To quantify water column chl a, we used particles from 450 ml of the epilimnion sample captured on a 0.7 μm glass fiber filter and we analyzed the sample using methanol extraction and fluorometry (Welschmeyer, 1994).

We conducted all statistical analyses in R (R Core Team, 2018) using the base and vegan (Oksanen et al., 2017) packages. To quantify variation in sediment microbial community composition across the 22 lakes, we first rarified every sample to 39,948 sequence reads to be comparable across sample sequencing depths. We then identified methanogen OTUs by the presence of “Methano” in taxonomic assignments, as a conservative means of identifying methanogens, and manually curated our list of putative methanogen genera for known methanogen taxa. Some methanogenic microbes were likely not identified based on this method (Kharitonov et al., 2021), but our approach did allow for estimation of relative sequence recovery of known methanogenic groups. We then conducted the following analyses for both the methanogen and non-methanogen community.

To quantify variation in microbial community composition across the 22 lakes, we standardized OTU matrices of both communities to reflect relative sequence recovery and generated pairwise Bray-Curtis dissimilarity matrices of each community using the vegdist function. We then generated principal coordinates for each community based on the Bray-Curtis dissimilarity matrices using the cmdscale function. The first principal coordinate axes (PCoA 1) represent 40.6 % and 29.2% of the variation in the methanogen and non-methanogen community, respectively. We also tested for effects of environmental conditions on community composition using permutational multivariate ANOVA (PERMANOVA) tests using the adonis function in R (Anderson, 2001). We evaluated significance of the PERMANOVA statistic using 1000 permutations of the dissimilarity matrices for each community. To quantify differences in microbial community diversity and richness, we also calculated the Shannon Diversity Index and community richness, or the sum of OTUs present in each sample. Thus, we used the Shannon Index, community richness, and the PCoA 1 scores for both communities as metrics of community composition in subsequent analyses.

To quantify the response of sediment CH₄ production to the OM additions across the 22 lakes, we calculated the difference in mean net CH₄ production between incubations with and without OM additions (∆CH₄ production) for every lake. We also quantified the maximum percent utilization of OM additions by calculating the total CH₄ produced through the length of the experiment divided by the amount of OM supplied. For the dried autochthonous OM supplied, we assumed a carbon content of 61.5% based on measurements from West et al. (2015). Finally, we determined significant correlations between the log-transformed ∆CH₄ production and metrics of microbial community composition and abundance using simple linear regressions. We used the log-transformation to abide by assumptions of normally distributed data. We also determined the correlation between log-transformed ∆CH₄ production and metrics of preexisting OM supply (sediment OM %, and water column chl a concentration). Data, metadata, and code for all analyses are available on GitHub at https://github.com/brittnibertolet/CH4response2OM.

Sampling along a gradient of pH resulted in sediment incubations with variable microbial communities. After rarifying, sequencing recovered 367 methanogen OTUs and 26,494 non-methanogen OTUs. Methanogen reads represented 1.5% of the total sequence reads and the relative sequencing reads of methanogens ranged from 0.7 % to 2.8% across all lakes. As hypothesized from previous work (Bertolet et al., 2019), the composition of both the methanogen and non-methanogen community was significantly related to pH and sediment OM (Figure 1). As such, PCoA Axis 1 of both communities was significantly correlated with pH (methanogen: R ² = 0.73, p < 0.001, non-methanogen: R ² = 0.69, p < 0.001) and OM (methanogen: R ² = 0.71, p < 0.001, non-methanogen R ² = 0.76, p < 0.001). Additionally, the relative abundances of different methanogen genera (Supplementary Figure S2) and non-methanogen phyla (Supplementary Figure S3) varied across lakes.

Microbial richness and diversity varied significantly across the lake samples as well as between the methanogen and non-methanogen communities. Methanogen richness ranged from 13 to 60 OTUs, while non-methanogen richness ranged from 1,404 to 6,059 OTUs. Methanogen Shannon Index ranged from 1.7 to 2.8, while non-methanogen ranged from 4.9 to 7.3. Additionally, neither methanogen richness nor diversity were significantly correlated with environmental variables (pH, sediment OM %, or chl a concentration). In contrast, non-methanogen richness and Shannon Index were both significantly correlated with both pH and sediment OM % (Supplementary Table S2).

Digital droplet PCR of the 16S rRNA and mcrA genes also indicated variation in the total bacterial + archaeal abundance and methanogen abundance across lakes (Supplementary Figure S4). Total bacterial + archaeal abundance ranged from 3.07 × 10⁵ to 1.22 × 10⁸ copies per gram wet sediment. However, total bacterial + archaeal abundance was not related to sediment pH nor sediment OM %. Instead, total bacterial + archaeal abundance was significantly correlated with water column chl a concentration (R ² = 0.24, p < 0.05; Supplementary Figure S5). Similarly, methanogen abundance ranged from 3.31 × 10⁵ to 1.79 × 10⁷ copies per gram wet sediment and was significantly correlated with 16S rRNA copy number (R ² = 0.26, p < 0.05) and water column chl a concentration (R ² = 0.19, p < 0.05; Supplementary Figure S4). However, the relative sequence reads of methanogens, derived from the 16S rRNA Illumina sequencing, was not related to ddPCR estimates of microbial abundance nor chl a concentration. Additionally, neither differences in community composition nor abundance were related to sample year (Supplementary Table S4).

Rates of CH₄ production were variable across lakes and were significantly influenced by autochthonous OM additions. OM additions increased the net CH₄ production rate across all lakes (Kruskal-Wallis test: χ ² = 73.2, p < 0.001; Figure 2). However, the magnitude of the response (∆CH₄ production) varied and ranged from 1.1 to 27.3 μmol CH₄ L sediment slurry⁻¹ day⁻¹. Additionally, maximum potential utilization of OM ranged from 0.28 % to 7.1% of algal carbon added to a slurry.

Variation in the ∆CH₄ production was related to differences in both metrics of sediment microbial community composition and pH across the sampled lake sediments (Figure 3). Specifically, the Shannon Diversity Index and richness of the non-methanogen explained the most amount of variation in ∆CH₄ production, and pH and metrics of methanogen community composition were also significantly correlated with the log-transformed ∆CH₄ production (Table 2). However, neither abundance metrics nor environmental variables related to OM supply (sediment OM % or water column chl a) explained variation in ∆CH₄ production, despite significant influences of sediment OM % on both methanogen and non-methanogen community composition (Figure 1). Additionally, the sample year did not influence patterns ∆CH₄ production (Supplementary Table S2).