Oxic methane production from methylphosphonate in a large oligotrophic lake: limitation by substrate and organic carbon supply

To evaluate the possibility of OMP occurring in Flathead Lake, we investigated the dynamics of methane and oxygen concentrations in Flathead Lake over a 3-year period (Fig. 1). Sampling was conducted at the long-term monitoring site Midlake Deep (MLD), one of the deepest points in the lake (~113 m). The lake was persistently oxic, with oxygen percent saturation rarely falling below 80%. Methane concentrations typically ranged from 20 to 200 nM and were consistently oversaturated with respect to atmospheric levels. Methane was generally highest and evenly distributed across all depths during the spring. As the summer progressed, a subsurface methane maximum developed. Concentrations were elevated throughout the upper 10–15 m and peaked within the upper thermocline in August and September, suggesting an in situ source of methane within the epilimnion.

We performed metagenomic sequencing to explore the potential for methane generation in Flathead Lake during 2018 (Table S1). The methyl coenzyme M reductase (mcrA) gene, a diagnostic marker for methanogenesis (44), and known methanogenic archaea were absent in our metagenomic and 16S rRNA amplicon data sets. In contrast, we found that up to 15% of the species in Flathead Lake had the capacity to cleave MPn, producing phosphate and methane, using the C-P lyase gene phnJ (Fig. 2). When accounting for the relative abundances of these organisms, they represented between 2% and 10% of the total overall community. The majority of the phnJ diversity was represented by sequences related to the Betaproteobacteria (order Burkholderiales; some databases now classify them as Gammaproteobacteria) and Actinomycetota (previously Actinobacteria or Actinobacteriota; class Ilumatobacteraceae) but also included members of the Alphaproteobacteria, Bacteroidota (previously Bacteroidetes), other Pseudomonadota (Proteobacteria), and Verrucomicrobiota. While the number of organisms with phnJ generally appeared higher within the epilimnion than in deeper waters, neither the percentage of species (Pearson correlation, R < 0.1, P > 0.9) or percentage of organisms (Pearson correlation, R < 0.1, P > 0.75) correlated with methane concentrations. These observations were consistent when looking at the distributions of other genes in the phn operon.

We screened the in situ metagenomes for the capacity for phosphonate production and alternative pathways of methane cycling. Members of the phylum Nitrospirota (previously Nitrospirae) and the family Rhodospirillales (Alphaproteobacteria) possessed the pepM gene for the first step in phosphonate biosynthesis. These genes were present in the deep waters of the lake. We also identified putative, divergent pepM genes related to members of the Dehalococcoideae, Ignavibacteriae, Bacteroidota, Bacillota, Woesarchaeota, Pseudomonadota (<70% similar to Pelagibacter), and Polynucleobacter. These sequences generally showed less than 70% similarity to all currently published sequences. The gene ppd (aepY), involved in the conversion of Pn-pyruvate to Pn-acetylaldehyde, was also divergent in our metagenomes, with genes often showing less than 50% similarity to sequences belonging to the Deltaproteobacteria, Alphaproteobacteria, Woesarchaeota, Nitrospirota, Verrucomicrobiota, Bacillota, and Flavobacteriaceae within the Bacteroidota. We were unable to find the mpnS gene, which encodes the protein for the final step of MPn biosynthesis, in any of the in situ lake metagenomes based on KEGG annotations. While not the focus of this study, we did see evidence that other types of phosphonates may also be important in freshwater systems; genes for the degradation of 2-aminoethylphosphonate (phnW, phnX) were found in abundant members of the Flathead Lake bacterial community, including the genus Limnohabitans, Chloroflexota (previously Chloroflexi), Nitrospirota, and Verrucomicrobiota. Notably, we identified a methane monooxygenase subunit A (pmoA) gene related to the Methylococcaceae, often used to identify methanotrophs, in April at multiple depths when methane concentrations were elevated. This was consistent with 16S rRNA gene amplicon sequencing which revealed that the methanotrophic family Methylococcaceae (<0.2% of all communities) was present at all depths in April but subsequently decreased to below detection the upper 15 m in the summer despite elevated methane concentrations (Fig. S1).

