Microbial drivers of methane emissions from unrestored industrial salt ponds

Wetlands are important carbon (C) sinks, yet many have been destroyed and converted to other uses over the past few centuries, including industrial salt making. A renewed focus on wetland ecosystem services (e.g., flood control, and habitat) has resulted in numerous restoration efforts whose effect on microbial communities is largely unexplored. We investigated the impact of restoration on microbial community composition, metabolic functional potential, and methane flux by analyzing sediment cores from two unrestored former industrial salt ponds, a restored former industrial salt pond, and a reference wetland. We observed elevated methane emissions from unrestored salt ponds compared to the restored and reference wetlands, which was positively correlated with salinity and sulfate across all samples. 16S rRNA gene amplicon and shotgun metagenomic data revealed that the restored salt pond harbored communities more phylogenetically and functionally similar to the reference wetland than to unrestored ponds. Archaeal methanogenesis genes were positively correlated with methane flux, as were genes encoding enzymes for bacterial methylphosphonate degradation, suggesting methane is generated both from bacterial methylphosphonate degradation and archaeal methanogenesis in these sites. These observations demonstrate that restoration effectively converted industrial salt pond microbial communities back to compositions more similar to reference wetlands and lowered salinities, sulfate concentrations, and methane emissions.

. (a) A bipartite network displaying indicator OTUs grouped at the level of family that are positively and significantly associated (p < 0.001 and r > 0.5) with one or more of the three habitats. The size of each node is proportional to mean OTU abundance (DESeq2 normalized counts) for each type of habitat. Only the 5 most abundant families per site type are shown. (b) The co-occurrence networks of all OTUs were constructed using the DESeq2 normalized counts and conducted Spearman rank correlations between OTUs and visualized with the Fruchterman-Reingold layout using NetworkX. Only strong correlations (Spearman's |r| > 0.9) were visualized in order to reduce complexity. Only indicator OTUs significant from both indicspecies and DESeq2 were labelled in different colors in order to represent their uniqueness to each or two types of sites. The red indicated historic wetlands, the green indicated restored salt ponds and the blue indicated industrial salt ponds. (c) The co-occurrence network of only OTUs that were defined as site-sensitive OTUs of industrial salt ponds was constructed using the matrix of DESeq2 normalized counts and conducted Spearman rank correlations between OTUs and visualized with the Fruchterman-Reingold layout using NetworkX. The node sizes were corresponding to means of DEseq2-normalized counts in unrestored salt ponds. Figure S3. Gene occurrence network constructed using 149 unique CNPS genes and visualized with the Fruchterman-Reingold layout using NetworkX. The significant and positive correlations between these genes with Spearman's r > 0.7 and PFDR < 0.05 were selected for visualization. The nodes representing each gene in the network are labeled with gene names and colored based on their pathway associations, with surrounding halos colored according to subnetwork assignment. The Girvan-Newman algorithm was applied for discovering subnetworks and the optimal assignment was selected based on the highest modularity calculated. Figure S4. The five optimal subnetworks that were selected using the Girvan-Newman algorithm (shown in Figure S5) and visualized separately with the Fruchterman-Reingold layout using NetworkX. The nodes are labeled and colored as in Figure S5. The shadows behind nodes are pie charts indicating relative abundance in historic wetlands and restored and unrestored salt ponds (blue, green and red respectively).    Figure S8. The abundance of phnJ (K06163) and mcrA (K00399) across three types of habitats. The estimated gene copies were normalized using MUSiCC, which corrected these biases of abundance based on universal single-copy genes and machine learning methods. Figure S9. a) Functional genes abundance (DESeq2 normalized counts with z-score normalization) of phosphonate pathway and methylphosphonate synthesis in both sediments and surface water samples. b) Boxplot of mpnS abundance (DESeq2 normalized and log2 transformed) averaged across sample type. Figure S10. The phylogenetic tree constructed based on 56 SCG single copy genes including only 310 bins with methane generation related genes, associated with the presence/absence of these genes, percentage of completeness/contamination, bins' correlation with methane (only bins with Spearman's |r| > 0.5) and MAG abundance (log2 and then z-score transformed). The taxonomy columns (Phylum and Family) were predicted and assigned using the Bin Annotation Tool, BAT (only more than 50% predicted ORFs assigned to the category). Figure S11. The gene occurrence network including 149 unique CNPS genes of 173 bacteria and 3318 archaea from IMG were subjected to phylogenetic profiling analysis. Correlations were calculated using the presence/absence matrix of the 149 unique CNPS genes in genomes. Only the cassettes of genes with Spearman's r > 0.8 and PFDR < 0.05 were considered to have significant co-occurrence within the majority of microbial genomes and selected for visualization. Figure S12. The phylogenetic tree constructed using 56 single copy genes shows that the 310 MAGs are classified into at least 24 phyla. The MAGs from this project are indicated using black dots.

