Pseudomonas spp. are key players in agricultural biogas substrate degradation

Anaerobic degradation (AD) of heterogeneous agricultural substrates is a complex process involving a diverse microbial community. While microbial community composition of a variety of biogas plants (BPs) is well described, little is known about metabolic processes and microbial interaction patterns. Here, we analyzed 16 large-scale BPs using metaproteomics. All metabolic steps of AD were observed in the metaproteome, and multivariate analyses indicated that they were shaped by temperature, pH, volatile fatty acid content and substrate types. Biogas plants could be subdivided into hydrogenotrophic, acetoclastic or a mixture of both methanogenic pathways based on their process parameters, taxonomic and functional metaproteome. Network analyses showed large differences in metabolic and microbial interaction patterns. Both, number of interactions and interaction partners were highly dependent on the prevalent methanogenic pathway for most species. Nevertheless, we observed a highly conserved metabolism of different abundant Pseudomonas spp. for all BPs indicating a key role during AD in carbohydrate hydrolysis irrespectively of variabilities in substrate input and process parameters. Thus, Pseudomonas spp. are of high importance for robust and versatile AD food webs, which highlight a large variety of downstream metabolic processes for their respective methanogenic pathways.

Methodologies for analyzing and interpreting omics data have rapidly changed in the last decade. Tools for network analyses such as MENA 20 , SparCC 21 or CoNet 22 are frequently used to predict interactions between microorganisms. Such network calculations are mainly based on 16S rRNA gene amplicon abundance data and therefore anaerobic degradation-based findings have to be interpreted with caution, as metabolic functions are difficult to predict. Other tools like STRING 23 focus on discovery of protein-protein interaction networks for explanation of microbial interactions. Unfortunately, the relatively low number of reliable database entries limits those tools for AD.
In contrast to most other studies, we analyzed not only the metaproteome of a single BP but sixteen large-scale BPs. Furthermore, we measured the metaproteome of five independent replicates for each BP (same time point) to produce robust results. Main goals of this study are (i) to identify most important parameters driving the AD on protein level, (ii) to group the BPs according to their metaproteome (taxonomic and functional profiles), (iii) to arrange the BPs corresponding to their prevalent methanogenesis pathway and (iv) to identify microbial key players and their interaction patterns by a metabolic and microbial network analysis. The overall aim was to gain a better understanding of the metabolic processes during AD with a focus on methanogenesis as well as explore possibilities of metaproteomics for practical applications.

Results
Methanogenic proteins are dominating in most biogas plants. In total, 5,854 protein groups from 2,178 different species were identified in all BPs, whereof between 3,238 (BP01) to 3,624 (BP11) were found in average in each BP. About 77% of the identified protein groups were affiliated to Bacteria and 23% to Archaea (for community composition on different levels see Supplementary Table S1), which were assigned to 4,894 protein groups (about 84%) and their corresponding molecular function(s) (K numbers) according to KEGG [24][25][26] (Kyoto Encyclopedia of Genes and Genomes).