Given the presence of the metabolic potential for OMP via MPn demethylation as revealed by metagenomics, we evaluated the functional capacity for this process using rate measurements of methane production under different nutrient amendments during the summers of 2017, 2018, and 2021 (Fig. 3). These experiments, using water collected from 5 m in July and August, were timed to coincide with periods when methane concentrations were elevated in the upper 10 m of the lake. Amendments included combinations of MPn, phosphate (P), nitrate (N), and glucose (C) under different light regimes (Table 1). In the unamended controls, net rates of methane production were low, averaging 0.55 ± 0.14 nmol CH₄ L⁻¹ d⁻¹ (n = 6 experiments). Treatments amended with only N + C demonstrated net methane consumption (averaging −0.75 ± 0.69 nmol CH₄ L⁻¹ d⁻¹; n = 4 experiments). Similarly, treatments amended with P + N + C yielded net rates of methane production that were lower than in the controls (0.15 ± 0.26 nmol CH₄ L⁻¹ d⁻¹; n = 6 experiments). None of these treatments were significantly different than the control (Dunnett’s test, P > 0.05). In contrast, Flathead Lake water amended with MPn + N + C consistently produced methane after a 2- to 3-day lag time, with rates averaging 116.6 ± 23.2 nmol CH₄ L⁻¹ d⁻¹ (n = 6 experiments; Dunnett’s test, P < 0.05). Methane production was correlated with the drawdown of total dissolved phosphorus (Fig. S2). In some treatments, methane concentrations plateaued at a maximum of ~0.5 µM methane despite remaining MPn. Some of these treatments also showed apparent remineralization of phosphorus over the course of the incubation, with total dissolved phosphorus (TDP) initially decreasing, but then increasing again later in the incubation period. These amendments had the lowest carbon:nitrogen:phosphorus (106:16:1 µM) concentrations added. In contrast, in treatments where carbon (250 µM) or nitrogen (30 µM) was added in excess, methane increased to final concentrations that approached the total added MPn. When samples were amended with both MPn + P, along with labile carbon and nitrogen sources, methane production rates were significantly repressed, averaging only 13.0 ± 10.1 nmol CH₄ L⁻¹ d⁻¹ (n = 2 experiments). No differences were observed between amendments incubated in the dark, in the light, or in a dockside incubator maintaining in situ light and temperature fluctuations (dark, 0.58 ± 0.20 nmol CH₄ L⁻¹ d⁻¹; constant or in situ light, 0.49 ± 0.23 nmol CH₄ L⁻¹ d⁻¹; t-test, P > 0.05).

To evaluate the community response to MPn, we measured a more exhaustive list of properties during an amendment using lake water from 5-m depth in Summer 2020 performed under simulated in situ lake light and temperature conditions. This included measurements of cell abundances, chlorophyll a (chl a), nutrient concentrations, and metagenomic sequencing. Methane concentrations were highest in the upper 15 m at the time of sampling for this experiment, with in situ concentrations in excess of 1,000% of atmospheric equilibrium (Fig. S3). The rate of methane production without amendment was 1.51 nmol CH₄ L⁻¹ d⁻¹ (Fig. 4). Methane production was not stimulated through the addition of P + N or P + N + C (1.11 and 1.32 nmol CH₄ L⁻¹ d⁻¹, respectively). The addition of MPn + N also did not substantially increase methane production (2.46 nmol CH₄ L⁻¹ d⁻¹). In contrast, as previously observed, the addition of MPn + N + C increased rates of net methane production to >50 nmol CH₄ L⁻¹ d⁻¹, with methane produced stoichiometrically with the amount of phosphorus consumed. Consistent with the previous experiments, the addition of P alongside MPn + N + C inhibited methane production. Greater than 75% of the amended P was consumed in the MPn + N + C, P + N + C, and P + N treatments where 1 µM P was added (Fig. 4; Fig. S4). In contrast, less than ~10% of the available phosphorus and nitrate were consumed in the MPn + N treatment. TDP also remained in treatments where 2 µM phosphorus was added, likely reflecting carbon or nitrogen limitation.

