Cross-Feedings, Competition, and Positive and Negative Synergies in a Four-Species Synthetic Community for Anaerobic Degradation of Cellulose to Methane

ABSTRACT Complex interactions exist among microorganisms in a community to carry out ecological processes and adapt to changing environments. Here, we constructed a quad-culture consisting of a cellulolytic bacterium (Ruminiclostridium cellulolyticum), a hydrogenotrophic methanogen (Methanospirillum hungatei), an acetoclastic methanogen (Methanosaeta concilii), and a sulfate-reducing bacterium (Desulfovibrio vulgaris). The four microorganisms in the quad-culture cooperated via cross-feeding to produce methane using cellulose as the only carbon source and electron donor. The community metabolism of the quad-culture was compared with those of the R. cellulolyticum-containing tri-cultures, bi-cultures, and mono-culture. Methane production was higher in the quad-culture than the sum of the increases in the tri-cultures, which was attributed to a positive synergy of four species. In contrast, cellulose degradation by the quad-culture was lower than the additive effects of the tri-cultures which represented a negative synergy. The community metabolism of the quad-culture was compared between a control condition and a treatment condition with sulfate addition using metaproteomics and metabolic profiling. Sulfate addition enhanced sulfate reduction and decreased methane and CO2 productions. The cross-feeding fluxes in the quad-culture in the two conditions were modeled using a community stoichiometric model. Sulfate addition strengthened metabolic handoffs from R. cellulolyticum to M. concilii and D. vulgaris and intensified substrate competition between M. hungatei and D. vulgaris. Overall, this study uncovered emergent properties of higher-order microbial interactions using a four-species synthetic community.

ABSTRACT Complex interactions exist among microorganisms in a community to carry out ecological processes and adapt to changing environments. Here, we constructed a quad-culture consisting of a cellulolytic bacterium (Ruminiclostridium cellulolyticum), a hydrogenotrophic methanogen (Methanospirillum hungatei), an acetoclastic methanogen (Methanosaeta concilii), and a sulfate-reducing bacterium (Desulfovibrio vulgaris). The four microorganisms in the quad-culture cooperated via cross-feeding to produce methane using cellulose as the only carbon source and electron donor. The community metabolism of the quad-culture was compared with those of the R. cellulolyticumcontaining tri-cultures, bi-cultures, and mono-culture. Methane production was higher in the quad-culture than the sum of the increases in the tri-cultures, which was attributed to a positive synergy of four species. In contrast, cellulose degradation by the quad-culture was lower than the additive effects of the tri-cultures which represented a negative synergy. The community metabolism of the quad-culture was compared between a control condition and a treatment condition with sulfate addition using metaproteomics and metabolic profiling. Sulfate addition enhanced sulfate reduction and decreased methane and CO 2 productions. The cross-feeding fluxes in the quad-culture in the two conditions were modeled using a community stoichiometric model. Sulfate addition strengthened metabolic handoffs from R. cellulolyticum to M. concilii and D. vulgaris and intensified substrate competition between M. hungatei and D. vulgaris. Overall, this study uncovered emergent properties of higher-order microbial interactions using a fourspecies synthetic community. IMPORTANCE A synthetic community was designed using four microbial species that together performed distinct key metabolic processes in the anaerobic degradation of cellulose to methane and CO 2 . The microorganisms exhibited expected interactions, such as cross-feeding of acetate from a cellulolytic bacterium to an acetoclastic methanogen and competition of H 2 between a sulfate reducing bacterium and a hydrogenotrophic methanogen. This validated our rational design of the interactions between microorganisms based on their metabolic roles. More interestingly, we also found positive and negative synergies as emergent properties of high-order microbial interactions among three or more microorganisms in cocultures. These microbial interactions can be quantitatively measured by adding and removing specific members. A community stoichiometric model was constructed to represent the fluxes in the community metabolic network. This study paved the way toward a more predictive understanding of the impact of environmental perturbations on microbial interactions sustaining geochemically significant processes in natural systems.
C omplex anaerobic microbial communities can degrade cellulose to methane, CO 2 , and other products in many anoxic ecosystems, including wetlands (1), rumen gut (2), waste treatment systems (3), and biogas digesters (4). In these communities, cellulolytic bacteria metabolize cellulosic materials primarily to CO 2 , hydrogen (H 2 ), and a variety of organic acids such as acetate, lactate, formate, and succinate (5). Then, H 2 and CO 2 are converted to CH 4 by hydrogenotrophic methanogens, and acetate is metabolized to CH 4 by acetoclastic methanogens. When sulfate is available, sulfate-reducing bacteria (SRB) may oxidize acetate and lactate to CO 2 and reduce sulfate to H 2 S. When sulfate is limited, some SRB grow as syntrophic oxidizers, fermenting lactate to acetate, H 2 , and CO 2 (6). The overall methane and CO 2 emissions are determined by the nutrient and redox potential of their local environments that eventually drive the interactions among these functional guilds (7)(8)(9). Increased sulfate availability caused by, for example, intrusion of sulfate-laden seawater to estuarine wetlands (10), will likely enable SRB to produce more H 2 S, utilizing H 2 and/or organic acids such as acetate and lactate as electron donors. The competition from SRB for H 2 and acetate may reduce methane production by hydrogenotrophic and acetoclastic methanogens (7,11). Characterization of these microbial interactions will provide insights into community metabolism under changing environmental conditions and allow modeling of metabolic flux rates among community members.
It is a challenge to mechanistically study complex microbial interactions in natural communities. Synthetic communities provide a complementary reductionist system to characterizing microbial interactions from microbiological, ecological, and biochemical perspectives. Synthetic communities can be directly manipulated, and the response of individual community members can be identified (12,13). Previous synthetic community studies have investigated mutualistic interactions between Methanococcus maripaludis and Desulfovibrio vulgaris in bi-cultures (9), cooperative interactions between Geobacter metallireducens and Geobacter sulfurreducens (14), cross-feeding between Escherichia coli and Rhodopseudomonas palustris (15), and competitive interaction in a synthetic community of Pseudomonas aeruginosa and Agrobacterium tumefaciens (16). Rationally designed synthetic communities are experimentally more tractable and enable identifying interactions between key functional guilds that are challenging to assess in environmental communities.
In this study, we developed a simplified, yet functionally complete, synthetic community to study the microbial interactions and carbon exchanges during anaerobic carbon degradation. Our synthetic microbial community was composed of the major functional guilds (cellulolytic fermenter, sulfate reducer, hydrogenotrophic methanogen, and acetoclastic methanogen) that mediate the anaerobic conversion of cellulosic biomass to CH 4 and CO 2 . The choice of a sulfate-reducing bacterium (Desulfovibrio vulgaris Hildenborough) introduced metabolic versatility and enabled investigations into the community response to sulfate intrusion. We characterized the biochemical and physiological response of each microorganism in the synthetic community, investigated the ecological and metabolic functions, and modeled the interspecies interactions governing carbon exchange among these four organisms. This study provides a foundation for the development of mechanistic metabolic models of anaerobic degradation of cellulose to methane, CO 2 , and H 2 S.