Biogas plants are grouped by taxonomy and function.
Principal component analysis (PCA) on both, protein and species level revealed nine very similar clusters, while on functional level only six clusters were obtained ( Fig. 2 and Supplementary Fig. S1). The nine different clusters (C1-9) on species level (Fig. 2B) were also found by microbial community composition analysis (Supplementary Table S1).
Bacterial families Clostridiaceae, Peptococcaceae and Pseudomonadaceae were prominent in all BPs ( Supplementary Fig. S3). C1 (BP06 to BP08 and BP12), C2 (BP15) and C3 (BP10, BP13) additionally showed higher abundances of Thermoanaerobacteraceae and Thermoanaerobacterales Family III. C4 (BP01) was comprised with higher abundances of Peptostreptococcaceae, while high relative abundances for members of Petrotogaceae were observed in C5 (BP16). Proteins from family Porphyromonadaceae were particularly  www.nature.com/scientificreports www.nature.com/scientificreports/ HyMe also showed highest number of modules (42), followed by BoMe (32) and AcMe (26). Based on their topological role, nodes acting as generalists were observed for all three networks. Most of these generalistic nodes were assigned to Firmicutes, irrespective of predominant methanogenic pathway. In addition, higher proportions of generalistic nodes were assigned to Proteobacteria and Actinobacteria (HyMe and BoMe) or Euryarchaeota (AcMe). Only two nodes were acting as a generalist in more than one network. S0024 (Peptococcaceae bacterium 1109) was a module hub in AcMe and HyMe, while S0271 (Defluviitoga tunisiensis) was a module hub in  Table S1). Green color of edges correspond to positive interactions.    Steps of hydrogenotrophic methanogenic pathway are colored in blue, steps of acetoclastic methanogenic pathway in red and steps of methylotrophic methanogenesis in black. Additionally, steps being part of hydrogenotrophic and acetoclastic methanogenesis are colored in grey, while steps being part of all three pathways are colored in yellow. Bar  www.nature.com/scientificreports www.nature.com/scientificreports/ Pseudomonas spp. carbohydrate metabolism was independent on methanogenic pathway.
Further analyses of those shared nodes and edges identified, among others, 15 different Pseudomonas spp. acting as key players in AD for all main methanogenic pathways. They were found to share 37 edges among themselves in all networks (Fig. 3). Functional analysis revealed that they were mainly active during hydrolysis of carbohydrates using the Entner-Doudoroff (ED) pathway ( Supplementary Fig. S5). Most of the high abundant proteins (e.g. succinate CoA ligase, enolase or ATP synthase) correspond to P. fluorescens. Expression of most protein-coding genes was consistent among all networks (Supplementary Fig. S6). Proteins upregulated in HyMe and BoMe were C3K6H6 (enolase) and A0A010SIZ3 (L-glaceraldehyde-3-phosphate reductase). In contrast, C3K613 (fatty acid oxidation complex subunit alpha) was upregulated in BoMe and AcMe, while A0A085VRFO (ornithine decarboxylase), A0A176V720 (glyceraldehyde-3-phopshate dehydrogenase) and A0A098T689 (GMP-synthase) were upregulated only in AcMe.
Methanogenic protein patterns allow assignment to one or more methanogenic pathways. Proteins of all methanogenic pathways were found in each BP but in very different relative abundances ( Fig. 4). High protein abundances were observed for Methanoculleus bourgensis (especially in HyMe) and M. barkeri (most AcMe) in all BPs. In addition, BP05 and BP13 showed high abundances of MCR from Methanothrix soehngenii ( Table 1). Most abundant proteins specific for HyMe (Fig. 4) were methylenetetrahydromethanopterin dehydrogenase, mainly from M. bourgensis, Methanospirillum hungatei and Methanoculleus marisnigri, and 5,10-methylenetetrahydromethanopterin reductase (5,10-Methylene-THMPT) from M. bourgensis and Methanosphaerula palustris (Table 1). In contrast, AcMe was governed by Acetyl-coenzyme A synthetase, mainly from M. soehngenii, M. mazei and in lower abundances from M. acetivorans. In addition, subunits of acetyl-CoA decarbonylase/synthase complex were abundant, which were mainly affiliated to M. soehngenii and M. mazei. Higher abundances for those proteins were also observed for SAOB S. schinkii (Table 1). Co-methyltransferase was the only representative protein for methylotrophic methanogenesis with high abundances in few BPs (BP02, BP14 and BP09), which was mainly affiliated to M. mazei and Methanosarcina sp. Ant1.
Link between process parameters, taxonomy and microbial functions. Process temperature was the most important driver for the metaproteomic profiles ( Table 2). Members of the family Porphyromonadaceae were negatively correlated to temperature and families Synergistaceae and Petrotogaceae were positively correlated (Supplementary Table S4). Correlation of Petrotogaceae was mainly linked to thermophilic D. tunisiensis and its Glyceraldehyde-3-phosphatase dehydrogenase, as well as pyruvate-phosphate dikinase. No correlation of temperature and Archaea was observed.
Similarly, bacterial proteins (carbohydrate metabolism and transporter proteins) of Peptococcaceae, Ruminococcaceae and Tissierellaceae were positively correlated with pH, but no archaeal proteins (Supplementary  Table S4). Proteins from family Methanomicrobiaceae (mainly M. bourgensis) were positively and from Methanosaetaceae and Methanocellaceae (mainly MCR subunits of M. soehngenii, Methanosaeta harundinacea and M. mazei) were negatively correlated to VFA concentrations.