We performed metagenomic sequencing to identify the microorganisms which responded to the Summer 2020 amendment which simulated in situ lake light and temperature conditions. Based on 16S rRNA gene composition, at the onset of the experiment the in situ lake community contained high abundances of members of the Actinomycetota (Fig. S5), including sequences similar (>99%) to those from the genera Nanopelagicus, Planktophila, and other members of the hgcI clade. The most abundant 16S rRNA gene was related to the clade SAR11, identical to that obtained from Lake Superior (45) and 99.8% similar to the isolate Fonsibacter ubiquis (46). Abundant sequences related to the Burkholderiales included those similar to Limnohabitans, Limnobacter, Polynucleobacter, Methylopumilus, and Rhodoferax. Other taxa of note included those related to the Methylacidiphilaceae within the phylum Verrucomicrobiota, the cyanobacterial genus Cyanobium, and the Flavobacterium within the Bacteroidota. These organisms are consistent with those reported from other freshwater lakes (47–50), indicating that our starting community reflected the diversity typically found in freshwater systems. At the end of the experiment, the community in the unamended control was similar to the starting sample.

Nutrient amendments strongly influenced community composition. Photosynthetic and eukaryotic taxa responded to the P + N treatment based on cell abundances, chl a concentrations, and predicted eukaryotic contigs (Fig. S6). Cyanobacteria belonging to the genus Cyanobium were present, along with Bacteroidota which are known to be associated with phytoplankton blooms [e.g., references (51–53)]. Eukaryotic 18S rRNA genes also became more prevalent in the P + N treatment, including rRNA genes deriving from diatoms Bacillariophyta (Thalassiosira, Skeletonema), the green algae Chlorophyta, and the ciliate Ciliophora. In contrast, the microbial community from the treatment amended with MPn + N remained largely similar to the control community, reflected by high abundances of Fonsibacter, the actinobacterial hgcI clade, and diverse members of the Burkholderiales. In treatments amended with P + N + C, we observed a strong response of members of the Bacteroidota. The most abundant sequences in this treatment were related to the Spirosomaceae and included the genus Flectobacillus (54), an organism which can change its cell size to avoid grazing (55–58). In addition, members of the Alphaproteobacteria, including sequences related to the genus Caulobacter and the groups Sphingomonadales and Rhodobacteraceae, and the phylum Actinomycetota (genera Rhodoluna and Aquiluna) all responded to the P + N + C treatment. Uncharacterized sequences within the Burkholderiales, which were distantly related (<96% similar) to Rhodoferax and Acidovorax, were also present.

In contrast to the above community responses, samples amended with MPn + N + C were dominated by members of the Burkholderiales. The most abundant 16S rRNA gene sequences belonged to the genus Acidovorax and were similar to Acidovorax sp. RAC01 (59). Sequences related to Rhodoferax, Curvibacter, Ideonella, and Hydrogenophaga were also present (60, 61). Altogether, members of the Betaproteobacteria represented more than 90% of the community following amendment with MPn + N + C. The alphaproteobacterial family Rhizobiaceae also made up 1% of these communities. Despite the presence of MPn, N, and C, the MPn + P + N + C amendment was comparable to samples amended with P + N + C, showing an enrichment of Bacteroidota, Microbacteriaceae, and Sphingomonadales.

To identify organisms which responded to our treatments that had the biosynthetic capacity for methylphosphonate use, we screened the amendments for the phnJ gene (Fig. 4). Approximately 7% of the species in the control had the phnJ gene, consistent with our analyses from 2018 (Fig. 2). These sequences were similar to those reported above, including members of the Betaproteobacteria (Burkholderiales; Rhodoferax, Curvibacter, Limnohabitans, and Methylopumilus), Actinomycetota (Actinomycetia and Ilumatobacteraceae), Alphaproteobacteria (Rhodospirilalles, Caulobacteraceae, and Sphingomonadales), and Bacteroidota. In contrast, over 40% of the species in the MPn + N + C amendments had phnJ, showing that MPn amendment in the presence of labile carbon and nitrogen selected for microorganisms capable of using MPn. These sequences were almost entirely related to the Burkholderiales and included the genera Acidovorax and Rhodoferax. While the in situ lake samples also had sequences distantly related to Acidovorax and Rhodoferax, the MPn + N + C enrichment phnJ were unique (sharing <90% amino acid identity), indicating distinct microorganisms in these samples.