RESULTS
Overview of the experimental design. A four-species synthetic community consisting of a cellulolytic bacterium (Ruminiclostridium cellulolyticum), a hydrogenotrophic methanogen (Methanospirillum hungatei), an acetoclastic methanogen (Methanosaeta concilii), and a sulfate-reducing bacterium (Desulfovibrio vulgaris) was constructed. To unravel how each species and their interactions contribute to the community metabolism, the growth of this quadculture was compared with those of a mono-culture, three bi-cultures, and three tri-cultures ( Fig. 1A; Table 1). All cultures were grown anaerobically in triplicate using the same medium with cellulose as the sole carbon source. To simulate the increased sulfate availability from seawater intrusion into a wetland soil community, the D. vulgaris-containing cultures received a daily sulfate addition after the first day of cultivation. The growth status of these multispecies cultures was measured over a week by daily analysis of substrate consumption and product accumulation. The quad-cultures were analyzed with metaproteomics at the end of experiment to characterize the community structure and metabolic activities. A stoichiometric model was built to analyze the cross-feeding fluxes through the four species ( Fig. 1A; see Table S1 and Table S2 in the supplemental material).
Metaproteomics analyses measured the relative proteome abundances of the organisms in the quad-cultures (Fig. 1B) (17,18). A total of 3,736 proteins were identified from the four microbial species (Table S4). R. cellulolyticum was dominant in the quad-cultures as the producer of the growth substrates for the other three microorganisms. Sulfate addition reduced the relative proteome abundance of R. cellulolyticum, on average, from 68% to 59% (P = 0.047), while it increased that of D. vulgaris, on average, from 19% to 28% (P = 0.0004). The relative proteome abundances of M. concilii and M. hungatei remained at 10% and 3%, respectively, in both the control condition and the sulfate addition condition. The fermentation products of the synthetic communities were tracked by carbon molarity (Fig. 1C) and electron molarity (Fig. 1D). The mono-culture of R. cellulolyticum hydrolyzed cellulose to reducing sugars and fermented these sugars to acetate, lactate, ethanol, H 2 , and CO 2 . The quad-cultures transferred 1.37 6 0.03 (29%) mmol of carbon (Fig. 1C) and 10.99 6 0.21 (51%) mmol of electrons ( Fig. 1D) to methane in the control condition. The quad-cultures transferred 9.36 6 0.04 (43%) mmol of electrons (Fig. 1C) to H 2 S in the sulfate addition condition. Sulfate addition reduced  the electrons and carbon transferred to methane by 1.28 6 0.02 mmol and 0.16 6 0.02 mmol, respectively, compared to the control condition.
Metabolism of cellulolytic bacterium Ruminiclostridium cellulolyticum. Cellulose was the only carbon source provided to all the synthetic communities. After 2 days of incubation, all cultures began to rapidly consume cellulose (Fig. S1). Over 7 days of incubation, the R. cellulolyticum mono-culture degraded 184.39 6 0.54 mg of cellulose. Cellulose degradation increased significantly when R. cellulolyticum was cultivated with the other three species (P = 0.0028) ( Table 1). The highest cellulose degradation was achieved by the quad-culture (210.02 6 2.10 mg) ( Fig. 2A).
Metaproteomic analysis of the quad-culture identified 64 R. cellulolyticum cellulase and 50 R. cellulolyticum ABC transporters related to uptake of monosaccharides, disaccharides, and oligosaccharides. Glucose and cellobiose produced by cellulose hydrolysis were detected at micromolar concentrations in the medium (Table 1). In the R. cellulolyticum mono-culture acetate and lactate accumulated in the medium (Fig. 2B), and H 2 and CO 2 accumulated in the headspace (Fig. 2C). These fermentation products from R. cellulolyticum can be used by other microorganisms as carbon and/or energy sources in multispecies cultures. Compared with the R. cellulolyticum mono-culture 29.5% less Accumulation of H 2 and CO 2 in the gas phase. Significant differences (*, P , 0.05; **, P , 0.01; ***, P , 0.001) indicate targeted biological hypotheses. (D) Protein expression profile of the major carbon metabolism pathways in R. cellulolyticum. Enzymes are color-coded based on their protein abundance changes induced by sulfate addition. Unidentified enzymes are not labeled. Sat, sugar-related ABC transporter; Pfo, pyruvate flavodoxin/ferredoxin oxidoreductase; Hyd, hydrogenase, Fe only; Nhd, flavin oxidoreductase/NADH oxidase; Ldh, L-lactate dehydrogenase; Aad, acetaldehyde dehydrogenase; Alr, alcohol dehydrogenase; Ads, AMP-dependent synthetase; Ack, acetate kinase; Pta, phosphate acetyltransferase.
Tertiary Interactions in a Four-Species Community mBio acetate (P = 0.0006), 51.9% less lactate (P = 0.0002), 4.4% less H 2 (P = 0.027), and 30.6% more CO 2 (P , 0.0001) accumulated in the quad-culture control condition without sulfate addition (Fig. 2C). Metaproteomics identified key enzymes in the fermentative pathways in R. cellulolyticum, including hydrogenase (Hyd) for H 2 generation, pyruvate flavodoxin/ferredoxin oxidoreductase (Pfo) for CO 2 production, L-lactate dehydrogenase (Ldh) for lactate generation, and phosphotransacetylase (Pta), AMP-forming acetyl-CoA synthetase (Ads), and acetate kinase (Ack) for acetate generation (Fig. 2D). The protein abundance of L-lactate dehydrogenase from R. cellulolyticum in the quad-culture increased by 3.7-fold in the sulfate addition condition compared to the control condition (q = 0.04). This may have contributed to the 1.8-fold higher lactate accumulation in the quad-cultures in the sulfate addition condition than the control condition (P = 0.003) (Fig. 2B).
Metabolism should be produced by R. cellulolyticum in the bi-culture, which was more than twice that in the R. cellulolyticum mono-culture (1.36 6 0.02 mmol). This indicates that M. hungatei promoted H 2 production by R. cellulolyticum through the alleviation of product inhibition. The methane production increased from 0.16 6 0.01 mmol in the tri-culture without M. hungatei to 1.37 6 0.03 mmol in the quad-culture with M. hungatei (Fig. 3C). Two key enzymes in hydrogenotrophic methanogenesis, formylmethanofuran dehydrogenase (Fmd) and formate dehydrogenase (Fdh), were identified in the quad-cultures by metaproteomics (Fig. 3D). This indicated active hydrogenotrophic methanogenesis by M. hungatei in the multispecies cultures.
M. concilii was added as a model acetoclastic methanogen. The net accumulation of acetate in the multispecies cultures containing M. concilii increased in the first 5 or 6 days and then decreased or leveled off in the last 1 or 2 days, which suggested a lag in acetate consumption (Fig. S2). In the cultures lacking M. concilii, the amount of acetate kept increasing during the 7-day growth. The accumulative amount of acetate in the bi-culture of M. concilii and R. cellulolyticum (0.48 6 0.01 mmol) was significantly lower than in the mono-culture of R. cellulolyticum alone (0.78 6 0.002 mmol) (P = 0.0006). The accumulation of acetate in the quad-culture (0.55 6 0.01 mmol) was also significantly lower than that in the tri-culture that lacked M. concilii (0.71 6 0.02 mmol) (P = 0.0005) (Fig. 3B). Approximately 0.48 mmol more methane accumulated in the quad-culture (1.37 6 0.03 mmol) than in the tri-culture (0.89 6 0.02 mmol) without M. concilii (Fig. 3C) (P , 0.0001). In the M. concilii acetoclastic methanogenesis pathway, acetyl-coenzyme A synthetase (Acs), CO dehydrogenase/acetyl-CoA synthase (Cds), tetrahydromethanopterin S-methyltransferase, methyl-coenzyme A reductase (Mtf), and CoB-CoM reductase (Mcr) were identified by metaproteomics (Fig. 3E). This indicated active acetoclastic methanogenesis by M. concilii in the multispecies cultures.
Sulfate addition reduced quad-culture methane production by 11.7% (P , 0.0001) (Fig. 3C). The protein abundance of formate dehydrogenase, involved in the hydrogenotrophic methanogenesis in M. hungatei, significantly decreased by 27.8-fold (q , 0.0001) (Fig. 3D). In addition, sulfate addition significantly decreased the protein abundances of enzymes involved in acetoclastic methanogenesis. Specifically, the protein abundances of CO dehydrogenase/ acetyl-CoA synthase (q = 0.03) and methyl-coenzyme M reductase (q = 0.001) significantly decreased by at least 4-fold (Fig. 3E). These decreases indicated that the sulfate addition suppressed both hydrogenotrophic methanogenesis and acetoclastic methanogenesis in the quad-culture.
Metabolism of the sulfate-reducing bacterium Desulfovibrio vulgaris. D. vulgaris can ferment lactate to acetate, H 2 , and CO 2 when sulfate is limited (Fig. 4F) or reduce sulfate with lactate or H 2 as electron donors (Fig. 4G) (19). In the control condition, while lactate in the cultures without D. vulgaris accumulated continuously over 7 days, lactate in the cultures with D. vulgaris increased in the first 3 days and then decreased slightly from day 4 onward. Over 7 days, 0.13 6 0.005 mmol lactate accumulated in the bi-culture of R. cellulolyticum and D. vulgaris, which was 51.9% lower than that in the mono-culture (0.27 6 0.001 mmol) (P = 0.0005). The lactate accumulation by the tri-culture that lacked D. vulgaris (0.29 6 0.02 mmol) was also higher than that in the quad-culture that included D. vulgaris (0.13 6 0.004 mmol) (P = 0.002) (Fig. 4A). Metaproteomics identified cytochrome c 3 (c 3 ), high-molecularweight cytochrome (Hmc), periplasmic [NiFe] hydrogenase (Hyd), membrane-bound Coo hydrogenase (Coo), lactate dehydrogenase (Ldh) and L-lactate permease (Lpe) in the hydrogenic lactate oxidation pathway of D. vulgaris (Fig. 4F). In the mono-culture of R. cellulolyticum, 1.36 6 0.02 mmol H 2 accumulated at the end of the cultivation, while 1.52 6 0.003 mmol H 2 accumulated in the bi-culture of R. cellulolyticum and D. vulgaris. These results suggest that D. vulgaris fermented lactate produced by R. cellulolyticum to H 2 . The addition of D. vulgaris to the bi-culture of R. cellulolyticum and M. concilii did not significantly increase the methane production. But the addition of D. vulgaris to the bi-culture of R. cellulolyticum and M. hungatei significantly increased the methane production by ;112% (from 0.42 6 0.007 mmol to 0.89 6 0.02 mmol) (Fig. 4D). D. vulgaris enhanced hydrogenotrophic methanogenesis, likely by producing additional hydrogen from lactate fermentation (9).
We found that the synergy between hydrogenotrophic and acetoclastic methanogenesis depended on the presence of D. vulgaris. Methane produced by the tri-culture of  Table 1). From another perspective to compare, the addition of D. vulgaris to the tri-culture of R. cellulolyticum, M. concilii, and M. hungatei increased total methane production by 0.75 mmol. This was 67% more than the additive increase of adding D. vulgaris to either a R. cellulolyticum-M. hungatei bi-culture (0.47 mmol increase) or to a R. cellulolyticum-M. concilii bi-culture (0.02 mmol decrease) (Fig. 4D). Both ways of comparison showed a positive synergistic relationship among methanogenesis and fermentation via D. vulgaris.
The mono-culture of R. cellulolyticum produced 0.98 6 0.01 mmol of CO 2 . Coculturing R. cellulolyticum with D. vulgaris did not significantly increase their CO 2 production (1.04 6 0.03 mmol). Adding M. hungatei to R. cellulolyticum did not significantly change the CO 2 production in their bi-culture (0.99 6 0.03 mmol) either. However, combining both D. vulgaris and M. hungatei with R. cellulolyticum significantly increased CO 2 production (1.22 6 0.03 mmol). Further, adding M. concilii to this tri-culture did not result in a significant increase of CO 2 production by the quad-culture (1.28 6 0.007 mmol) (Fig. 4E). These results suggested a positive synergy between M. hungatei and D. vulgaris on CO 2 production.
Sulfate addition enhanced sulfate reduction by D. vulgaris. At the end of cultivation, the quad-culture consumed 1.17 6 0.01 mmol of sulfate in the sulfate addition condition (Table 1). Sulfate addition increased the relative abundance of sulfite reductase by 2.0-fold (q = 0.003) (Fig. 4G). H 2 accumulation significantly decreased in the sulfate addition condition (1.07 6 0.02 mmol) compared with the control condition (1.30 6 0.02 mmol) (P = 0.003) (Fig. 4C). Sulfate reduction promoted sulfidogenic lactate oxidation as evidenced by a 1.8-fold increase of the relative abundance of FeS cluster assembly ATPase (Fca) (q = 0.009) in the sulfate addition condition (Fig. 4G). However, lactate accumulated in the first 3 days and then slightly decreased in the following 4 days in the control condition, but lactate continuously accumulated over the 7 days in the sulfate addition condition (Fig. 4A). The increased accumulation of lactate in the sulfate addition condition, despite the enhanced D. vulgaris sulfidogenic lactate oxidation, can be attributed to the increased R. cellulolyticum lactate fermentation stimulated by sulfate addition as described above (Fig. 2D).
Sulfate addition significantly decreased the CO 2 accumulation in the two tri-cultures containing D. vulgaris and the quad-culture, but not in the tri-culture without D. vulgaris ( Fig. 4E; Table 1). Sulfate addition also significantly reduced the acetate accumulation in the two tri-cultures containing D. vulgaris and the quad-culture, but not in the tri-culture without D. vulgaris ( Fig. 4B; Table 1). These observations can be attributed to the suppression of the D. vulgaris fermentation of lactate to acetate and CO 2 by the sulfate addition.
Modeling of carbon and energy flow through the synthetic community. A stoichiometric model was constructed to simulate fluxes through the primary energy production pathways in the four organisms (Tables S1 and S2). For simplicity, extracellular cellulose hydrolysis was omitted, and R. cellulolyticum was considered to ferment glucose. The overall model has seven reactions, including two for R. cellulolyticum, three for D. vulgaris, and one for each methanogen (Tables S1 and S2). Reactions that were not expected to be active were removed based on either genetic or nutrient availability in the system (i.e., sulfidogenic reactions were removed in the absence of sulfate, and acetoclastic methanogenesis was removed in the absence of M. concilii). The contributions of each reaction were fitted to the accumulative millimoles of glucose, acetate, lactate, H 2 , CH 4 , and H 2 S produced or consumed over the 7 days of incubation. To access redundancy of metabolic pathways in this stoichiometric model, each reaction pathway was sequentially excluded before fitting, and the impact on model fit was evaluated. These simulations identified two and three possible flux distributions through the quad-culture in the control and sulfate addition conditions (Tables S1 to Table S3) with equal goodness of fit to the experimental data. The correlation coefficients, R 2 , exceeded 0.9 in both conditions. The normalized root mean square errors (NRMSEs) were 0.11 6 0.006 in the control condition and 0.25 6 0.01 in the sulfate addition condition. The consistency between the experimental data and the modeling results indicated a reliable estimation of the carbon and energy flows (Fig. 5A and B).
R. cellulolyticum in the cultures was modeled by lactate fermentation and hydrogenic acetogenesis (Table S1). Sulfate addition stimulated the production of acetate and lactate in R. cellulolyticum. D. vulgaris was modeled to have three metabolic possibilities, lactate oxidation to acetate with proton reduction to hydrogen (hydrogenic lactate oxidation) (Fig. 4F), lactate oxidation to acetate with sulfate reduction to sulfide (sulfidogenic lactate oxidation) (Fig. 4G), and hydrogen oxidation to protons with sulfate reduction to sulfide (sulfidogenic hydrogen oxidation) (Fig. 4G). As expected, sulfate addition increased the estimated flux through sulfate reduction from 0 to 3.09 6 0.03 mmol (P , 0.0001). The estimated flux of carbon and electrons through D. vulgaris was poorly constrained due to redundant pathways. However, the flux did increase within related possible flux distributions. For example, the maximum flux of hydrogen production by D. vulgaris increased from 3.16 6 0.02 to 4.06 6 0.02 mmol with the addition of sulfate, which agrees with the increased relative abundance of the D. vulgaris proteome in the community metaproteome (Fig. 1B). M. hungatei and M. concilii were represented by hydrogenotrophic methanogenesis and acetoclastic methanogenesis, respectively. Compared to the control condition, the addition of sulfate halted hydrogenotrophic methanogenesis in M. hungatei and increased the flux through acetoclastic methanogenesis in M. concilii by 1.42-fold (P , 0.0001).
The model estimated the cross-feeding fluxes from R. cellulolyticum to the other microorganisms (Fig. 5C). The fluxes of acetate and CO 2 from R. cellulolyticum to M. concilii and M. hungatei were each 0.97 mmol in the control condition. In the sulfate addition condition, the flux of H 2 consumption by M. hungatei fully stopped, suggesting that substrate competition for H 2 by D. vulgaris was the main cause of the lower methane yields in the sulfate-amended quadculture. Indeed, the H 2 consumption flux of D. vulgaris increased from 0 in the control condition to 3.09 mmol under the sulfate addition condition. For comparison, sulfate addition increased H 2 production by R. cellulolyticum from 3.16 to 4.63 mmol. Overall, this analysis demonstrates that sulfate addition deeply restructures carbon fluxes and metabolic handoffs through the four-species synthetic community, ultimately affecting overall community metabolism.