Discussion
We observed for almost all metabolic steps of AD high protein abundances for transport, glycolysis/gluconeogenesis and methanogenesis. Proteins for all different methanogenic pathways were present in each BP, indicating that a mixture of methanogenic pathways simultaneously convert agricultural substrates to biogas in large-scale BPs. Nevertheless, majority of BPs was dominated by either hydrogenotrophic or acetoclastic methanogenesis as prevalent pathway. In comparison with other approaches like stable isotope probing combined with a nucleic acid approach (DNA-SIP or RNA-SIP), metaproteomics enable accurate access to microbial phylogeny, function and its abundances in any scale of the bioreactors. In contrast, RNA-or DNA-SIP approaches are well established for small-scale bioreactors [27][28][29][30] , but expensive and artificial as results cannot easily be upscaled. Robustness and reproducibility of our approach is supported by the results for the five replicates of a BP: in only two cases one replicate was assigned to another main methanogenic pathway compared to remaining four replicates. This approach could also be attractive for plant operators. If they know the main methanogenic pathway, they could adapt relevant process parameters to optimize the AD processes in their digester to increase biogas production. For instance pH adjustment is biogas rate limiting, if hydrogenotrophic methanogenesis prevails 31 . Nevertheless, for further critical review of our approach, more metaproteome data sets and its associated process parameters from large-scale BPs have to be analyzed.
Numbers and compositions of BP clusters differed on taxonomic and functional level (Fig. 2), suggesting a highly adaptive metabolic network of various active key members in a complex microbial community. Even if similar microbial communities were present, the members fulfilled different ecological processes resulting in defined interaction patters. This in turn suggests a high grade of specialization, which is dependent on a variety of factors, such as process parameters, presence and absence of potential interaction partners or bioavailability of substrates. Therefore, to gain a preferably comprehensive understanding of the AD processes, the need for holistic approaches like metaproteomics are profitable. Results of microbial community composition on 16S rRNA gene for instance has the risk of underestimating the importance and activity of especially methanogens compared to proteomics based approaches 10 .
Environmental variables, such as substrate or temperature, are known to affect community composition during AD as revealed by nucleic acid based analyses 6,10,32 . Similarly, a variety of parameters (e.g. temperature, pH, feedstock, VFA) were found to significantly affect the metaproteome profiles of BPs what is in accordance with other studies 5,33 . Some already suggested suitable marker organisms for different types of BPs, as well as potential biomarkers for AD in BPs could be approved 5 . In addition, some potential new biomarkers were described here. In our study is D. tunisiensis a promising candidate for a marker organism of thermophilic BPs. This species is www.nature.com/scientificreports www.nature.com/scientificreports/ known to be a key player for hydrolysis in BPs, and its genome encode a variety of genes associated with complex polysaccharide degradation 34 . Positive correlations of glycolytic proteins for D. tunisiensis indicate a high metabolic activity of Petrotogaceae during degradation of complex carbohydrates under thermophilic conditions, as observed in other studies 35 .
Peptococcaceae bacterium 1109 was highly abundant in most BPs and positively correlated to pH. Therefore, a decrease in pH (as observed during acidification) could possibly be detected by a decrease of proteins from this species. In addition, observed correlations of pH and different species were solely positive. This is surprising as all BPs were single-stage fermenters and negative correlations between hydrolyzing organisms and pH were expected, due to their lower pH optima. This findings indicate a highly adapted bacterial community, which is able to hydrolyze substrates efficiently even at high pH ≥ 7.7 (see Supplementary Table S2) and possibly overcome bottlenecks during rate-limiting hydrolysis. Many different pathways seem to be influenced by VFA concentration (Supplementary  Table S4), which is in line with previous results as high VFA concentrations inhibit both, hydrolysis and methanogenesis 36 . Nevertheless, it has been reported that high acetate (i.e. about 2400 mg/L) and butyrate concentrations (about 1,800 mg/L) seem to have no effect on methanogenesis, while high propionate concentrations (900 mg/L) are more critical 37 . Higher proportions of propionate have been observed for three BPs (BP 10, BP 12 and BP 16), that were assigned to HyMe. Although there are more hydrogenotrophic BPs, it is noticeable that these three BPs showed lower abundances of known SAOB, such as T. phaeum or S. schinkii (Supplementary Table S1). Acetate oxidation by syntrophs can be a rate limiting step of hydrogenotrophic methanogenesis 38 , and therefore a lack of SAOB during hydrogenotrophic methanogenesis is likely to lead to an accumulation of VFAs, which in turn can inhibit methanogenic activities. Therefore, monitoring the abundance of proteins from SAOBs seems to be a suitable way to check the stability and performance of hydrogenotrophic BPs. In addition, our results suggest that MCR from Methanosarcinales is as promising biomarker candidate for acidification. Acetoclastic organisms in combination with the presence of mixed acid fermentation enabled a fast metabolization of different VFAs. Acetoclastic microorganisms have to metabolize more substrate to obtain the same energy as hydrogenotrophic microorganisms 39,40 . As differences in process parameters as well as used substrates could also strongly influence VFA concentrations 41,42 , future studies should evaluate whether VFA monitoring is useful for all types of BPs.
Network analyses approach revealed distinct metabolic and microbial interaction patterns for each methanogenic pathway. Different number of nodes and edges for each network indicated highly complex microbial interactions patterns, with HyMe being the most complex one. In contrast, BoNe showed the fewest number of nodes and edges and could be considered to be the least complex network. Possibly, the prevalence of two methanogenesis pathways prevent the organisms from a too deep specialization and consequently, less organisms and interactions are necessary for a stable AD process. In addition, network analyses support the importance of Defluviitoga tunisiensis, as well as P. bacterium 1109 for AD. As no common interactions of both species with other microorganisms could be observed among all networks, it can be assumed that their high flexibility for different interaction partners, may lead to their generalistic behavior. In contrast, highly specialized organisms can substantially contribute to AD such as 15 different Pseudomonas spp., which shared 37 edges among themselves over all networks and BPs. Stable expression for most proteins of the different Pseudomonas spp. (Supplementary Fig. S6), indicating a highly conserved carbohydrate metabolism of those species, even if process parameters are different ( Supplementary Fig. S2). Additionally, Pseudomonas spp. seems to prefer the Entner-Doudoroff (ED) pathway instead of the Embden-Meyerhof (EM) pathway ( Supplementary Fig. S5) for glucose degradation, even if ED is energetically unfavorable. This finding may be linked to the lack of phosphofructokinase for most Pseudomonas spp., which is the key enzyme of EM 43 . Different Pseudomonas spp. metabolize glucose through ED which have been previously shown in lab based experiments of glucose metabolism 44,45 . Results of this study can be used to optimize process conditions to facilitate metabolic activity of Pseudomonas spp. in agricultural BPs. Even if no direct interactions between SAOBs and hydrogenotrophic methanogens could be observed, the numbers of interactions of SAOBs in each of the three networks clearly indicate the great importance of SAOBs for HyMe, while they seem to play a minor role during AcMe. During AcMe a big proportion of produced acetate is metabolized by acetoclastic methanogens, and therefore SAOBs have to compete with them. These findings were supported by higher mean protein abundances for T. phaeum (AcMe: 0.2%, HyMe: 1.0%, BoMe: 0.6%), as well as S. schinkii (AcMe: 0.3%, HyMe: 4.7%, BoMe: 1.8%).