To further explore the metabolic potential within each amendment, we performed genome binning (Table 2). We report metagenome-assembled genomes (MAGs) from a wide group of phyla that reflect the most abundant organisms in each treatment (Fig. 5; Fig. S7). This includes two MAGs belonging to the genus Acidovorax which were the most abundant genomes in the MPn + N + C and 2MPn + N + C amendments. These genomes shared 98.3% average nucleotide identity, higher than the 95% average nucleotide identity (ANI) demarcation used for identifying distinct species (62, 63). To see if the response of this organism was consistent across not only replicates but also across years, we sequenced one 2018 MPn + N + C enrichment. We discovered that the 2018 enrichment also contained an abundant Acidovorax MAG which shared >98.5% ANI to those from 2020. All three MAGs shared >98% ANI with strain UKL202B from Klamath Lake (64). A further Acidovorax genome was obtained from the MPn + N enrichment; however, this MAG shared less than 83% ANI with the other three Acidovorax genomes. Similarly, related yet distinct MAGs belonging to the genus Rhodoferax showed differential read recruitment between amendments (Fig. S7). Altogether, these findings suggest the presence of many closely related members of the Burkholderiales in Flathead Lake that may respond in distinct but reproducible ways to the availability of nutrients and organic matter.

Given the consistent response of members of Acidovorax to MPn enrichment, we explored the distribution of genes involved in phosphorus use within representative members of this genus (Fig. 5). Genes involved in phosphate transport (pstCABS) were present in every genome. The phn operon (phnFGHIJKLMNP) was widely distributed, indicating that phosphonate use is not specific to species in Flathead Lake. Genomes of various members of the Acidovorax that do not possess the capacity for MPn degradation appear associated with hosts and generally fall along the lines of a recently proposed genus division (65, 66). While rare, these organisms instead have the capacity for phosphonate production, including pepM and ppD/aepY to convert phosphoenolpyruvate to Pn-acetylaldehyde and its subsequent degradation to either acetate or acetaldehyde using phnYA or phnX. None of the genomes had mpnS, the gene responsible for the final step in the production of MPn. Ultimately, we conclude that nearly every organism in the genus Acidovorax has some capacity for phosphonate cycling, although mechanisms of phosphonate production and MPn demethylation appear present in discrete species.

To evaluate the distributions of phnJ-containing Acidovorax, Rhodoferax, and Allorhizobium in Flathead Lake, we used metagenomic read recruitment and 16S rRNA amplicon sequencing. Metagenomic read recruitment against these MAGs showed that these organisms were rare in the in situ lake community, likely representing less than 1% of the total population (Fig. S8). Amplicon sequencing revealed one amplified sequence variant (ASV) identical to the dominant Acidovorax 16S rRNA gene sequence from our 2020 MPn + N + C amendment and shared 99.71% similarity with the 16S rRNA gene sequence from the 2018 MPn + N + C amendment. This ASV was a rare member of the Flathead Lake community, found only in the epilimnion during the summer and reached a maximum relative abundance of 0.30% (Fig. 5).

Methane concentrations in Flathead Lake were measured approximately monthly over a 3-year period in 2018–2021. Water was placed in glass serum bottles crimp-sealed with gray chlorobutyl rubber stoppers and injection of NaOH to a final concentration of 0.1 M was used to stop biological activity. A headspace was introduced into each sample using ultra-high pure N₂ gas followed by agitation to allow gas concentrations in the headspace to equilibrate with the liquid sample. Methane was analyzed by gas chromatography via headspace (20 mL) gas introduction (Model 8610C, SRI Instruments). Methane and oxygen concentrations were visualized and contoured using Ocean Data View (86).

To document the microbial diversity and potential for phosphonate use in Flathead Lake, we performed metagenomic sequencing from 16 samples collected at MLD in 2018. Depths ranged from 5 to 90 m. Approximately 1 L was serially filtered onto both a 3-µm, 25-mm GTTP polycarbonate filter (EMD Millipore, MA, USA) and a 25-mm, 0.2-µm polyethersulfone filter (SUPOR, Pall Co., NY, USA). Samples were stored at −80°C prior to DNA extraction. Genomic DNA from the <3.0 to >0.2 µm fraction was extracted using a MasterPure DNA purification kit (Lucigen, WI, USA). DNA libraries were sequenced on a Novaseq (NEB Ultra II DNA library prep kit; Novogene, Sacramento, CA), a NextSeq 2000 (Illumina DNA Prep kit; Microbial Genome Sequencing Center, Pittsburgh, PA), or a MiSeq (Nextera; Genomics Core Facility, University of Montana).

