Introduction

Salt marshes are globally important to carbon storage and have the potential to influence both the changing climate and the stability of coastal systems. Salt marsh sediments contain between 400 – 6500 Tg of organic carbon (Duarte et al. 2013; Mcleod et al. 2011) and store an additional 10.2 - 44.6 Mt annually, the equivalent of 0.5-1.0% of anthropogenic C emissions (Ouyang and Lee 2014). However, carbon storage in salt marshes is declining and estimates indicate that over a recent nine-year period, salt marsh habitat loss resulted in emission of 16.3 Tg of CO₂ and a reduction of 0.045 Tg of CO₂ burial per year (Campbell et al. 2022).

Rapidly increasing sea level rise and warming climate could destabilize current stocks of organic carbon with unpredictable consequences (Fontaine et al. 2007; Spivak et al. 2019). Microbial communities within salt marsh sediments are compositionally and metabolically diverse and their interactions influence pools of soil organic carbon (SOC) (Bulseco et al. 2019; Erb et al. 2024). However, there are important knowledge gaps in our understanding of the metabolic capacity and distribution of microbial communities in 2this system often lead to the identification of previously unidentified carbon cycling among novel taxa (Payne et al. 2025; Seitz et al. 2016; Vineis et al. 2023). Identifying how the pools of carbon are influenced by microbial communities will enhance our understanding of the fate of blue carbon stocks, how they are formed, and how they change over time.

The decomposition of organic material is strongly influenced by plant communities that blanket the surface of the marsh in nearly homogenous stands. The regions of the marsh inhabited by plants are highly dynamic zones that contain hot spots of microbial activity due to plant-microbial interactions and feedbacks that are sensitive to seasonal and tidal influence (Kuzyakov and Blagodatskaya 2015). Spartina patens is a highly productive salt marsh macrophyte with roots that release organic carbon and aerate the sediment, thereby influencing carbon sequestration of salt marshes through the “microbial carbon pump” (Liang et al. 2017; Spivak et al. 2019). Live roots of S. patens are dominantly 0-10cm below the sediment surface (Gallagher and Plumley 1979) and plant inputs could influence the SOC pool up to 100 cm (Liu et al. 2017).

The patchy distribution of resources and dynamic processes in marsh sediments present challenges for unraveling the underlying rules for community assembly and metabolism. Identifying the distribution and putative functions of the microbial genomes within the dynamic zone influenced by plants and tidal inputs remains one of the most challenging but critical aspects to understanding the carbon and nutrient cycling in salt marsh systems. While decomposition within the salt marsh is conducted primarily by sulfate reducers (Howes et al. 1984), most of the alternative pathways for anaerobic respiration, fermentation, and carbon fixation pathways have also been identified in marsh sediments (Baker et al. 2015; Vineis et al. 2023; Bulseco et al. 2020) and we are only beginning to understand how many of the carbon cycling processes are linked among organisms and within individual genomes.

Below the rooting zone and beyond the influence of living plants, salt marshes have accumulated vast stores of organic peat that date back, in some cases, millennia (Kirwan et al. 2023). These older sediments are likely to contain microbes adapted to survive in the deep, low energy environment where prior oxidation of organic material has led to a decline in the bioavailability of organic carbon with sediment age and depth (Arndt et al. 2013). Microbial communities below the influence of living plants therefore might mirror other marine low energy sediment communities supported by complex organic matter derived from compounds such as lignin, cellulose, and phenolic polymers (Jørgensen and Boetius 2007). Decomposition of these compounds within energy limited sediment is often dependent on syntrophic interactions among microbial consortia (Peng et al. 2023) The specific mechanism of microbial decomposition has important consequences for carbon stability because microbial turnover can result in the production of carbon that is protected from further decomposition (Liang et al. 2017). Therefore, identifying the capacity for microbial decomposition is important to understanding the susceptibility of long-term carbon storage within salt marsh sediments.

The physical and biogeochemical differences between the dynamic upper layers of the sediment, the energy poor deeper layers, and the transition between the two environments, are likely to support distinctly different microbial communities. We hypothesize that communities within dynamic upper layers will contain genomes with diverse functional potential. Additionally, we hypothesize that the energy limited deeper regions will contain microbial consortia assembled to efficiently metabolize complex carbon. We aim to identify the boundaries of both metabolically flexible genomes and those relying on metabolic handoffs. We also consider the possibility that compositional changes in these low energy sediments will mirror the burial processes observed in marine environments that contain similar energetic limitations (Kirkpatrick et al. 2019; Walsh et al. 2016).

To test these hypotheses, we evaluated the structure and functional potential of microbial communities within 5-240 cm depth fractions from two S. patens salt marshes, using a metagenomics and genome reconstruction. Initially, we evaluated the phylogeny, distribution, and co-occurrence of reconstructed draft genomes. We also characterized the potential genomic capacity to decompose organic carbon. Our initial observations of co-occurrence and functional potential prompted us to employ genome-scale metabolic models (GEMs) among co-occurring microbial genomes to identify potential functional complementarity important to carbon cycling.