conclusions
Metaproteome analyses of 16 agricultural large-scale BPs enable a deeper insight into AD. Temperature, pH, substrate and VFA concentration were identified as main drivers for metaproteomic profiles. Comprehensive correlation analyses enabled the identification of potential marker organisms for defined process conditions, such as Petrotogaceae for high temperatures. Moreover, monitoring the MCR of Methanosarcinales could be a suitable biomarker to recognize ongoing acidification and avoid process failures.
BPs clustered similar on both, species and protein level but differed on functional level, indicating a high resilience and flexibility of the microbial community. BPs were classified to acetoclastic, hydrogenotrophic or a mixture of both pathways, while methylotrophic methanogenesis was of minor importance. Such classifications could be meaningful for upcoming studies.
Network analysis of each methanogenic pathway revealed that microbial interaction patterns widely differed among BPs reflecting large differences in metabolic processes in the prevalent methanogenic pathway. However, all methanogenic pathways have in common that Pseudomonas spp. were main drivers in hydrolytic processes indicating their versatile metabolism in wide-ranging process conditions and substrate variabilities. Although no interactions between hydrogenotrophic methanogens and SAOBs were observed in the networks, our data emphasize the importance of syntrophic acetate-oxidation for hydrogenotrophic methanogenesis. In addition, network analysis underline the role of D. tunisiensis and P. bacterium 1109 during AD of various agricultural substrates. Further research is required to obtain a deeper understanding of anaerobic degradation and to include this knowledge in control technology to make biogas production more flexible and efficient.