Raw reads were quality trimmed using Trimmomatic v0.39 (87) and assembled with MEGAHIT (88) using default parameters. The depth of coverage of contigs was estimated using Bowtie 2 v2.3.5.1 (89) and SAMtools v1.10 (90). Open reading frames were identified with Prodigal V2.6.3 (91). Gene annotation was performed using GhostKOALA (92). C-P lyase genes (phnJ; K06163) involved in MPn cleavage were classified taxonomically against the nr database using blastP (cut off 75%) (93) and through the construction of phylogenetic trees using reference sequences. Trees were created by alignment with MUSCLE (94), built using FastTree (95), and visualized using iTOL (96). We classified sequences at the phylum or class level. We estimated the percentage of species able to cleave MPn by dividing the number of phnJ genes by the total number of organisms estimated using four single-copy marker genes (recA, K03553; gyrB, K02470; atpD, K02112; tufA, K02358) as previously performed (23, 30, 68). We also estimated the relative abundance of the total community represented by organisms with phnJ by using the depth of coverage of all phnJ and single-copy marker genes within each metagenome. We further identified marker genes for phosphonate production (pepM, K01841 and K23999; ppd, K09459; mpnS, K18049), the degradation of 2-aminoethylphosphonate (phnW, K03430; phnX, K05306), anaerobic methanogenesis (mcrA; K00399), and methanotrophy (pmoA, K10944; mmoX, K16157).

To determine methane production rates under different nutrient concentrations, amendments were performed during the summers of 2017, 2018, 2020, and 2021. For all years, water was collected from MLD at a depth of 5 m, placed in 20-L polycarbonate carboys, and returned to the Flathead Lake Biological Station. Carboys were amended with MPn, P, N, and C at different final concentrations and molar ratios and incubated under specific light and temperature conditions depending on the experiment (Table 1). Carboys were subsampled into sealed borosilicate glass serum vials for subsequent incubations. Amendments from 2017 and 2018 were maintained in a laboratory incubator at ~22.5°C. While most of these amendments were maintained in the dark, one experiment (24 July 2018) was performed under both constant dark and constant light. The incubator provided photosynthetically active radiation (400–700 nm) of ~400 µmol photon m⁻² s⁻¹. To mimic in situ light and temperature conditions more closely, amendments in 2020 and 2021 were maintained in a dockside incubator. The incubator was plumbed to continuously supply lake water from ~2-m depth with minimal shading (at midday; ~16°C, ~70% incident photosynthetically active radiation, ~900 µmol photon m⁻² s⁻¹). A blue plastic layer made of plexiglass served to approximate the underwater spectral quality (97, 98). Amendments from 2017, 2018, and 2021 were performed using triplicate 160-mL glass serum vials with gray chlorobutyl rubber septa. All July 2020 amendments were performed with one replicate using individual 1-L glass serum bottles sealed with a blue butyl rubber septum.

Methane concentrations were followed over time to calculate rates of production. All experiments except July 2020 sacrificed whole samples at each time point. In the case of the 2020 experiment, 15 mL of sample was withdrawn and replaced with 15 mL of zero air lacking hydrocarbons. Samples were preserved with a final concentration of 0.1 M NaOH. Methane concentrations were determined as described above. Rates of methane production were calculated between equilibrated T₀ methane concentrations and methane concentrations following 5–7 days of incubation. Because methane production showed a lag period following nutrient amendment, these rates do not necessarily reflect the maximum rates between any two time points but allow for comparison across all experiments.

For determination of nitrate + nitrite (NO_x) and TDP, samples were filtered through MilliQ and lake water rinsed, 47-mm diameter, 0.45-µm pore size mixed cellulose ester filters, and frozen at −20°C until analysis on an Astoria A2 segmented flow analyzer (Astoria-Pacific, OR, USA). TDP was measured colorimetrically following wet chemical oxidation using an alkaline potassium persulfate digestion and treatment with heteropoly-molybdenum blue. For NO_x determinations, nitrate was converted to nitrite via cadmium reduction and quantified colorimetrically using Greiss chemistry.

A more comprehensive set of measurements were performed during the amendment on 20 July 2020. These measurements are described below.