Appendix 1: Functional guild assignments for CH4, N, S, and Fe cycling from 16S rRNA taxa
Microbial functional guilds were assigned based on 16S rRNA taxonomic identities where biogeochemical functions could be ascertained from taxa names. While the authors recognize the inherent limitations of this approach, these functions were assigned where literature reviews indicated functions are relatively monophyletic, or where taxonomic naming conventions are applied consistently. Guild assignments were applied based on taxonomic strings in OTU abundance tables, filtering only organisms matching search terms using grep searches for particular taxa names. Search terms and supporting literature for each microbial guild are described in detail below, and the code to derive these guild abundances from OTU tables is available on github (functions in section 3 of this code).

Methane production and oxidation
Methanogenesis is a monophyletic function carried out exclusively by archaeal taxonomic groups. Abundances of methanogens were obtained by filtering OTU tables according to assigned taxonomy, using the search term "Methano" across taxonomic ranks. This naming prefix appears present in dominant methanogen classes (Methanobacteria, Methanococci, Methanopyri), orders, families, and genera (Hedderich andWhitman 2013, Nazaries et al. 2013). However, some newly discovered methanogen groups would not be recovered by this search term, like the candidate phylum Bathyarchaeota (Bräuer et al. 2020). Within results for methanogenic taxa, a further search was conducted for the order "Methanosarcinales" to approximate the abundances of "acetoclastic" methanogens, while other groups were labeled as "hydrogenotrophic." This classification scheme is an approximation, recognizing that some acetoclastic organisms may also be capable of hydrogenotrophy or methylotrophy (e.g. Methanosarcinacea), while other Methanosarcinales do not use acetate or hydrogen (Hedderich and Whitman, 2013). However, this scheme separates exclusively hydrogenotrophic methanogens from the "acetoclastic" group, while recognizing that methylotrophic methanogenesis does not appear to be monophyletic (Nazaries et al. 2013, Bräuer et al. 2020. Methanotrophic microbial groups were identified from published reviews, and divided into functional guilds based on monophyletic taxonomic groupings (Hedderich andWhitman 2013, Nazaries et al. 2013). Nearly all genera of methanotrophs have names beginning with the string "Methylo" (Knief 2015, Khadka et al. 2018, which was used to filter (grep) taxonomic abundance tables for methanotrophs. A notable exception to this naming convention are Verrucomicrobia organsims, named with the prefix "Methyla," which would not be captured by the initial search term. Guilds of methanotrophic taxa were delineated by further taxonomic subdivision, with search filtering for class-level taxa "Gammaproteobacteria" used to identify Type I methanotrophs, and "Alphaproteobacteria" used to filter Type II methanotrophs, which were further subdivided to separate Type IIa methanotrophs using the search term "Methylocystaceae" on family-level taxonomic assignments (Knief 2015, Khadka et al. 2018. A further class level filter for "Betaproteobacteria" was used to remove methylotrophs from the original "Methylo" search, as this class contains true methylotrophs (e.g. utilizing methanol) like those in the family Methylophilaceae (Chistoserdova 2015, Chistoserdova et al. 2009, Smith and Wrighton 2019. Finally, methane consumption may also be carried out by anaerobic methane-oxidizing archaea, which follow a taxonomic naming convention of "ANME" (Wang et al. 2014) which was used as an additional search filter.

Nitrifiers and anammox
Some microbial nitrogen cycling functions like nitrification and anamox (anaerobic ammonia oxidation) appear to be phylogenetically constrained. Nitrifying bacteria can be classified into functional guilds including ammonia oxidizers (including bacteria and archaea, or AOB and AOA), nitrite oxidizers (NOB), and anammox bacteria (in the Planctomycetes) (Bouskill et al., 2012). Significantly, nitrifying bacteria follow a taxonomic convention where nitrite oxidizers are named with the prefix "Nitro" and ammonia oxidizers with the prefix "Nitroso", except for anammox bacteria. This naming convention is broadly applicable to AOA (Alves et al., 2018), AOB, and NOB (Bouskill et al., 2012), although it may not capture certain Actinobacterial AOB (Khadka et al. 2018). The search term "Nitros" was first used to obtain nitrifying bacteria and archaea, which were then split using an additional search for "Nitroso" to separate ammonia oxidizers. An additional "Archaea" search term at kingdom level taxonomy was applied to separate AOB from AOA. Anammox bacteria were obtained by searching for genus names of the few confirmed genera, including "Kuenenia", "Anammoxoglobus", "Scalindua", "Brocadia", and "Jettenia" (Cai et al. 2020).
Syntrophic bacteria are closely related to sulfate reducers, though some may have lost the capability for sulfate reduction (Plugge et al. 2011). The search term "Syntroph" was used to identify these organisms, given many are found in the order Syntrophobacterales, with three families (Plugge et al. 2011). However, other groups of syntrophic bacteria may exist which would not be obtained by this search (Sieber et al. 2012, Worm et al. 2014.