Methods

1. Sampling

Six sediment cores were taken using a Russian peat corer from the salt marsh platform at two sites in Rowley MA, USA (42.759 N, 70.891 W) primarily vegetated by Spartina patens. The sites are part of the Plum Island Ecosystem (PIE) Long-Term Ecological Research site (LTER). Three cores were collected from the high marsh Spartina patens platform at two adjacent creeks, West and Sweeney Creek (Fig. 1A). The two creeks are approximately 2 km apart and cores collected within each creek were separated by less than 50 meters. Cores were sampled up to the point of resistance, reaching a maximum depth of 240 cm. Once extracted, subsections of the core were collected approximately every 10 cm and a 2 cm sub-section was homogenized in a sterile 50 mL falcon tube, flash frozen using liquid nitrogen and then stored at -80°C for nucleic acid analysis. A detailed description of sample handling and analysis, including biogeochemistry of the cores has been previously described (Bulseco 2018).

Figure 1.

2. Metagenomic Sequencing

The input DNA for sequencing library construction was extracted from homogenized sediment from each 2 cm thick fraction (0.25 g wet weight) using a Qiagen Power Soil DNA extraction kit (Qiagen, Germantown, MD) according to the manufacturer’s recommendations. We evaluated the quality and the quantity of DNA using a Nanodrop and Quant-IT DNA assay (Invitrogen, USA) respectively. Purified DNA was shipped overnight to the Joint Genome Institute (JGI) on dry ice. At JGI, libraries were created using paired-end (PE) library construction with an insert size of 250 bp. Each library was uniquely barcoded and sequenced according to PE 2x151 chemistry on an Illumina NovaSeq. Sequences were quality filtered using the default settings of BBTools v38.26 at JGI (https://jgi.doe.gov/data-and-tools/bbtools/). Quality filtered sequences of each sample were independently assembled using SPAdes assembler 3.12.0 (Bankevich et al. 2012).

3. MAG reconstruction, dereplication, phylogenomics and metagenomic read recruitment

We used assemblies and quality filtered sequences generated by JGI/IMG for each of the 55 samples to reconstruct metagenome-assembled genomes (MAGs) for each independent sample (Fig. 2B). The assembled MAGs represent a population of cells. Therefore, we refer to MAGs as populations or population genomes. Genome reconstruction efforts can be improved by co-assembly of multiple samples in a gradient or time series, especially when read counts are not sufficient to provide the minimum amount of coverage necessary for assembly (Stewart et al. 2018). However, the detection of redundant genome populations across samples requires that genome reconstruction be carried out for each sample independently. Additionally, the number of high-quality genomes recovered can be improved by reconstructing genomes from individual samples compared to co-assembly of multiple samples (Olm et al. 2017).

Figure 2.

We reconstructed draft genomes for each sample using CONCOCT (Alneberg et al. 2014) as an automated guide for genome binning, followed by manual curation of each CONCOCT MAG in Anvi’o, using contig tetranucleotide frequency, coverage, taxonomic assignment, completion, and contamination to guide this process (Eren et al. 2015; 2020). In cases where the placement of contigs could not be resolved among a group of similar MAGs, we used anvi-refine to provide increased resolution on this subset of similar MAGs. This allowed for the separation of closely related MAGs and/or the improvement of draft genome bins, as determined by Anvi’o real-time completion and contamination scores.

Dereplication of the MAGs and selection of the representative MAG from the dereplicated collections were conducted using the default settings for dRep (Olm et al. 2017). Default settings in dRep include dereplication of genomes with an average nucleotide identity (ANI) of 0.95 and a minimum overlap of 10% between the two genomes. dRep incorporates additional essential methods, including CheckM (Parks et al. 2015) and FastANI (Jain et al. 2018). dRep uses completion and contamination based on checkM, which were never at or above the 75% threshold for any of the Patescibacteria. Patescibacteria MAGs were dereplicated separately without a checkM quality score filter. We used the GTDB-Tk to assign taxonomy to each of the MAGs recovered using default settings (Chaumeil et al. 2020).

To identify the phylogenomic relationships among the collection of dereplicated MAGs, we first exported the amino acid sequence of the best hit for 71 single copy genes (Bacteria_71). The amino acid sequences for these genes were aligned using MUSCLE and the alignment was trimmed using trimAL v1.4 (Capella-Gutiérrez et al. 2009). We reconstructed the phylogenetic tree using the concatenated sequence alignment according to anvi-gen-phylgenomic-tree which employs the FastTree tree building algorithm (Price et al. 2009). The fasta file of single copy gene sequences, alignment, tree, and additional details of phylogenetic reconstruction are available in this github repository https://github.com/jvineis/DEEP-CORE-MANUSCRIPT.git

We mapped the metagenomic reads from each sample to the collection of dereplicated MAGs using bowtie2 v2.5.3 (Langmead and Salzberg 2012). Unmapped reads were filtered from the results and samtools v1.19.2 (Li et al. 2009) was used to tabulate the number of reads recruited to each MAG. The relative abundance of each MAG was estimated by dividing the number of reads mapped to the MAG by the number of reads recruited to the non-redundant collection of MAGs.

4. Identifying co-occurring MAGs using network reconstruction

To identify groups of co-occurring genome populations that may be indicative of a shared niche space and/or burial trajectory, we constructed a molecular ecological network (MEN) using random matrix theory (RMT) according to the procedures recommended by iNAP (Feng et al. 2022). RMT is more commonly used for functional or phylogenetic marker genes but is broadly designed to identify co-occurrence patterns among complex datasets that contain large numbers of variables (species occurrence, microarray data, marker genes) (Gibson et al. 2013). The count matrix of reads recruited to each of the 377 dereplicated MAGs was used as the input for RMT MEN construction. Following log transformation of the count matrix, an adjacency matrix was calculated according to Pearson correlation coefficients. To separate significant co-occurrence from noise, RMT estimates the nearest neighbor spacing distribution (NNSD) within the entire correlation matrix. A random matrix of NNSD values will follow a Gaussian distribution and a non-random matrix will follow a Poisson distribution (Luo et al. 2007). To identify the threshold where the correlation estimates transition from Gaussian to Poisson, we ran iterations of successively smaller correlations starting at 1 with a step size of 0.01. A Chi-square test was used to evaluate transition from a Poisson to a Gaussian distribution at each step and iNAP created a plot of resulting p-values to evaluate each step. Visual inspection of the cutoff indicated that correlations below 0.92 were more likely due to random noise (Fig. S1). All correlations below this 0.92 cutoff were set to zero and a network was constructed. RMT network properties were calculated for the resulting filtered network, including R-square of the power-law relationship. A network containing few MAGs with large numbers of connections with other MAGs will have connectivity that follows a power-law distribution. If the connectivity doesn’t fit the power-law distribution, this indicates that there are many MAGs with high connectivity to other MAGs (a very dense network). Module separation and identification of module hubs were calculated according to a greedy modularity optimization (Newman 2004). Within-module (Z) and among-module connectivity (P) were used to classify each MAG as a network hub (Z > 2.5 and P > 0.62), module hub (Z > 2.5 and P < 0.62), connector (Z < and P > 0.62), or peripheral (Z < 2.5 and P < 0.62) (Deng et al. 2012). Based on previous ecological interpretation of networks, the peripheral, connector hub, module hub, and network hub might represent a continuum from specialist to generalist in that order (Olesen et al. 2006). For our purposes of identifying co-occurring genomes, we focused on the identification of modules, which are interlinked subsets of MAGs, and connector hubs that could be important to the stability of the different modules (Olesen et al. 2006).

5. Metabolic potential and distribution of MAGs within identified network modules

To estimate the differences in functional potential among the RMT modules, we constructed a matrix composed of KEGG functional IDs obtained using anvi-estimate-metabolism. MAGs that were included in modules containing fewer than four members or not assigned to a module were also removed from the matrix. We constructed an Atchinson distance matrix and built an NMDS ordination using the vegdist and metaMDS functions respectively. We tested the influence of RMT module and taxonomic classification at the phylum level on the community turnover using the envfit function. All diversity analysis was conducted in R (R Core Team 2024) with functions contained in the vegan package (Oksanen et al. 2024). Relative abundance of MAGs within each module were visually inspected to identify their distribution within the vertical sediment profile.

6. Connector hub (Bathyarchaeia BA1) subnetwork genome scale metabolic networks

Visualization of the complete RMT network (described above) using Cytoscape (Shannon et al. 2003) led to the identification of four modules with similar taxonomy, functional profiles, and high relative abundance within deeper sediment layers. The MAGs identified as connector hubs common to all four modules, Bathyarchaeia BA1, was used to identify a subnetwork of MAGs that were positively connected and therefore contained the potential for metabolic complementarity. Focusing on MAGs in this subnetwork, we estimated the metabolic potential for key pathways related to decomposition of organic compounds, aromatics, acetate, hydrogenases, and other pathways important to carbon cycling within low energy sediment environments using METABOLIC according to the default parameters (Zhou et al. 2022). We reconstructed genome-scale metabolic networks (GSMNs) using gapseq (Zimmermann et al. 2021) and miscoto (Frioux et al. 2018) was employed to screen GSMNs for key members of minimal communities required to produce the metabolites of the entire community. Analysis of GSMNs was employed using Metag2Metabo (M2M) (Belcour et al. 2020). This approach identifies putative metabolites produced by the co-occurring members of the network module, the key microbial taxa required to produce these metabolites, and interchangeable (redundant) members of the network module.

Results

1. Diverse MAG populations were broadly distributed within salt marsh sedimentary layers and both sites

We generated 645.9 Gbp sequences from 55 samples collected up to ∼240 cm of salt marsh sediment from six cores (Fig. 1B). The manual genome binning and refinement effort for each of the 55 samples recovered a total of 2,346 MAGs. Dereplication resulted in a collection of 377 MAGs that included 44 phyla, with novel MAGs identified at the class, family, order and genus level (Table S1). The dereplicated MAGs were an average of 87% of expected single copy genes (completion) and 2.7% of these genes found more than once (contaminated) (Fig. 2, Table S1).

Most of the MAGs recovered represented novel families within Bacteria and Archaea and were well distributed within the layers of sediment (Fig. 2) and dereplication revealed that many of the MAGs were independently recovered from multiple metagenomic samples (Fig. 3). Phyla containing the largest number of MAGs included Chloroflexota, Planctomycetota, Alphaproteobacteria, Gammaproteobacteria, and Desulfobacterota. A large fraction of the Chloroflexota MAGs were Dehalococcoidia, which were primarily detected in sediment layers below 70 cm and Anaerolineae which had a shallower distribution (Fig. 2). Nearly all 41 MAGs in the Planctomycetota were primarily located below 40 cm. MAGs classified as Desulfobacterota were primarily found in shallow samples, but a distinct clade was detected primarily within shallow sediment. The 39 Proteobacteria MAGs were split evenly between Alphaproteobacteria with MAGs that span the entire depth of sampling and Gammaproteobacteria that were mostly detected in sediment above 40 cm (Fig. 2). Similarly, the 28 MAGs within the phylum Bacteroidota were partitioned according to depth, although not as strictly as Alpha and Gammaproteobacteria. Within the Bacteroidota, families Vadin HA17, UBA10428, and UBA5072 were generally found below 30-40 cm while Ignavibacteriaceae were found at 10 cm.

Figure 3.

The number of MAGs contained in each dereplicated cluster ranged from 1 to 30 (Fig. 3). Some examples of genome populations represented by a large cluster of dereplicated genomes include 25 independently reconstructed Ca. Aminicenans that were reduced to two unique MAGs, 63 Actinobacteriota reduced to 4 MAGs, 52 Bipolaricaulia reduced to 6 MAGs, 85 WOR3-B reduced to 11 MAGs, and 24 TA06 reduced to 1 unique MAG (Fig. 3). Many of the MAGs were recovered independently from more than one depth and often from different cores at locations adjacent to different creeks (Fig. 3). In some instances, MAGs from the same dereplicated cluster were independently reconstructed from the top 10-20 cm and below 200 cm. For example, a dereplicated WOR3 MAG contained nearly 30 independently assembled MAGs from samples that spanned the entire 10-240cm vertical profile (Fig. 3). The Desulfobacterota, Bacteroidota, Actinobacteriota, Planctomycetota, and Chloroflexota each contained several MAGs that were recovered from all depths. Sparsely detected Proteobacteria, Gemmatimonadota, and Asgardarchaeota were recovered predominantly from shallow sediment layers, 90cm and above (Fig. 3).

The proportion of metagenomic sequences recruited (mapped) to our collection of nonredundant MAGs was generally greater in sediment samples collected below 40 cm (Fig. S2). Read recruitment below 40 cm was at least 35% and reached nearly 60% in most of the samples below 90 cm. Conversely, metagenomic sequences from sediment samples collected from above 20 cm recruited less than 20% of sequences per sample (Fig. S2) and the number of MAGs recovered was generally lower.

2. Co-occurring MAGs

2.1 Strong positive correlations exist among groups of functionally similar MAGs in the deeper sediment layers

Based on MAG relative abundance, we identified a connected network containing 1547 positive connections (edges) among 226 of the 377 MAGs and zero negative edges (Fig. 4A). Most of MAGs were peripheral nodes in the network, but one MAG (Dehalococcoidales – UBA5760) was identified as a module hub and 15 MAGs were classified as connector hubs (Fig. S3). Among the connector hubs were two Bathyarchaeia classified as BA1, two Dehalococcoidales, Desulfatiglandales, and GIF9 MAGs. The R² of the power-law was 0.97, indicating a scale-free network where most MAGs in the network have few connections. We chose to focus our functional analysis on the top seven modules identified in the network because the remaining 11 modules contained fewer than 5 MAGs (Table S1).

Figure 4.

Modules 1 and 2 contained the largest number of MAGs with a high degree of connectivity within and between the two modules (Fig. 4A). Modules 1 and 2 were also diverse, containing MAGs from 14 and 23 phyla respectively. MAGs within module 7 contained many positive links with modules 1 and 2 and the functional potential among MAGs within these three linked modules was very similar (Fig. 4B, Fig. S4). However, modules 1 and 2 contained larger dispersion in the ordination than module 7 due to the presence of Patescibacteria, which were functionally distinct from other MAGs, largely due to their limited functional potential. Most of the Patescibacteria contained fewer than 800 genes (Table S1). A relatively large proportion of MAGs in modules 1, 2, and 7 (> 25%) contained pathways for complex carbon decomposition, fermentation, amino acid metabolism, and several C1 pathways (Fig. S4). However, pathways for carbon fixation, nitrogen and sulfur reduction, and oxidative phosphorylation were identified in a relatively lower proportion relative to modules 3 and 6 (Fig. S4). Connectivity among modules 4 and 5 was sparse in comparison to the large number of positive connections among modules 1, 2, and 7. However, the overall functional potential was similar among all five modules (1,2,4,5,7) (Fig. 4B, S4).

2.2 Modules containing widespread and abundant MAGs in sediment layers

The relative abundance of MAGs at the phylum level within a module was consistent across sediment layers (Fig. 5). For example, the top four most abundant phyla in module 1 included Chloroflexota, Planctomycetota, Bacteroidota, and Actinobacteriota. In module 2, the cumulative relative abundance of Choloroflexota, Planctomycetota, TA06, and WOR3-B were also among the most abundant phyla within this module and were found in similar relative proportions in the sediment layers.

Figure 5.

Despite the diversity of MAGs assigned to a particular phylum within a module, a small number of MAGs were often the most abundant at all sediment depths. For example, 12 of the 70 MAGs in module 1 represented the majority of reads recruited (Fig. S5). The abundant MAGs in the module were distributed across 6 of the 13 phyla in the module. Within module 1, 5 Chloroflexota MAGs were among the most abundant (Fig. S5). Other phyla within the module were low diversity, such as TA06 in module 2, where a single MAG was the sole contributor to the phylum level abundance patterns.

Although often detected in the entire depth profile, MAGs within modules 1, 2, 4 and 7 were especially abundant below 20 cm (Fig. 5, S4). Modules 1 and 2 contained several Chloroflexota MAGs that collectively recruited at least 15% of the mapped reads. Average coverage of these MAGs was over 50X in the majority of samples below 40 cm. Module 2 contained two of the most abundant MAGs below 40 cm, TA06 and Bipolaricaulota, which recruited nearly 15% of mapped reads in many samples below 100 cm. The peak relative abundance of MAGs in modules 1, 2, 5, and 7 were distinct, with module 1 peaking in abundance at 70-110 cm, module 2 peaking in all samples below 100 cm, module 5 peaking at 30 - 40 cm, and module 7 peaking at 150 – 250 cm (Fig S4).

2.3 Unconnected MAGs are dominant in the upper 40 cm and potentially have a niche in the deepest sediment

In addition to MAGs that were associated with the modules described above, there were numerous unconnected MAGs that were predominantly located in sediments at or above 40 cm (Fig. S4). Proteobacteria (Pseudolabrys and Methyloceanibacter) and Crenarchaeota (Nitrososphaeraceae and Bathyarchaeia BA1) MAGs were among the most abundant unconnected MAGs. However, we also observed MAGs within many sediment layers that were abundant both in the upper 40 cm and below 200 cm within the core collected from West Creek (WE3) (Fig. S6). Although they maintained a patchy distribution in the surface, the relative abundance was consistent at 200-240 cm (Fig. S6). Methyloceanibacter MAGs identified in this “bimodal” group are putative bacterial methylotrophs (Vekeman et al. 2016), but none of the mags contained the methane monooxygenase or co-methyltransferase genes. Bimodal Nitrososphaeraceae MAGs contained ammonia oxidation potential (amo) and the Pseudolabrys contained sulfur oxidation (sox), sulfur reduction (dsr, apr, sat), and nitrate reduction potential (narG) (Fig. S6). These functions suggest a mix of aerobic and anaerobic respiration potential within the deepest layers sampled in the West creek sediment core.

2.4 Bathyarchaeia BA1 are potentially important facilitators of abundant deep sediment communities

Because of the potentially important ecological role of Bathyarchaeia BA1 MAGs in the network as connector hubs (Fig. S3) and global significance of this family to carbon cycling (Hou et al. 2023), we focused on these MAGs and the connections to other MAGs in the network to understand the potential relevance to carbon cycling. Our network analysis identified 2 separate Bathyarchaeia BA1 MAGs that were connector hubs in the RMT network (Fig. S3). To explore the potential role of these Bathyarchaeia BA1 we created a subnetwork of MAGs that are directly connected to the Bathyarchaeia MAGs, resulting in a total of 33 MAGs in the subnetwork (Fig 5A). Among the 12 different orders contained in the subnetwork were members of Bipolaricaulia, GIF9, WOR3, Desulfatiglandales, and Caldatribacteriota and each of the four modules (1, 2, 5, and 7) were represented (Fig. 6A). The relative abundance of members at the phylum level in this subnetwork was also consistent across many of the samples (Fig. 6B).

Figure 6.

3. Putative functional complementarity among Bathyarchaeia subnetwork MAGs suggests potential decomposition of hydrocarbon and aromatics

3.1 Consistent relative abundance among members of the Bathyarchaeia subnetwork

Within the Bathyarchaeia subnetwork, GIF9 MAGs were the most abundant, totaling 8 MAGs, all of which are putative acetogens based on their functional profile. Six GIF9 MAGs were located in the same module (module 2), indicating a close ecological connection and yet they were not identified as functionally redundant members of the community according to the M2M analysis. The collective GIF9 MAG relative abundance reached 15% in many of the samples below 20 cm, but like all MAGs in the subnetwork, they were not excluded from the upper 20 cm (Fig. 6B). Desulfatiglandales, Dehalococcoidales, and Planctomycetota were also among the most abundant MAGs in the subnetwork and occurred at similar cumulative relative abundances. The high correlation among these MAGs detected in the network is reflected in the consistent relative abundance of the group members within a core and generally across all six cores (Fig. 6B).

3.2 Bathyarchaeia BA1 subnetwork MAGs use cooperative metabolism to degrade aromatics

M2M identified the putative functions of each community member and the added value metabolites produced by minimal communities. An added value metabolite is one produced through complementary metabolism among the microbes in a minimal community, and a minimal community is defined as a cohort of MAGs required to produce the entire set of metabolites in the collection of genomes. We found 82 metabolites that were producible by each of the 33 MAGs in the subnetwork (core functions) and 1098 total metabolites could be produced by the individual organisms all together (collective functions). The added value of cooperation within the minimal communities required 22 bacteria and 20 essential symbionts to produce an additional 228 compounds. These 228 compounds were only possible when all members of a minimal community had the potential to exchange metabolites. Redundant members of the community that were identified in at least one, but not all minimal communities included Dehalococcoidia – RBG-16-60-22, Caldatribacteriota JS1, Dehalococcoidia AB-539-J10, and Dehalococcoidia UBA5760.

Among the 228 compounds that could be formed solely if organisms exchanged compounds, key metabolites involved in anaerobic decomposition of aromatic compounds were identified and the presence of the genes conferring capacity for these putative functions were confirmed by METABOLIC. Benzoyl-CoA reductase (bcr) was identified in the genomes of three Desulfatiglandales, a novel order of Cloroflexota, and Bathyarchaeia. Genes to produce benzoyl-CoA from phenol (ubiX/bsdC) occurred in two of the Desulfatiflandales, the Bathyarchaeia and Bacteroidales MAGs (Fig. 6C). Phenylacetyl-CoA was also identified as a metabolite putatively produced through functional complementarity within the subnetwork.

3.3 Genomic potential for acetogenesis and metabolic versatility is common among the Bathyarchaeia BA1 subnetwork

Genes for heterotrophic and autotrophic acetogenesis, primarily through the Wood-Ljungdahl (WL) pathway, were detected in 25 of the 33 genomes in the Bathyarchaeia B26-1 subnetwork (Fig. 6C). These genes included acetyl-CoA dehydrogenase (acdA), acetate kinase (ack), and phosphate acetyltransferase (pta). acdA was the most commonly detected gene indicating acetogenic metabolism, and acs (Acetyl-CoA synthetase) was commonly detected within the same MAG (Table S2). The co-occurrence of genes encoding a putative WL pathway within most genomes and a formate oxidation gene (fae) required to convert CO₂ into formate provides additional evidence for synthesis of acetate from H₂ and CO₂. The genomic potential for the production of acetate enabled by carbon fixation via the WL pathway was observed in Desulfatiglandales, Bathyarchaeia, Dehalococcoidales, Chloroflexota, GIF9, Bipolaricaulia, and Planctomycetota MAGs in the subnetwork (Fig. 6C). A WL pathway was not detected in Caldatribacteriota, WOR3, and Bacteroidales. Instead, these MAGs might produce acetate through fermentation of pyruvate or lactate. Most acetogens in the network contain genes indicative of a flexible metabolic capacity through autotrophic and/or fermentative acetate production (Fig. 6C).

We detected several additional pathways potentially related to energy conservation among the members of the subnetwork. The presence of methane monooxygenase regulatory protein (mmoBD) in Bipolaricaulia and Desulfatiglandales (Fig. 6C) suggests methane oxidation potential. Additionally, aerobic carbon monoxide dehydrogenase coxSM was detected in members of Bipolaricaulia, Caldatribacteriota, Bathyarchaeia TCS64, and Desulfatiglandales. Additional pathways include selenate or arsenate reduction potential in most of the subnetwork MAGs and one of the GIF9 MAGs contained genes for both pathways.

Discussion

Salt marsh sediments contain diverse microbial communities with vast numbers of uncultivated microbial organisms that are important to carbon storage and biogeochemical cycling. Within the complex microbial communities inhabiting these sediments, we identified metabolic complementarity important to cycling pools of aromatic carbon among co-occurring groups of microbes. Bathyarchaeia BA1 was a key member of each group of MAGs that were uniquely distributed within the sediment profile, and they were primarily located below the dynamic upper layers of the sediment. Another key finding was that similar genome populations exist throughout the entire depth profile of multiple cores, spanning nearly 3000 years of sediment accumulation (Bulseco 2018). This pattern may be the result of a burial process followed by environmental selection and slow growth rates. Conversely, the upper layers of sediment contain additional diverse MAGs with greater functional potential representing the canonical pathways commonly associated with salt marsh microbial communities (Baker et al. 2015; Tobias and Neubauer 2019). Our results are important for understanding how the pools of organic carbon are shaped by microbial communities within the less well characterized subsurface because we identify which microbes are co-located, what forms of carbon may be made accessible due to metabolic complementarity, and the underlying biogeochemical processes that preserve these relationships.

1. Potential metabolic flexibility is a key feature among MAGs in the zone of plant and tidal influence

Plant roots and rhizomes along with bioturbating organism and other structural features can increase permeability in the upper layers of the sediment (Guimond et al. 2020; Moffett et al. 2012; Xiao et al. 2019) and these dynamic conditions potentially lead to observed high diversity of microbial communities in the upper sediment layers. Redox conditions and carbon availability can change rapidly in the upper sediment layers, creating conditions where metabolic flexibility would be favored. As a result, there is a decreased need for metabolic handoffs among members of the community.

Most of the MAGs that were unconnected in our network analysis were localized in the dynamic upper layers of the sediment and did not occur below 90 cm. The unconnected MAGs contained striking genomic versatility including the capacity to utilize multiple electron acceptors for oxidation of diverse carbon sources, conduct fermentation, and fix inorganic carbon. Many of the unconnected Proteobacteria and Desulfobacterota MAGs contain the genomic potential for each of these metabolic pathways in a single organism which was observed previously in similar sediments (Vineis et al. 2023). This metabolic flexibility allows each of the members of the community to inhabit a variety of conditions and this strategy makes sense given the dynamic nature of these sediments. Members of these communities contain canonical pathways that are commonly associated with carbon cycling in salt marsh sediments including fermentation, denitrification, and sulfate reduction (Baker et al. 2015).

2. Organic carbon decomposition in deeper sediment layers is potentially facilitated by co-occurring community members with complimentary metabolic pathways

The consistency in relative abundance among members within the core Bathyarchaeia subnetwork is a potential result of metabolic handoffs. Of particular significance is the detection of benzoyl-CoA as a metabolite that could only be produced by Bathyarchaeia, Dehalococcoidia, and Desulfatiglandales through metabolic handoffs in the community. Benzoyl-CoA is a known intermediate involved in the anaerobic decomposition of aromatic compounds toluene, phenol, ethylbenzene and benzoate (Carmona et al. 2009; Porter and Young 2014). Under anaerobic conditions, reduction of benzoyl-CoA leads to the production of acetyl-CoA which can enter into central metabolic pathways for growth and energy conservation (Porter and Young 2014). Phenylacetyl-CoA was also putatively produced through metabolic complementarity and is similar to benzoyl-CoA, leading to the potential production acetyl-CoA/succinyl-CoA. These metabolic pathways, central to aromatic decomposition are only possible due to the potential metabolic handoffs among the subnetwork members. Although it is possible that other combinations of MAGs could lead to the decomposition of aromatics, co-occurrence provides additional support that handoffs among the genome populations identified here are likely to occur in the environment.

Syntrophic interactions between a bacterial fermenter and CO₂ reduction are often cited as important to the decomposition of aromatic carbon (Dai et al. 2024; Qiu et al. 2006). The collective functional profile of the Bathyarchaeia subnetwork suggests that the initial fermentation of complex organic carbon is conducted by members of the Planctomycetota, WOR3, JS1 and others, followed by secondary decomposition of aromatics by Bathyarchaeia, Deslufatiglandales, Chloroflexota. The products of secondary fermentation, including H₂ and formate, usually used by methanogens, are instead potentially being used by a diverse group of acetogens. Acetogens grow on organic materials derived from cellulose, lignin, alcohols, hemicellulose, and chitin within energy depleted sediments (Lever 2012; Yu et al. 2018), and there is evidence for the use of lignin as an energy source in Bathyarchaeota (Yu et al. 2018).

These results suggest a central mechanism for carbon cycling among closely coupled consortia within the sediment layers of salt marsh sediment. While this process may be very slow, we have evidence from amplicon analysis of cDNA, that these communities are active (Bulseco 2018). Conversion of these aromatics into microbial biomass suggests that there may be a relevant cycling and/or recycling of carbon stabilized by microbial activity even within these energy depleted sediments.

3. Burial as a mechanism for aggregating putative consortia

Burial as a mechanism for co-occurrence patterns has been suggested in energy starved deep marine sediments where selection on the microbial community is thought to result in low diversity (Kirkpatrick et al. 2019; Walsh et al. 2016). Previous studies of microbial turnover within salt marshes described seasonal and transect community compositional changes (Tebbe et al. 2022; Yao et al. 2019) or succession along a chronosequence (Dini-Andreote et al. 2014; Wang et al. 2017). However, the composition and stability of microbial communities within salt marsh sediments remains relatively uncharacterized at the depths examine here.

Our study of the vertical sediment profile recovered two pieces of evidence that suggests burial is a likely mechanism for the distribution of closely associated consortia. The first observed pattern suggestive of burial was the recovery of MAGs assembled independently from multiple depths in the same core that were dereplicated into a single genome population. One of the strongest examples of this was observed in the candidate phylum TA06, where a single population was recovered from 24 of the 55 samples, including every depth fraction except 10 cm. For dereplication to occur at this scale (> 95% identity among MAGs in the dereplicated cluster), genome populations would have experienced limited evolution over the nearly 3000 yrs of separation that exists between the surface and bottom of the core. The second observed pattern was the consistent relative abundance of MAGs within each module, relative to each other at the phylum level (Fig. 5). The consistent relative abundance patterns among MAGs within a module was observed both within the vertical sampling of an individual core and across multiple cores. This is not to say that the relative abundance of members within a module are the same at every depth across all cores. However, when the MAGs were identified in a sample, they generally appear in similar ratios. The ratios of microbial taxa might reflect the proportions required to support syntrophic interactions.

Similar compositions of sediment communities in marine systems have been recovered from geographically distant communities (although at high taxonomic levels) and were presumed to be subseafloor specialists (Inagaki et al. 2006; Petro et al. 2017; 2019). However, our data suggest that the same populations were buried within the sediment, essentially preserving the population through high selection coupled with slow growth.

The low porosity, small particle size, and low circulation observed in salt marsh sediments on the East Coast of the US (Guimond and Tamborski 2021) provide additional support to co-occurrence shaped by burial that we observed in our genomic data. The turnover of water within the sediment can require 30-90 years (Guimond and Tamborski 2021), depending on the marsh system, indicating that water flow from groundwater or surface sources may be limited. These physical characteristics present significant barriers to the dispersal of subsurface communities, leaving burial as a likely mechanism. One caveat to the low porosity and water flow are the observations in the deepest sediment layers collected from West Creek where we observed members of the community that require oxygen or terminal electron acceptors commonly associated with anaerobic shallow sediments. We speculate that the groundwater may potentially “well up” from below the bedrock or sand located beneath the salt marsh platform at this site.

Conclusion

The metabolic connections identified among correlated members of the microbial community provide key insights into how carbon is cycled by microbial communities within salt marsh sediments. Additionally, these results provide evidence for the mechanisms important to understanding the assembly and distribution of these communities. Our analysis identifies syntrophy among members of a putative core consortia that provides a potential mechanistic understanding of organic matter decomposition. The process begins with decomposers of complex carbon, followed by aromatic decomposition, coupled with acetogenic energy conservation and finally oxidation of acetate through arsenate or selenate reduction. These findings represent a significant step forward in our understanding of how microbial communities survive within these energy limited systems and contribute to the stability of soil organic carbon in this blue carbon system. The persistent conversion of soil organic carbon by the microbial community indicates that the pools of carbon potentially originating from plant material undergo decomposition and transformation even under energy depleted conditions.

Data Availability

The data described here are described in the U.S. Department of Energy (DOE) Facilities Integrating Collaborations for User Science (FICUS) grant titled “Combining high resolution organic matter characterization and microbial meta-omics to assess the effects of nutrient loading on salt marsh carbon sequestration” (proposal ID = 503576, https://genome.jgi.doe.gov/portal/Comhiguestration/Comhiguestration.info.html). Supplemental data can be accessed using this link to a publicly available github page https://github.com/jvineis/DEEP-CORE-MANUSCRIPT.git

Acknowledgements

Samples were collected from the Plum Island Ecosystems LTER, which is supported by NSF (OCE 2224608) and from the NSF funded TIDE project (DEB 1902712). Sequencing was conducted at the Joint Genome Institute through a funded Facilities Integrating Collaborations for User Science (FICUS) grant (proposal ID 503576). Additional funding for this project was provided by an NSF PRFB to ANB (1907285), Simons Foundation (Grant #1247132 to ZGC, and Grant #01248016 to JLB), and Department of Energy (DOE DE-SC0024270).

Additional files

Supplemental

Additional information

Funding

United States Department of Energy (503576)

United States Department of Energy (DE-SC0024270)

National Science Foundation (OCE 2224608)

National Science Foundation (DEB 1902712)

National Science Foundation (PRFB 1907285)

Simons Foundation (1247132)

Simons Foundation (01248016)