Total cell abundances and chl a-containing phototrophic organisms were determined using an Attune Acoustic Focusing flow cytometer (Thermo Fisher, MA, USA). Samples (2 mL) were fixed with paraformaldehyde (0.8% final concentration) and frozen at −80°C until analysis. For total cell abundances, cells were stained with the nucleic acid stain SYBR Green I for 15 minutes, excited using a 20-mW blue (488 nm) excitation laser, and detected with a 530-nm emission filter. Phototrophic cells containing chl a and phycoerythrin were excited using a blue (488 nm) excitation laser at 20 mW. Small phototrophic eukaryotes were detected based on chl a emission using a 640-nm longpass filter, while cyanobacteria containing phycoerythrin were distinguished using an R-phycoerythrin emission filter (574 nm). Concentrations of chl a were determined following filtration of lake water (100 mL) onto a 25-mm diameter 0.45-µm pore size mixed cellulose ester filter (Millipore). Filters were extracted at −20°C in a 90% acetone solution overnight and the fluorescence quantified using a Turner 10-AU fluorometer (Turner Designs, CA, USA), including phaeophytin correction (99).

We performed metagenomic sequencing to determine community shifts in response to the experimental amendments. Approximately 100–500 mL of water was filtered onto a 25-mm diameter, 0.2-µm polyethersulfone filter (SUPOR, Pall Co., NY, USA) and stored at −80°C. Genomic DNA was extracted using a MasterPure DNA purification kit (Lucigen, WI, USA). DNA libraries were prepared using an Illumina DNA Prep kit (Illumina, San Diego, CA) and 150-bp paired-end reads were sequenced on a NextSeq 2000 at the Microbial Genome Sequencing Center (MiGS; Pittsburgh, PA). To see if the community response to MPn amendment was consistent across years, we performed further metagenomic sequencing on an MPn + N + C amendment from 2018 (July 24, dark incubation). The amendment was serially filtered through a 3-µm, 25-mm GTTP polycarbonate filter (EMD Millipore, MA, USA) onto a 25-mm, 0.2-µm polyethersulfone filter (SUPOR, Pall Co., NY, USA). Sample storage and DNA extraction were performed on the 0.2-µm filter as described above. Library preparation and shotgun 150-bp paired-end sequencing (150 bp) were performed on a Novaseq (Novogene, Sacramento, CA).

Read clean up, assembly, and functional annotation were performed as described for the in situ communities. To evaluate microbial community composition, 16S ribosomal RNA genes were identified using barrnap (100) and classified using the SINA Aligner (101) against the SILVA database (102). Only genes identified as prokaryotic were retained. The depth of coverage of each 16S rRNA gene was used to estimate the percent abundance of each organism within the community. We further compared metagenomes against one another based on functional potential. The coverage of genes with the same KEGG annotation was summed within each metagenome. KEGG abundances within each metagenome were rarefied to equal sampling depth and compared using NMDS ordinations based on Bray-Curtis dissimilarity. Permutational analysis of variance, based on the type of phosphorus added (MPn or phosphate), was performed using adonis in the package vegan (103). To explore the response of eukaryotic organisms, contigs predicted to belong to eukaryotes were identified using EukRep (104). The percentage of sequencing depth attributable to eukaryotic organisms in each metagenome was estimated based on the total reads that mapped to these contigs relative to the total number of reads.

Metagenome-assembled genomes were obtained using MetaBAT 2 v2.11.1 (105). Contigs >5 kb in length were retained. The size and quality of each genome bin were evaluated using QUAST v5.0.2 (106) and CheckM v1.0.13 (107). We report genome bins more than 60% complete and less than 10% contaminated, representing medium quality draft genomes (108). Genomes were taxonomically classified with GTDB-tk v1.6.0 (109) using KBase (110). Genome annotation was performed using Prokka v1.14.6 (111) and GhostKoala. Whole-genome trees of members of the genus Acidovorax were created using concatenated single-copy marker genes identified using CheckM, constructed using FastTree, and visualized using iTOL. Outgroups for Fig. 5B were members of the genus Hydrogenophaga. ANI comparisons were performed using OrthoANI (112). The relative representation of each bin within each metagenome was estimated by the number of reads mapped per kilobase million reads.

To further evaluate the distribution of Acidovorax in Flathead Lake, we performed 16S rRNA gene amplicon sequencing at MLD in 2018. Water collection, DNA extraction, and sequencing were performed as previously described (113). Sequence data were processed using the QIIME2 platform (114). ASVs were identified using DADA2 (115) and classified against the SILVA 138 database (116). Further processing and visualization were performed using phyloseq (117) in R (118). Sequences related to chloroplasts and eukaryotes were removed. We used blastn to identify ASVs identical to the most abundant 16S rRNA gene related to Acidovorax in the 2020 MPn + N + C amendment. We also used this data set to identify potential methanogenic archaea and methanotrophic organisms, including the class Methylococcaceae.