Material and Methods
Biogas plant sampling. Primary digesters of 16 large-scale BPs were sampled in Northern Bavaria in August 2016 as described previously 10 . Briefly, digester content was mixed thoroughly and all pipelines were flushed prior to sampling. Finally, approximately 300 mL of digester sludge were obtained from each BP and subsequently shock-frozen in liquid nitrogen. Samples were thereafter transferred in liquid nitrogen and stored at −80 °C for further investigations. Some process and plant parameters were collected from the plant operators, i.e. process temperature, pH, as well as types of main and additional substrates.
Determination of VfA concentration, VfA/tAc and c, n, S content. For determination of total volatile fatty acids (VFA) content and ratio of VFAs and alkalinity (VFA/TAC), samples were centrifuged for 15 min and 4,696 × g at 4 °C. 10 mL of the liquid phase were then titrated with 0.05 M sulfuric acid (Carl Roth GmbH, Karlsruhe, Germany) to pH values of 5.00, 4.40, 4.30 and 4.00. VFA as well as VFA/TAC were measured in triplicates and calculated as described elsewhere 46 . Composition of VFA was analyzed by atres Analytik (München, Germany) using an in-house gaschromatographic approach. Contents of following VFAs were measured: acetic acid, propionic acid, isobutyric acid, butyric acid, isovaleric acid, valeric acid, hexanoic acid and heptanoic acid.
For analyzing elemental composition (C, N, S) samples were dried at 50 °C for 48 h. Afterwards samples were grinded using mortar and pestle. 10 mg of each sample were oxidized in the combustion tube of the elemental analyzer Vario EL II (Elementar Analysensysteme GmbH, Langenselbold, Germany) at 1,150 °C according to manufacturer's instructions. Sulphanilic acid (Sigma-Aldrich Chemie GmbH, Taufkirchen, Germany) was used as standard. All analyses were performed in triplicate measurements. protein extraction and sample preparation for mass spectrometric analysis. Proteins were extracted according to Heyer and colleagues with modifications 33 . Shortly, 0.5 g biomass were mixed with 300 µL of a 2 M sucrose solution (Carl Roth GmbH, Karlsruhe, Germany), 800 µL of liquid phenol (Carl Roth GmbH, Karlsruhe, Germany) and 1 g silica beads with a diameter of 0.7 µm (Carl Roth). Samples were homogenized at 5 ms -1 for two minutes using a ball mill (FastPrep ® -24, MP Biomedicals, CA, USA), followed by centrifugation for 10 min at 10,000 × g for phase separation. 333 µL of the upper phenol phase, containing the proteins, were transferred and mixed with the same volume of a 1 M sucrose solution. After incubation for 10 min at 1,000 rpm, samples were centrifuged for complete phase separation using same conditions as mentioned above. 150 µL of the upper phenol phase were then mixed with fivefold volume of 0.1 M ammonium acetate (Merck KGaA, Darmstadt, Germany) solution in methanol (Carl Roth) and stored for one hour at −20 °C for protein precipitation. After precipitation, samples were centrifuged for 20 min at 10,000 × g and 4 °C. Protein pellets were washed twice with 500 µL of 80% acetone (Carl Roth) and 70% ethanol (Carl Roth). Each washing step was followed by incubation at −20 °C for 15 min and subsequent centrifugation for 10 min at 12,000 × g and 4 °C. Pellets were resuspended in 500 µL of a resuspension solution containing 7 M urea (Merck KGaA, Darmstadt, Germany), 2 M thiourea (Carl Roth) and 0.01 g mL −1 dithiothreitol (Carl Roth). Five independent replicates were taken from each of the 16 BPs, resulting in a total number of 80 metaproteome measurements.
Protein concentrations were determined using amido black assay modified from Schweikl et al. 47 . Briefly, 150 µL of amido black dye solution (0.26 mg mL −1 dissolved in 10% acetic acid and 90% methanol (all from Carl Roth)) was added to 25 µL of protein solution. After 5 min of incubation at room temperature, samples were centrifuged for 10 min at 16,000 × g. Supernatant was discarded, and 250 µL of washing solution, containing 10% acetic acid and 90% methanol, were added. Samples were centrifuged as mentioned above. This step was repeated twice. Protein pellet was thereafter resuspended in 175 µL of 0.1 M sodium hydroxide solution (Carl Roth). Protein concentration was determined by measuring absorbance at 615 nm using µDrop-plate (Thermo Fisher Scientific, Waltham, USA) and bovine serum albumin (Carl Roth) as calibration standard according to manufacturer's instructions.
For each sample, the volume containing 50 µg of protein were mixed with fivefold the volume of absolute acetone and incubated for one hour at −20 °C. After incubation, samples were centrifuged for 20 min at 12,000 × g at room temperature and supernatant was discarded. Remaining protein pellets were then dried for one hour at 35 °C using a rotary evaporator. Pellets were resuspended in denaturing buffer (7 M urea, 2 M thiourea and 0.01 g/ mL dithiothreitol) and mixed with Roti ® -Load 1 (Carl Roth). Polyacrylamide gel electrophoresis was carried out with a 12% acrylamide separating gel and a 4% stacking gel for seven minutes at 200 V. Protein bands were visualized by using Rotiphorese ® Blue R (Carl Roth). After destaining, gels were scanned with Molecular Imager ® Gel Doc ™ XR + (Bio-Rad, Hercules, USA).
All proteins were excised as one band with a length of approximately 10 mm from the stained gel lanes. Each piece was cutted into six to eight smaller slices. All slices of a sample were transferred to a new 1.5 mL reaction tube for subsequent tryptic digestion modified from Shevchenko et al. 48 . Peptide lysates were completely dried for about 90 min at 35 °C. Samples were desalted using C18-ZipTips (Sigma-Aldrich Chemie GmbH, Taufkirchen, Germany) following the manufacturer instructions.
Mass spectrometric analysis was performed on a Q Exactive HF mass spectrometer (Thermo Fisher Scientific, Waltham, MA, USA) with a TriVersa NanoMate (Advion, Ltd., Harlow, UK) source in LC-chip coupling mode as described elsewhere 49 . Briefly, the mass spectrometer was set on loop count of 20 using for MS/MS scans with higher energy collision dissociation (HCD) at normalized collision energy of 30%. MS scans were measured at a resolution of 120,000 in the scan range of 350-1,600 m/z. MS ion count target was set to 3 × 10 6 at an injection time of 80 ms. Ions for MS/MS scans were isolated in the quadrupole with an isolation window of 1.6 Da and were measured with a resolution of 15,000 in the scan range of 200-2,000 m/z. The dynamic exclusion duration was set to 30 s with a 10 ppm tolerance. Automatic gain control target was set to 2 × 10 5 with an injection time of 120 ms using the underfill ratio of 1%.
Proteome Discoverer (v2.2, Thermo Scientific) was used for protein identification and the MS/MS spectra acquired were searched with Sequest HT against the protein-coding sequences of Bacteria and Archaea of the UniProt database (release 08/2018). Enzyme specificity was selected as trypsin with up to two missed cleavages allowed, using 10 ppm peptide ion tolerance and 0.05 Da MS/MS tolerances. Oxidation at methionines as the variable modifications and carbamidomethylation at cysteines as the static modification were selected. Only peptides with a false discovery rate (FDR) < 1% calculated by Percolator 50 were considered as identified. Identified proteins were grouped by applying the strict parsimony principle, in which protein hits were reported as the minimum set that accounts for all observable peptides. Protein abundances were calculated based on the top3 approach implemented in Proteome Discoverer v2.3.