DISCUSSION
Synthetic communities have been shown to be a powerful approach to studying the interactions among major microbial guilds in the degradation of cellulosic biomass to methane and CO 2 . A bi-culture of D. vulgaris and Methanococcus maripaludis showed the requirement of hydrogen transfer for their syntrophic growth coupled with methane production (20). Cellulose degradation by Clostridium cellulovorans was enhanced by coculturing with Methanosarcina barkeri or Methanosarcina mazei (21). Coculturing R. cellulolyticum with D. vulgaris was also shown to accelerate cellulose degradation (22). However, these bi-culture can only reveal binary interactions between two organisms. By using a four-species synthetic community and comparing that to its constitutive parts, we investigated tertiary interactions (i.e., interactions of interactions) among four microorganisms involved in methanogenesis from cellulose. A tertiary interaction has positive synergy if the effect of combining three or more organisms on a phenotype of the community is greater than the addictive effects of combining two of these organisms separately. In contrast, a tertiary interaction has negative synergy when the effect of combining three or more organisms on a phenotype is less than the sum of the separate effects of combining the subsets. Here, we demonstrated that the bi-culture of R. cellulolyticum with any of the other 3 species increased cellulose degradation relative to the R. cellulolyticum mono-culture. This provides a baseline to determine the type of tertiary interactions. The quad-cultures further improved the cellulose degradation over any of the bi-cultures, but the increase of cellulose degradation by the quad-culture over the mono-culture (25.63 mg) was less than the sum of the increases by the 3 bi-cultures over the mono-culture (13.38 mg, 13.78 mg, and 19.34 mg). Such a negative synergy suggests that the putative beneficial effects of pH stabilization and cellulose colonization (21) provided by the methanogens and D. vulgaris may be antagonistic for cellulose degradation by R. cellulolyticum (23). The methane production by the tri-culture of R. cellulolyticum, M. concilii, and M. hungatei (0.62 6 0.02 mmol) was approximately equal to the combined methane production by R. cellulolyticum-M. concilii bi-culture (0.18 6 0.003 mmol) and R. cellulolyticum-M. hungatei bi-culture (0.42 6 0.007 mmol), while the methane production of the quad-culture (1.37 6 0.03 mmol) was higher than the combined methane production by the R. cellulolyticum, M. hungatei, and D. vulgaris tri-culture (0.89 6 0.02 mmol) and the R. cellulolyticum, M. concilii, and D. vulgaris tri-culture (0.16 6 0.001 mmol), hence indicating a positive synergy. Mechanistically, it may be that D. vulgaris enhances methanogenesis by producing the substrates CO 2 , acetate, and H 2 and removing lactate (24,25). The two methanogens can enhance each other's methanogenesis by also consuming the substrates. The whole of these positively synergistic effects was greater than the sum of the individual effects on methanogenesis.
Synthetic communities can be modeled to estimate their intracellular metabolic fluxes and cross-feeding fluxes, which can be used to infer competitive and cooperative interactions among the organisms. Genome-scale flux balance analysis (FBA) models of individual organisms have been combined to construct FBA models of their cocultures (25)(26)(27). However, genome-scale FBA is vastly underparameterized and requires strong assumptions on the objective function for optimization of the fluxes in an organism (28), This is further complicated in a synthetic community with different organisms possibly requiring divergent objective functions to maximize the total fitness of the community. In addition, reliable FBA models require labor-intensive curation of the genome-scale metabolic networks of all members of a synthetic community. As a result, the previous studies focused on the synthetic communities of microorganisms with preexisting FBA models (29). Alternatively, a kinetic model of anaerobic digestion, ADM1 (30), models the major metabolic processes in the degradation of complex biomass (e.g., cellulose) to methane and CO 2 . However, the lack of knowledge of kinetic parameters in a study would limit the prediction accuracy of the ADM1 modeling (31). Here, we constructed a simple stochiometric model for synthetic communities. Similar to ADM1, the stochiometric modeling approach considers only the major metabolic processes of the organisms, but it enables quantifying metabolic fluxes without knowledge of substrate concentrations or kinetic parameters. Unlike a genome-scale FBA model, our simple stoichiometric model can be fitted just using the experimental data without assuming any objective function (32). These complementary modeling approaches could be combined to build a multiscale model for synthetic communities. Specifically, the organism-level boxes in a communitylevel stoichiometric model (Fig. 5) can be expanded to individual genome-scale FBA models for selected organisms of interest. The estimated fluxes by the multiscale flux balancing analyses can then be integrated with a parameterized kinetic model for each of the organisms to predict synthetic community dynamics and its responses to disturbances.
The stoichiometric model and metaproteomics analysis characterized the metabolic activities of microorganisms in two different ways, so we compared the fluxes through the modeled reactions with the protein abundances of marker enzymes in those reactions (see Table S3 in the supplemental material). The stoichiometric model showed that sulfate addition zeroed out the flux through hydrogenotrophic methanogenesis. This change was also reflected in the significantly decreased abundance of M. hungatei formate dehydrogenase, a marker enzyme for hydrogenotrophic methanogenesis. However, sulfate addition significantly increased the flux of acetoclastic methanogenesis in M. concilii, while the relative abundance of M. concilii methyl-coenzyme M reductase significantly decreased. In addition to the protein abundance of the corresponding enzyme, a reaction may also be regulated by a variety of other mechanisms, such as posttranslational modifications and allosteric regulations (17). Compared with the control condition, the accumulation of CH 4 was reduced by only 12% in the sulfate addition condition, while the flux of hydrogenotrophic methanogenesis decreased to 0 mmol, so a 1.42-fold increase in acetoclastic methanogenesis flux should be reasonable.
The quad-culture and its model were used to study the effects of increased sulfate availability on the microbial communities. Such a perturbation may occur in estuarine wetland soil communities as a result of more prevalent sulfate-laden seawater intrusion from rising seawater levels. The impact of seawater intrusion has been studied using natural wetland ecosystems (33,34), demonstrating that seawater intrusion reduced the CH 4 and CO 2 fluxes from the wetland. However, it was not clear how microbial interactions responded to seawater intrusion because of the great complexity of natural communities. In this study, we showed that sulfate addition reduced CH 4 and CO 2 emissions, consistent with results obtained with natural ecosystems (33,35). The main reason for decreased CH 4 production was the sharp decrease in H 2 available for the hydrogenotrophic methanogen. The competition for H 2 use between methanogens and sulfate reducers has been proposed to have a thermodynamic basis, as hydrogenotrophic sulfate reduction is more energetically favorable than hydrogenotrophic methanogenesis (36). Therefore, such competitions exist in reactor systems and environments with transient sulfate (37,38). In addition, possible interactions can be estimated by modeling what is consumed and produced by the individual species. Our result demonstrated that sulfate addition increased both the production of acetate and lactate by R. cellulolyticum and their consumption by M. concilii and D. vulgaris. This suggested that sulfate addition intensified the carbon transfer from R. cellulolyticum to M. concilii and D. vulgaris. Although the flux of H 2 production by R. cellulolyticum increased, the flux of H 2 consumed to produce H 2 S in D. vulgaris increased even more. This, in turn, further intensified the competitive interactions between M. hungatei and D. vulgaris for H 2 as the substrate in the presence of sulfate.
Many syntrophic communities have been shown to develop around the exchange of H 2 and organic acids (9,39). In this study, our stoichiometric analysis was unable to constrain the possible exchange of lactate or H 2 between R. cellulolyticum and D. vulgaris, even in the absence of sulfate (Table S3). However, limited changes in hydrogenase abundance were observed in both R. cellulolyticum and D. vulgaris, indicating a likely exchange of lactate or some other oxidizable carbon intermediate. While the oxidation of lactate to acetate and production of H 2 is expected to be unfavorable under the H 2 partial pressures studied, alternative organic compounds could be intermediates, including pyruvate or alanine, which have been demonstrated to produce H 2 in sulfatereducing bacterium at high H 2 partial pressures (9). Furthermore, D. vulgaris was observed in both the sulfate-replete and -deplete conditions, indicating growth under both conditions and, therefore, at least a partial exchange of lactate or another metabolite to support D. vulgaris growth. In the presence of sulfate, the lactate dehydrogenase in R. cellulolyticum also increased, indicating a possible increased participation of the lactate exchange in the flux distributions (Table S3). While D. vulgaris should outcompete M. hungatei for H 2 in the sulfate addition condition, the relative total protein abundance of M. hungatei remained at 3% in both conditions. This could be explained if D. vulgaris had a long lag phase in the quad-culture and M. hungatei entered the stationary phase before being outcompeted by D. vulgaris. Alternatively, M. hungatei might employ direct electron transfer using its electrogenic pili (40) or utilize other compounds, such as formate (41), as the electron donor.
We anticipate that this four-species community, as well as the future construction of additional multispecies synthetic communities, introduces the level of complexity needed to more realistically model carbon cycling in natural anaerobic communities. By combining metabolic profiling, metaproteomics, and modeling, we observed many community behaviors as designed, especially the adaptive cross-feeding fluxes and the shifting substrate competitions under changing environments. More interestingly, the multispecies cultures also revealed some emergent properties of microbial communities, including synergistic epistasis and antagonistic epistasis of tertiary interactions. By quantifying the expected competitive/cooperative behaviors and revealing the emergent interaction properties, the synthetic communities may pave the way toward more accurate modeling of natural communities.
Cultures, growth conditions, and experimental design. Four strains were revived by respective media and conditions commonly published in previous studies. Briefly, modified VM medium supplemented with 2.0g/L yeast extract and LS4D medium were used for R. cellulolyticum H10 and Desulfovibrio vulgaris Hildenborough, respectively, and both were cultured at 34°C (43,44). M. hungatei JF1 and M. concilii were grown in RST medium with H 2 /CO 2 (80:20) in the headspace and 50 mM acetate at 37°C (45). M. concilii required about 5 weeks before noticeable growth and acclimation of acetate occurred. Subsequently, active inocula (optical density at 600 nm [OD 600 ] = 0.5) were transferred to the modified VM medium with 10 g/L cellulose as the sole carbon source for experiments (44).
All fermentation experiments were carried out in 160 mL serum bottles with 20 mL of modified VM medium containing 10 g/L cellulose; the chamber gas contained 3% H 2 and 97% N 2 . The modified VM medium contained 1.25 mg/L of FeSO 4 Á6H 2 O, which introduced approximately 0.14 mmol sulfate to the control condition. Revived cells were counted by flow cytometry; the number of inoculated cells per species was about 5 Â 10 7 for all mixed cultures. All cultures were incubated at 34°C with shaking at 200 rpm for 7 days and sampled each day.
During the sampling, 2 mL of well-shaken culture liquid was sampled, centrifuged at 14,000 rpm for 15 min to separate the cells and cellulose from the supernatant, and stored at 280°C for later analysis. We collected 5 mL gas phase products in vacuumed Labco exetainer glass vials (catalog no. 837w) for later analysis.
After sampling, 2 mL fresh medium with 10 g/L cellulose was added to the cultures, pH was adjusted to 7.3, and produced gases were vented to restore serum bottle headspace to atmospheric pressure. To mimic seawater intrusion, sulfate-treated experiments received 2 mL of 100 mM K 2 SO 4 solution every day from the first day of cultivation, and the control group received the same volume of sterile water.
Chemical analysis. Cellulose in the medium was measured by the phenol-sulfuric acid method as described previously (46). Headspace pressure was measured daily with a pressure gauge (Cole Parmer; catalog no. SK-68900-24) before sampling liquid and gas phases. Gas composition (H 2 , CO 2 , and CH 4 ) was analyzed from the stored Labco vacuum containers with a gas chromatograph (GC-8A; Shimadzu, Japan). The absolute gas composition in each bottle was calculated with the ideal gas law. The remaining sulfate in the medium was measured by ion chromatography (IC) and used to estimate H 2 S production in the cultures. The fermentation products (lactate and acetate) in liquid phase and soluble sugars (cellobiose and glucose) were analyzed by high-performance liquid chromatography (HPLC) as reported previously (47).
Protein identification and quantification. Samples from quad-cultures were selected for global metaproteomic analysis. Protein extraction was carried out according to a previously reported method (48,49). Briefly, harvested cell pellets were washed twice with Nanopure water, cells were suspended in lysis buffer (containing 10 mM Tris-HCl, 1% SDS, and 0.1 M dithiothreitol [DTT]), and incubated at 60°C for 1 h, and then the Tertiary Interactions in a Four-Species Community mBio supernatant was collected after centrifugation. Proteins were then precipitated by trichloroacetic acid overnight at 4°C and pelleted by centrifugation. Protein pellets were washed with ice-cold acetone three times and resuspended in guanidine buffer. Protein concentrations were quantified by bicinchoninic acid assay (50). Twenty-milligram samples of proteins were processed using filter-aided sample preparation and digested using trypsin-LysC mixture. Each sample was analyzed using two-dimensional liquid chromatography-tandem mass spectrometry (2D-LC-MS/MS) on an Orbitrap Fusion Tribrid mass spectrometer (Thermo Fisher Scientific, USA) at the IDeA National Resource for Quantitative Proteomics. Tryptic peptides were separated into 46 fractions on a 100-by 1.0-mm Acquity BEH C 18 column (Waters) with a 50-min gradient from 99:1 to 60:40 buffer ratio under basic pH conditions and then consolidated into 8 superfractions. Each superfraction was then separated by reverse phase XSelect CSH C 18 2.5-mm resin (Waters) on an in-line 120-by 0.075-mm column using an UltiMate 3000 RSLCnano system (Thermo Fisher Scientific, USA). Peptides were eluted using a 60-min gradient from 98:2 to 65:35 buffer ratio. Eluted peptides were ionized by electrospray (2.4 kV) followed by mass spectrometric analysis on an Orbitrap Fusion Tribrid mass spectrometer (Thermo Fisher Scientific, USA). MS data were acquired using a Fourier transform MS (FTMS) analyzer in profile mode at a resolution of 240,000 over a range of 375 to 1,500 m/z. Following high-energy collisional dissociation (HCD) activation, MS/MS data were acquired using the ion trap analyzer in centroid mode and normal mass range with normalized collision energy of 28 to 31% depending on charge state and precursor selection range. Mass spectrometry spectra were searched using Sipros Ensemble against a matched protein database constructed from the genomes of the four microorganisms making up the synthetic community (51,52). Raw search results were filtered to achieve a 1% false-discovery rate (FDR) at the peptide level, estimated by the target-decoy approach. Peptide identifications are assigned to proteins or protein groups following the parsimonious rule (51). A protein group included multiple proteins sharing a set of common peptides, at least one of which cannot be attributed to any other protein or protein group. A minimum of one unique peptide was required for each identified protein or protein group. To avoid ambiguity in data analysis, only proteins were used for biological analysis, while protein groups were not evaluated. Intensity-based labelfree quantification for protein was performed using ProRata (53,54). The protein abundance was represented by the total peak height of all quantified unique peptides from a protein (17).
Construction of the stoichiometric model. A set of overall reactions representing the documented metabolisms of the individual species were fit into the daily measured amounts of metabolites in the constructed synthetic communities (Tables S1 and S2). The metabolism of R. cellulolyticum was modeled to consume glucose to eliminate artifacts from cellulase production. The overall metabolism of R. cellulolyticum was modeled as lactate fermentation, hydrogenic acetogenesis. The low observed ethanol in ethanol/acetate fermentations, accounting for less than 0.36% of the carbon flux, led us to remove this pathway from analysis. D. vulgaris was modeled to oxidize lactate to acetate, coupled with reduction of protons or sulfate, or oxidize H 2 coupled to sulfate reduction. The methanogens were represented by hydrogenotrophic methanogenesis and acetoclastic methanogenesis. The code used for the fitting can be found at GitHub (https://github.com/thepanlab/Community_stoichiometric _model). The resulting fits indicate the relative contribution of each microorganism to the overall observed transformations and provide a basis for quantifying metabolic contributions. Deviations between the inferred metabolic transformations and the measured transformations indicate uncertainty in our understanding of the metabolic transformations or limitations in analytical capacity.
Statistical analysis. Metabolic measurement data were compared between different conditions using Student's t test. The differences in protein abundance in metaproteome were analyzed by the DeSeq R package (55), and the P values were adjusted to q values using the Benjamini-Hochberg method for multiplecomparison correction (55). Proteins with q values of ,0.05 were regarded as statistically significant.
Data availability. Proteomic data are available at the ProteomeXchange Consortium via the Proteomics Identification Database (PRIDE) partner repository with the data set identifier PXD035759.

SUPPLEMENTAL MATERIAL
Supplemental material is available online only.