Statistical data analysis.
Proteins that were not observed in at least three out of five biological replicates of a BP were excluded from further data analyses. Principal component analysis (PCA) was done on the log2-median transformed protein intensities using the FactoMineR package 51 in R. Cluster analysis was based on the Euclidian distance between the BPs (using first, second and third dimension of PCA) using the function hclust included in the package fastcluster 52 with the ward.D2-method.
Effects of environmental variables (Supplementary Table S1) on protein profiles were evaluated using CCA using function cca of R-package vegan 53 . Statistical significance was calculated using an ANOVA-like permutation test with 9,999 permutations included in the same package. Only variables with a variance inflation factor smaller than 10 were considered for CCA 54 . Significant correlations of different taxonomic and functional levels with plant and process parameters were analyzed by spearman´s rho with a significance level of p ≤ 0.01 using R-package psych 55 .
To reduce data complexity only entries of the appropriate level with a minimum mean relative abundance of at least 0.1% were used. Differences in protein abundance levels were calculated using R-package limma 56 .

Determination of dominant methanogenic pathway. Mean values of protein abundances of each
BP were normalized to 100%. Each protein group was assigned to the corresponding KEGG orthologous (KO) group 25 . A list of KO groups being unique for different methanogenic pathways (KEGG module 00567: hydrogenotrophic methanogenesis, KEGG module 00357: acetoclastic methanogenesis and KEGG module 00356: methanogenesis from methanol) was complied. All proteins were assigned to their corresponding Enzyme Commission number (EC) 57 . Protein abundances of all KO-Terms unique for acetoclastic and hydrogenotrophic methanogenesis were summed up. A factor F was generated by dividing the relative abundance of acetoclastic proteins through the relative abundance of hydrogenotrophic proteins. All samples with F ≥ 2.5 were assumed mainly acetoclastic, while F ≤ 0.4 indicates mainly hydrogenotrophic methanogenesis. F values between those values suggest that both pathways contribute substantially to methanogenesis. network analyses. Network analyses were carried out separately for each methanogenic pathway using Cytoscape plugin CoNet 22 . To reduce data complexity, relative abundances of all proteins were summarized up to species level and a unique identifier was given for each species (Supplementary Table 2). For more valid network calculation, no mean values for each BP were calculated, but factor F was generated for each replicate. Based on their dominant methanogenic pathway, 40 replicates were used as input for the HyMe, while 25 replicates were used for AcMe calculation. Network for replicates with no clearly dominant methanogenic pathway (BoMe) was based on remaining 15 replicates. Network calculation itself was carried out with following parameters: Pearson, Kendall and Spearman correlation, as well as Bray-Curtis and Kullback-Leibler dissimilarity as methods. Minimal occurrence of observations for each set of replicates was set to at least 60% (i.e. 24, 15 and 9 observations for the different networks). Threshold was set to 2,500 top and bottom edges, so that each correlation and dissimilarity method contributed 2,500 positive and 2,500 negative edges to the initial network. MinSupport was selected to be three, so only edges supported by at least three of the five methods were kept. The method-specific p-values were computed by using the mean and standard deviation of the bootstrap distribution (100 iterations) as a parameter of the normal distribution. Method-specific p-values were then merged using the method of Brown 58 . Only edges with p < 0.05 were kept after multiple-testing correction of Benjamini and Hochberg 59 . Nodes of the final networks were assigned to modules by using GLay community algorithm 60 . Finally, within-module connectivity (z) and among-module connectivity (P i ) were calculated as described by Guimera and Amaral 61 with an automated in-house excel sheet. Peripheral nodes (specialists) were defined by z ≤ 2.5 and P i ≤ 0.62, connectors by z ≤ 2.5 and P i > 0.62, module hubs by z > 2.5 and P i ≤ 0.62 and network hubs by z > 2.5 and P i > 0.

Data Availability
The datasets generated during and/or analysed during the current study are available in the PRIDE repository, with the dataset identifier PXD014605.