Microbiome signature and diversity regulates the level of energy production under anaerobic condition

The microbiome of the anaerobic digester (AD) regulates the level of energy production. To assess the microbiome diversity and composition in different stages of anaerobic digestion, we collected 16 samples from the AD of cow dung (CD) origin. The samples were categorized into four groups (Group-I, Group-II, Group-III and Group-IV) based on the level of energy production (CH4%), and sequenced through whole metagenome sequencing (WMS). Group-I (n = 2) belonged to initial time of energy production whereas Group-II (n = 5), Group-III (n = 5), and Group-IV (n = 4) had 21–34%, 47–58% and 71–74% of CH4, respectively. The physicochemical analysis revealed that level of energy production (CH4%) had significant positive correlation with digester pH (r = 0.92, p < 0.001), O2 level (%) (r = 0.54, p < 0.05), and environmental temperature (°C) (r = 0.57, p < 0.05). The WMS data mapped to 2800 distinct bacterial, archaeal and viral genomes through PathoScope (PS) and MG-RAST (MR) analyses. We detected 768, 1421, 1819 and 1774 bacterial strains in Group-I, Group-II, Group-III and Group-IV, respectively through PS analysis which were represented by Firmicutes, Bacteroidetes, Proteobacteria, Actinobacteria, Spirochaetes and Fibrobacteres phyla (> 93.0% of the total abundances). Simultaneously, 343 archaeal strains were detected, of which 95.90% strains shared across four metagenomes. We identified 43 dominant species including 31 bacterial and 12 archaeal species in AD microbiomes, of which only archaea showed positive correlation with digester pH, CH4 concentration, pressure and temperature (Spearman correlation; r > 0.6, p < 0.01). The indicator species analysis showed that the species Methanosarcina vacuolate, Dehalococcoides mccartyi, Methanosarcina sp. Kolksee and Methanosarcina barkeri were highly specific for energy production. The correlation network analysis showed that different strains of Euryarcheota and Firmicutes phyla exhibited significant correlation (p = 0.021, Kruskal–Wallis test; with a cutoff of 1.0) with the highest level (74.1%) of energy production (Group-IV). In addition, top CH4 producing microbiomes showed increased genomic functional activities related to one carbon and biotin metabolism, oxidative stress, proteolytic pathways, membrane-type-1-matrix-metalloproteinase (MT1-MMP) pericellular network, acetyl-CoA production, motility and chemotaxis. Importantly, the physicochemical properties of the AD including pH, CH4 concentration (%), pressure, temperature and environmental temperature were found to be positively correlated with these genomic functional potentials and distribution of ARGs and metal resistance pathways (Spearman correlation; r > 0.5, p < 0.01). This study reveals distinct changes in composition and diversity of the AD microbiomes including different indicator species, and their genomic features that are highly specific for energy production.

www.nature.com/scientificreports/ volume = 375 kg and lowest input volume = 35 kg) (Table S1). Periodic increments of the raw CD resulted in increased biogas production (Fig. 1A, Table S1). This trend was observed until Day 35, after which biogas production began to decline gradually, although the CD loading was continued. The log phase of methane (CH 4 ) production started around the second day of the experiment, and reached its maximum percentage (74.1%) at Day 35 of the digestion process. It should be pointed out that after Day 35, the CH 4 percentage started to decrease, reaching 59.2% on Day 44 (Table 1). Average concentration (%) CO 2 was observed 39.52 (minimum = 27.7, maximum = 56) throughout 44 days of the digestion process. Concentration of H 2 S was maximum (938 ppm) at Day 3, and later on the concentration fluctuated based on feeding (Fig. 1A, Table S1). The overall environmental temperature, AD temperature, AD pressure and humidity were 34 Table S2). On Day 35 of the digestion, when maximum methanogenesis was observed, the concentration of organic carbon (OC) and total nitrogen (TN) in the fermentation pulp were 15.48% and 1.22%, respectively, whereas the concentration of OC and TN in the slurry (CD + seed sludge; Day-0) were 34.39% and 1.96%, respectively ( Table 2). The overall C/N ratio of the feedstock also gradually decreased with the advent of anaerobic digestion process, and found lowest (12.7:1) at Day 35. Similarly, the amount of non-metallic element (phosphorus and sulfur) and heavy metals (chromium, lead and nickel) content significantly decreased at the Day 35 of the digestion process (Table 2). However, the amount of zinc and copper did not vary significantly throughout the digestion period ( Table 2). As shown in Fig. 2, number of operational taxonomic units (OTUs) were significantly positively correlated with AD bacterial α-diversity (Observed, Chao1 and Shannon indices; r were 1.00, 1.00, and 0.57, respectively, p < 0.05). Similarly, we found significant positive correlation between CH 4 concentration (%) and pH (r = 0.92, p < 0.001), CH 4 concentration (%) and O 2 level (%) (r = 0.54, p < 0.05), and CH 4 concentration (%) and environmental temperature (°C) (r = 0.57, p < 0.05) (Fig. 2). We also found positive correlation between the AD CO 2 (%) and H 2 S (%), environmental temperature and digester temperature, environmental temperature and digester pressure, digester temperature and digester pressure (r = 0.64, p < 0.01; r = 0.71, p < 0.01; r = 0.72, p < 0.01; r = 0.89, p < 0.01, respectively). Interestingly, we also found a negative correlation between CO 2 and humidity, O 2 and humidity, environmental temperature and humidity, and digester temperature and humidity (Fig. 2).
Microbiome composition and diversity in the AD. The samples (n = 16) were categorized into four groups (Group-I, Group-II, Group-III and Group-IV) based on the level of energy production (CH 4 %) ( Table 1).

Figure 1.
Dynamic changes in the physicochemical parameters of the anaerobic digester (AD) over the study period. The R package, ggplot2 (https:// cran.r-proje ct. org/ web/ packa ges/ ggplo t2/ index. html) was used to visualize the line graphs.    Table 2. Physicochemical properties of raw cow dung, slurry and active sludges. CD: raw semi-solid cow dung which was mixed with water; Slurry: mixture of CD and seed sludge from previous biogas plant; AS: slurry from the AD when the gas production rate was the highest (i.e. 74.1% CH 4 concentration); and BDL: below detection limit.

Parameters Cow dung (CD; Day-0) Slurry (CD + seed sludge; Day-0) Active sludge (AS; Day-35, highest CH 4 concentration)
Moisture (%) 23 www.nature.com/scientificreports/ coverage depth for most samples was sufficient to capture the entire microbial diversity. We also observed significant differences in the microbial community structure among the four metagenome groups (i.e., beta diversity analysis). Principal coordinate analysis (PCoA) at the strain level (Fig. 3C), showed a distinct separation of samples by the experimental groups. Besides, we found significant (p = 0.032, Kruskal Wallis test) differences in the abundance of ARGs and metabolic functional genes/pathways (DataS2) which could strongly modulate the level of energy production through microbiome dysbiosis in the AD. In this study, on an average 0.43% WMS reads (assigned for r RNA genes) mapped to 28, 110 and 552 bacterial phyla, orders and genera respectively, and relative abundance of the microbiome differed significantly (p = 0.034, Kruskal-Wallis test) across the metagenome groups (Data S1). We observed significant shifts/dysbiosis in the microbiome composition at strain level. The PS analysis detected 2,513 bacterial strains across the four metagenomes, of which 768, 1421, 1819 and 1774 strains were found in Group-I, Group-II, Group-III and Group-IV metagenomes, respectively. Only, 18.34% bacterial strains were found to be shared across the four energy producing metagenomes (Fig. 3D, Data S1). The archaeal fraction of the AD microbiomes was represented by 5, 17, 61 and 343 archaeal phyla, orders, genera and strains, respectively, and the relative abundance of these microbial taxa also varied significantly among the four metagenome groups (Fig. 3E). Remarkably, 95.90% (329/343) of the detected archaeal strains shared across these metagenomes (Fig. 3E, Data S1). In addition, 472, 536, 535 and Observed, Chao1, Shannon, and Simpson are the α-diversity indices. *Significant level (■ p < 0.1; *p < 0.05; **p < 0.01; ***p < 0.001). The R package, chart. Correlation function of PerformanceAnalytics (https:// cran.rproje ct. org/ web/ packa ges/ Perfo rmanc eAnal ytics/ index. html) was used to analyze and visualize the plot. www.nature.com/scientificreports/ 536 strains of bacterial viruses (bacteriophages) were identified in Group-I, Group-II, Group-III and Group = IV metagenomes, respectively (Data S1).
Microbial community dynamically changed over time in the AD. Significant changes in the abundances of core microbial groups were observed under anaerobic condition of the AD. At phylum level, the AD metagenome was dominated by Firmicutes, Bacteroidetes, Proteobacteria, Actinobacteria, Spirochaetes and Fibrobacteres comprising > 93.0% of the total bacterial abundances. Among these phyla, Firmicutes was the most abundant phylum with a relative abundance of 41.94%, 37.99%, 40.40% and 38.96% in Group-1, Group-II, Group-III and Group-IV, respectively. The relative abundance of Bacteroidetes (from 37.87% in Group-I to 22.40% in Group-IV) and Actinobacteria (from 3.94% in Group-I to 3.30% in Group-IV) gradually decreased with the advance of AD digestion time. Conversely, relative abundance of Proteobacteria (from 8.08% in Group-I to 18.92% in Group-IV) and Spirochaetes (from 1.28% in Group-I to 3.70% in Group-IV) gradually increased with the increase of anaerobic digestion time in AD. The rest of phyla also differed significantly across these four groups keeping comparatively higher relative abundances during highest CH 4 producing stage (Group-IV) of the AD (Fig. 3F). The relative abundances of archaeal phylum Euryarchaeota, steadily increased with the increasing demand of energy (lowest relative abundance in Group-I and highest relative abundance in Group-IV) (Fig. 3F).
Similarly, Clostridiales and Bacteroidales were identified as the top abundant order in Group-1, Group-II, Group-III and Group-IV with a relative abundance of 32.37%, 27.81%, 29.22% and 27.87%, and 32.49%, 27.42%, 27.53% and 14.94%, respectively (Data S1). The structure and relative abundances of the bacteria at the genus level also showed significant differences (p = 0.031, Kruskal-Wallis test) across the study groups. In Group-I, Group-II and Group-III metagenomes, Bacteroides was the most abundant bacteria with a relative abundance of 18.10%, 14.90% and 15.16%, respectively, but remained lower (8.31%) in Group-IV samples. Clostridium was found as the second most predominant bacterial genus, and the relative abundance of this bacterium was 11.92%, 11.13%, 11.73% and 12.15% in Group-1, Group-II, Group-III and Group-IV, respectively. The relative abundance of Ruminococcus, Eubacterium, Parabacteroides, Fibrobacter, Paludibacter, Porphyromonas and Bifidobacterium gradually decreased with the increase of energy (CH 4 ) production rate, and remained lowest in Group-IV. Conversely, Candidatus, Bacillus, Treponema and Geobacter showed an increasing trend in their relative abundances gradually with the advance of digestion time and remained lowest in relative abundances in Group-IV. The rest of the bacterial genera had lower relative abundances in four metagenomes of the AD (Fig. S3, Data S1). In our present study, Methanosarcina was the most abundant archaeal genus, and the relative abundance of this genus remained two-fold higher in Group-III (35.84%) and Group-IV (36.53%) compared to Group-I (17.52%) and Group-II (18.32%). Notably, the relative abundance of Methanoculleus was found higher in Group-II (11.59%) and Group-IV (13.80%) and lowest in Group-I (3.46%). Likewise, Methanobrevibacter was predominantly abundant at the initial phage of digestion (highest in Group-I; 19.35%) and remained lowest in abundance in the top CH 4 producing metagenome (Group-IV; 5.01%). Besides these genera, Methanothermobacter (5.30%), Methanosaeta (5.16%), Methanococcus (4.74%), Thermococcus (2.96%), Methanocaldococcus (2.53%), Pyrococcus (2.35%), Methanosphaera (2.32%) Methanococcoides (2.10%) and Archaeoglobus (2.01%) were the predominantly abundant archaeal genera in Group-I samples and their relative abundances gradually decreased with the increase of energy production (Fig.  S4, Data S1). On the other hand, Methanoregula (6.43%), Methanosphaerula (2.99%), Methanoplanus (2.37%) and Methanohalophilus (1.39%) were the most abundant archaeal genera in Group-IV metagenome. The rest of the genera remained much lower (< 1.0%) in relative abundances but varied significantly across the four metagenomes ( Fig. S4, Data S1).

Identification of potential indicator species and their co-occurrence.
To identify microbial taxa (bacteria and archaea) that could discriminate across the four metagenome groups of the AD in terms of energy production (% CH 4 ), the indicator species analysis (ISA) was performed both in individual group and combination basis, as shown in (Fig. 5). Indicator species were those which were significantly more abundant and present in all samples belonging to one group, and also absent or low abundance in the other group ( Fig. 5, Data S1). The core taxa were selected based on their relative frequency (≥ 75% occurrence in each of the four groups) (Data S1). Although, 26, 3 and 19 indicator species were found in Group-I, Group-II and Group-IV, respectively, and no indicator species were identified in Group-III ( Fig. 5A, Data S1  5A; Data S1). Considering the combined group effects of the indicator species associated with energy production, our analysis revealed that Methanosarcina vacuolate, Dehalococcoides mccartyi, Methanosarcina sp. Kolksee and Methanosarcina barkeri in Group-III + Group-IV (top CH 4 producing groups) having IVs of 0.88, 0.887, 0.879 and 0.879, respectively were highly specific for energy production (Fig. 5B, Data S1). All of the indicator phylotypes displayed reduced abundance in the initial stage of biogas production (Group-I and Group-II, lower CH 4 production rate) compared to their increased relative abundance up to Day 35 of the experiment (in Group-III and Group-IV) (Data S1). We then visualized networks within each metagenome group of the AD for both positive and negative cooccurrence relationships (Fig. 6, Data S1). The correlation networks analysis was performed based on the significantly altered species (n = 106) in different groups as revealed by indispecies analysis. This network analysis explored significant association (p = 0.021, Kruskal-Wallis test) in the co-occurrence patterns of the energy producing microbial taxa (species and/or strains) based on their relative abundances in four metagenome groups. In the correlation network of four metagenomes; Group-I to Group-IV), Firmicutes and Bacteroidetes exhibited strongest relation. The resultant network consists of 106 nodes (17 in Group-I, 58 in Group-II, 5 in Group-III and 26 in Group-IV) which were clearly separated into four modules/clusters (Fig. 6). Taxa in the same group may co-occur under the same AD conditions (temperature, O 2 and H 2 S percentage, pressure and humidity). Across . Microbial abundance and the correlation between different physicochemical parameters and microbial relative abundance. (A) The species-level taxonomic abundance of microbiome. Stacked bar plots showing the relative abundance and distribution of the dominant abundant species (43 taxa), with ranks ordered from bottom to top by their increasing proportion among the four metagenomics groups. Each stacked bar plot represents the abundance of each strain in each sample of the corresponding category. The relative abundances of archaeal species (red color legends) steadily improved as energy demand increased (lowest relative abundance in Group-I and highest relative abundance in Group-IV). In contrast, the relative abundances of most of the known bacterial species gradually decreased with the passage of time (increased energy production), and mostly remained higher in Group-I and lower in Group-IV. (B) Spearman's correlation analysis between different physicochemical parameters [pH, CH 4 (%), CO 2 (%), O 2 (%), Others (%), H2S (ppm), Env_Temp (°C) (Environmental_Temperature), Dig_Temp (°C) (Digester_Temperature), Dig_pre (mb) (Digester_pressure), Humidity (%)] and dominant microbial relative abundance at species level. The numbers display the Spearman's correlation coefficient (r). Blue and red colors indicate positive and negative correlation, respectively. The color density, circle size, and numbers reflect the scale of correlation. *Significant level (*p < 0.05; **p < 0.01; ***p < 0.001). The R packages, Hmisc (https:// cran.r-proje ct. org/ web/ packa ges/ Hmisc/ index. html) and corrplot (https:// cran.r-proje ct. org/ web/ packa ges/ corrp lot/ vigne ttes/ corrp lot-intro. html) were used respectively to analyze and visualize the data.  6). However, when moving down to the species-level in microbiome co-occurrence in the AD, keystone taxa were much more consistent between networks with different correlation cutoffs. These results reveal that applying the same conditions in the AD for energy production, network elements must happen under careful consideration of the parameters used to delineate co-occurrence relationships. The positive correlations between Group-I and Group-II were observed among the microbiomes of the AD while Group-III and Group-IV showed negative correlations in terms of energy production with the microbial taxa of other two groups (Fig. 6). These findings therefore suggest that different strains of Euryarcheota and Firmicutes phyla were negatively correlated but associated with highest level of energy production (highest % of CH 4 ; Group-IV).
Genomic functional potentials of the anaerobic microbiomes. In this study, there was a broad variation in the diversity and composition of the antimicrobial resistance genes (ARGs) (Fig. 7, Data S2). The results of the present study revealed significant correlation (p = 0.0411, Kruskal-Wallis test) between the relative abundances of the detected ARGs and the relative abundance of the associated bacteria found in four metagenomes (Data S2). ResFinder identified 45 ARGs belonged to eight antibiotic classes distributed in 2,513 bacterial strains (Data S2). The Group-III microbiomes harbored the highest number of ARGs (42), followed by Group-II (38), Group-IV (29) and Group-I (22) microbes (Fig. 7, Data S2). The tetracyclines (doxycycline and tetracycline) resistant gene, tetQ had the highest relative abundance (23.81%) in Group-I associated bacteria followed by Group-II (22.85%), Group-III (16.49%) and Group-IV (6.73%)-microbes. Macrolides (erythromycin and streptogramin B) resistant genes such as mefA (16.80%), mefB (15.32%) and msrD (11.10%) had higher rela-  ((circle size) is proportional to the mean relative abundance in that group of AD and p-values (circle color) were plotted. Circles colored in grey indicate insignificant taxa. Higher indicator values (IVs) suggested better performances in the microbial signature of the assigned taxa. The R package, Indicspecies (https:// cran.r-proje ct. org/ web/ packa ges/ indic speci es/ index. html) was used to analyze the data, and ggplot2 (https:// cran.r-proje ct. org/ web/ packa ges/ ggplo t2/ index. html) was used to draw these bubble groups. www.nature.com/scientificreports/ tive abundances in highest CH 4 producing metagenome compared to other metagenome groups (Fig. 7A). The broad-spectrum beta-lactams resistant gene, cfxA2-6 was found as the common ARG among the microbiomes of four metagenomes, displaying the highest relative abundance (35.58%) in inoculum (Group-I) microbiota followed by Group-II (23.09%), Group-III (8.02%) and Group-IV (0.14%) microbiomes. The rest of the ARGs also varied in their expression levels across the four metagenomes, being more prevalent in the Group-III microbiomes (Fig. 7A, Data S2).
Most of the dominant archaeal species had stronger positive correlation with sul4, tet(W), tet (44), tet(C), Inu(B), tet (32), catQ, tetM) and Inu(D) (Spearman correlation; r > 0.5, p < 0.01), and were negatively correlated with tet(Q), tet (40), Inu(C), cfxA, cfxA2, cfxA3, cfxA4, cfxA5 and cfxA6 (Spearman correlation; r > 0.6, p < 0.01) (Fig. 7B, Data S2, Fig. S6). In contrast, the dominating bacterial species of the AD had significant positive correlation with cfxA, cfxA2, cfxA3, cfxA4, cfxA5, cfxA6, tet(Q), blaACI-1, tet (40) and Inu(C) (Spearman correlation; r > 0.6, p < 0.01) (Fig. 7B, Data S2, Fig. S6). Analyzing the correlation between physicochemical parameters and ARGs classes, we found that the digester pH, CH 4 concentration, pressure and temperature were positively correlated with macrolide-lincosamide-streptogramin (MLS), tetracyclines and phenicol while most of the physicochemical parameters were negatively correlated with sulfonamides class of antibiotics (Spearman correlation; r > 0.6, p < 0.01) (Fig. 7C, Data S2, Fig. S7A). Likewise, the digester pH, CH 4 concentration, pressure and temperature showed strongest positive correlation with MLS resistance MFS and ABC efflux pumps associated mechanisms (Spearman correlation; r > 0.5, p < 0.01) followed by chloramphenicol acetyltransferases and tetracycline resistance ribosomal protection proteins related functions (Fig. 7D, Data S2). Conversely, Class A betalactamases and sulfonamide resistance dihydropteroate synthases were the negatively correlated mechanisms with digester physicochemical parameters like pH and CH 4 (Fig. 7D, Data S2). The mechanisms of the antimicrobial resistance also varied between the dominant bacterial and archaeal species (Fig. S7B). For Figure 6. Microbiome co-occurrence in the AD within the four metagenomic groups. Correlation networks of significantly altered species based on the indicator species analysis. Spearman's absolute correlation coefficients ≥ 0.5 and p-values < 0.05 were retained. The node size is proportional to the mean abundance of the species. Nodes are colored by taxonomy with labelled genera names. The positive correlation is represented by the green line, while the negative correlation is represented by the red line. The nodes were grouped according to their higher abundance in each group (It means that a species can be found in all or at least one group. But they were kept in such a group where it was highly abundant). The microbiomes of the AD showed positive associations between Groups I and II, while Groups III and IV showed negative correlations in terms of energy production with the microbial taxa of the other two groups. The R packages, Hmisc (https:// cran.r-proje ct. org/ web/ packa ges/ Hmisc/ index. html) package was used to analyze correlation and Gephi (https:// gephi. org/) software was utilized to visualize the network plot. www.nature.com/scientificreports/ instance, most of the bacterial species showed highest positive correlation Class A betalactamases (Spearman correlation; r > 0.5, p < 0.01) while the enriched archaeal species showed negative correlation with this functional mechanism of antibiotics. The methanogenic archaeal species predominantly identified however had significant positive correlation with tetracycline resistance ribosomal protection proteins related functions, 23S rRNA methyltransferases and MLS resistance MFS efflux pumps (Fig. S7B, Data S2). In addition to these ARGs, the highest CH 4 producing microbiomes were enriched with the higher relative abundance of genes coding for cobalt-zinc-cadmium resistance (18.85%), resistance to chromium compounds (12.17%), arsenic (6.29%), zinc (4.96%) and cadmium (3.26%) resistance compared to the microbes of other three metagenomes (Fig. S8, Data S2). By comparing the possible mechanisms of the detected ARGs, we found that antibiotic efflux pumps associated resistance had the highest level of expression in the anaerobic microbiomes of the AD followed by antibiotic inactivation, enzymatic inactivation and modification, antibiotic target protection/alteration, and folate pathway antagonist-attributed resistance mechanisms (Fig. S8, Data S2). Notably, the physicochemical properties of the AD including pH, CH 4 concentration (%), pressure, temperature and environmental temperature were found to be positively correlated with metal such as arsenic, copper, cadmium, mercury, chromium compounds and zinc resistance pathways (Spearman correlation; r > 0.5, p < 0.01). (Fig. 8A, Data S2). Remarkably, most of the www.nature.com/scientificreports/ methane producing dominant archaeal species showed significant positive correlation with these metal resistance pathways (Spearman correlation; r > 0.7, p < 0.01). (Fig. 8B, Data S2) while the predominant bacterial species to be associated with methane production showed stronger negative correlation with most of the metal resistance pathways (Spearman correlation; r > 0.7, p < 0.01). (Fig. 8B, Data S2). Functional metabolic profiling of the gene families of the same KEGG pathway for AD microbiomes revealed significant differences (p = 0.012, Kruskal-Wallis test) in their relative abundances, and positive correlation with level of energy production (Fig. 9, Data S2). Among the detected KO modules, genes coding for CHO metabolism and genetic information and processing were top abundant, however did not vary significantly across the metagenome groups. Remarkably, the relative abundance of genes coding for energy metabolism, xenobiotics biodegradation and metabolism, butanoate metabolism, citrate synthase (gltA), succinyl-CoA synthetase subunits (sucC/D), pyruvate carboxylase subunits (pycA) and nitrogen metabolism gradually increased with the increasing rate of CH 4 production, and had several-fold over expression among the microbiomes of Group-IV. Conversely, fumarate hydratase (fumA/B), malate dehydrogenase (mdh) and bacterial secretion system associated genes were predominantly overexpressed in Group-I related microbiomes which gradually decreased with advance of digestion process, and remained more than two-fold lower expressed in the peak level of CH 4 production (lowest in Group-IV) (Fig. 9A, Data S2). We also detected 41 statistically different (p = 0.033, Kruskal-Wallis test) SEED functions in the AD microbiomes. Overall, the top CH 4 producing microbiomes (Group-III and Group-IV) had higher relative abundances of these SEED functions compared lower CH 4 producing microbiomes (Group-I and Group-II), except for regulation of virulence (highest in Group-I microbes; 17.08%), gluconeogenesis (highest in Group-I microbes; 16.27%) and transposable elements (highest in Group-I microbes; 17.28%) (Fig. 9B, Data S2). The Group-IV-microbiomes (highest CH 4 producing) were enriched in genes coding for tetrapyrroles (17.42%), one carbon (10.29%) and biotin (4.55%) metabolism, oxidative (18.76%) and osmotic (9.94%) stress, proteolytic pathway (7.74%), MT1-MMP pericellular network (6.45%), acetyl-CoA production (5.33%) and motility and chemotaxis (3.13%) compared to the microbes of the other metagenomes. The Group-I microbiomes however had a higher abundance of SEED functions involved in protection from ROS (16.28%), heat shock (18.31%) and NAD and NADP (19.03%) (Fig. 9B, Data S2). Analyzing correlation between AD physicochemical parameters and microbial genomic functions, we found that SEED subsystems such as cofactors, vitamins, prosthetic groups, www.nature.com/scientificreports/ pigments, membrane transport, motility and chemotaxis, respiration, stress response, dehydrogenase_complexes, dehydrogenase_ kinases, glycolysis and gluconeogenesis including archaeal enzymes, glyoxylate bypass, pentose phosphate pathway, pyruvate:ferredoxin oxidoreductase, pyruvate metabolism I and II had significant positive correlation with digester pH and CH 4 concentration (%) (Spearman correlation; r > 0.6, p < 0.01). (Fig. 10A, Data S2). Simultaneously, the Kos like ascorbate and aldarate metabolism, citrate synthase (gltA), energy, glyoxylate, dicarboxylate, methane and nitrogen metabolism, xenobiotics biodegradation and metabolism, succinyl-CoA synthetase subunits (sucC/D) and oxidative phosphorylation were significantly positively correlated with digester pH, CH 4 concentration (%), pressure and temperature (Spearman correlation; r > 0.6, p < 0.01). (Fig. 10B, Data  S2).

Discussion
This study is the first ever approach to reveal the dynamic shifts in microbiome composition and abundances in different levels of biogas production under the anaerobic digestion system using the state-of-the-art WMS technology along with analysis of the physicochemical parameters in Bangladesh. Anaerobic digestion of organic wastes is favored by the metabolic activities of different types of microorganisms including bacteria and archaea 30 .
The physicochemical parameters assessment of the AD before and after digestion revealed that biogas production was in increasing trend up to Day 35 with periodic loading of slurry (Day 1, 2, 3, 8, 10, 14, 16, 24 and 35) at 1:1 ratio of raw CD and active sludge (Table S1). After 30 days of incubation, we found that cow wastes produced highest amount of biogas (on Day 35, CH 4 ; 74.1%), and corresponds to an average pH of 7.01. Earlier studies showed that optimum pH range in an AD is 6.8 to 7.2 supporting our present findings 31 . The CH 4 production level thereafter decrease gradually reaching 59.2% on Day 44 of processing (Fig. 1). Moreover, with the gradual increase of AD temperature and pH value (~ 7; neutral), the CO 2 and CH 4 concentration positively increased up to Day 35 which further correlated with the higher abundance of methanogenic bacteria. Biogas production chiefly depends on the content and chemical nature of biodegradable matter. Correlation analysis between the relative abundances of specific microbial taxa and digesters' physicochemical parameters provides a qualitative analysis method to explore the possible roles of and the interactions among these microbes. Although, most of the dominant bacterial species had negative correlation with digester pH, CH 4 concentration, humidity, temperature and pressure, the methanogenic archaea showed strongest positive correlation with these parameters. The biochemical parameter of cow waste (slurry) reflects the presence of high content of readily biodegradable organic matter in the first phases (up to Day 35) of anaerobic digestion 30 . The CO 2 concentration (%) found to be varied throughout the digestion process keeping an average value of 39.52%. Our analysis revealed that OC and TN content was higher at the time of loading of slurry in the AD compared to that of highest biogas production stage (at Day 35). The amount of carbon available of the substrate determines the maximum amount of CH 4 and CO 2 that can be formed by anaerobic digestion 9 . Conversely, C/N ratio remained lowest in this peak stage of CH 4 production. The OC is essential for bacterial growth, and determining of the C/N ratio is essential for optimal biogas production 32 . Moreover, the total content of phosphorus, Sulphur and heavy metals (chromium, lead and nickel) also remained lowest at this highest stage of biogas production (Day 35). The lowest chromium, lead and nickel concentration during highest CH 4 producing stage might be associated with their higher abundances of heavy metal resistance genes, and small concentrations of these metals found in the process are essential for microbial maintenance 33,34 . Certain specific metals such as Sulphur, cobalt and nickel serve as cofactors in the enzymes involved in the formation of CH 4 during anaerobic processing 35 . However, the minerals (e.g. zinc and copper) content of the AD did not vary throughout the digestion process revealing their important roles in various metabolic pathways of anaerobic digestion 33 . Since, CH 4 gradually declined after Day 35 of digestion irrespective differences in other physicochemical parameters, we collected samples for microbiome analysis only up to peak level CH 4 production (Day 35). The limitation of the study is that a complete correlation between heavy metal concentration and microbial community in terms of CH 4 is missing, so we actually do not know the effect of the heavy metals on the single microorganism, but to perform this kind of analysis, the methods should be also changed significantly. Energy production in the AD is carried out by a series of biochemical reactions including hydrolysis, acidogenesis, and acetogenesis performed by fast-growing bacteria through which complex organic macromolecules into volatile fatty acids (VFA), CO 2 , and H 2 . However, in many cases, these reactions occur at low pH, while methanogenesis is performed by slower-growing bacteria that thrive best at neutral pH 36 .
Remarkably, the pH, CH 4 concentration, pressure and temperature were the main driving factors for variations in microbial community structure in the AD as also reported in previous studies 25,30,36 . The within (alpha) and between (beta) sample diversity of the AD microbiomes showed that that microbial dysbiosis in the AD is closely linked to different levels of biogas production. Compared to loading phase of AD (Group-I), increased microbial diversity and species richness was observed in the later phases (Group-II, Group-III and Group-IV) of anaerobic digestion. Beta diversity also revealed a substantial microbial disparity in different levels of biogas production, and segregated the samples accordingly. Despite having higher taxonomic resolutions, the microbiomes of the AD remained inconsistent and fluctuates more in Group-II, Group-III and Group-IV than those of Group-I metagenome. The taxonomic annotations of the four groups of AD showed that they were a reservoir of bacteria, followed by archaea and viruses, which corroborated the findings of other studies 37,38 . Among the identified domains, bacteria dominated in abundance, comprising 81.80% of the total microbial populations, followed by archaea (15.43%), while viruses (2.77%) comprised the least abundant population. The observed high bacterial abundances suggest their crucial metabolic roles in biomass conversion and other reactions within the reactor systems 38 . The identified affiliates of archaea were mostly consumers of smaller substrates that were generated by the bacterial taxa. The archaeal species are able to use different methanogenic routes to convert the substrates into methane gas. Nevertheless, the main roles of the identified less abundant bacteriophages were unclear, though the strains could have been active in degrading other microbial cells in the AD systems 37 www.nature.com/scientificreports/ In this study, the most of the dominant methanogenic archaeal species (n = 12) were belonged to phylum Euryarchaeota and class Methanomicrobia and these findings are consistent with several previous studies 17, 38 on biogas plants where phylum Euryarchaeota (class Methanomicrobia) was reported as the dominant archaea contributing to the energy production.
The methanogenic archaea are obligate anaerobes that can reduce simple substrates to methane (CH 4 ) to produce cellular energy. Most of the dominant archaeal species belonged to Methanosarcinales and Methanobacteriales orders which are syntrophic acetate oxidization digesters 39 . M. vacuolata can play roles in energy production by utilizing methylamines, and resistance to harsh conditions 40 . Methanoculleus sp. are acetate-oxidising bacteria, which convert acetate directly to methane and of considerable importance in ammonia-rich engineered biogas processes 41 . Several members of Firmicutes and Bacteroidetes phyla are well known as fermentative and syntrophic bacteria, and can ferment cellulose and cellobiose to produce acetic acid, ethanol, CO 2 /H 2 and CH 4 42 . Moreover, the enriched consortia of Fibrobacteres, Spirochaetes, Actinobacteria, Proteobacteria, and Thermotogae phyla had negative correlations with physicochemical parameters, indicating some members working on acetate and butyrate syntrophic oxidation might be dominant in the AD 42 . The Candidatus Cloacimonetes phylum is reported to play a role in hydrolysis of cellulose and/or fermentation of hydrolysis products in AD 43 . In this study, Firmicutes, Bacteroidetes, Proteobacteria, Actinobacteria, Spirochaetes and Fibrobacteres were the most abundant bacterial phyla, and their relative abundance also varied according to the level of energy production in the corresponding sample groups. The bacterial phyla Firmicutes and Bacteroidetes appeared to dominate biogas communities in varying abundances depending on the apparent process conditions corroborating with earlier studies 43 . The observed community composition at phylum level, with dominance of Firmicutes, Bacteroidetes and Proteobacteria, is also in line with previous findings for biogas reactors 34 . The first three phases (hydrolysis/ cellulolysis, acidogenesis/fermentation and acetogenesis) of the anaerobic digestion are solely performed by fermentative Bacteria. In this study, the bacterial phyla Firmicutes (37.0-42.0%) and Bacteroidetes (22.0-38.0%) appeared to dominate biogas communities in varying abundances depending on the apparent process conditions. The initial phase of the AD digestion process involves hydrolytic reactions that convert large macromolecules into smaller substrates 37 . In addition, the nature and composition of the substrates, availability of nutrients and ammonium/ammonia contents can affect both the composition and diversity of the methanogenic archaea 37 . Only certain methanogenic archaea are able to synthesize CH 4 from the end products of bacterial fermentation. The performance and efficiency of these processes depend to a large extent on the presence of appropriate and adequate microorganisms along with the physicochemical conditions of the digester and quality of the substrates (organic materials) etc. The degree and rate of degradation (hydrolysis, fermentation, acetogenesis, and methanogenesis) and the biogas yield depend not only on the chemical and physical characteristics of the substrates, but also on the chosen process parameter such as temperature, humidity and retention time, that shape the composition of different microbial groups and communities that active in the process 34 . During the anaerobic process Clostridiales and Bacteroidales were identified as the top abundant order in all of the four metagenome groups. The large number and proportion of members of the Clostridiales order are indicative of the important role in the proper functioning of the microbial community in an AD fed with complex substrates. Numerous members of the Clostridiales order, and members of the Clostridiaceae family are capable of performing diverse fermentation pathways and known to use lactate and acetate as electron sources for the reduction of iron and cobalt under anaerobic condition. They primarily ferment sugars to organic acids through the reductive acetyl-CoA pathway, which is typical in acetogenic bacteria and in some Archaea 4 . In this process, CO 2 is reduced to CO and then converted to acetyl-CoA, H 2 serving as electron donor. It should additionally be noted that a large number of Clostridia actively produce H 2 , an important substrate for the hydrogenotrophic methanogens 4,44 . In anaerobic environments, Clostridiales has been reported as the main cellulose degrader 45 , and play important role in the hydrolysis step 34 . These findings are in line with many of the previous reports 37,46 who reported the potential roles of these bacterial taxa in the efficient and increased production of biogas under anaerobic condition 47 . In addition, bacteria belonging to the order Bacteroidales have been suggested to be involved in the degradation of lignocellulose materials, such as straw and hay, the chief component of cattle feed 47,48 . The identified affiliates of Bacteroidales belong to Bacteroidetes, and are majorly known to ferment carbohydrates and proteins, concomitantly releasing H 2 33 . Moreover, Bacteroides, the predominating genus in the AD coming from Bacteroidales order were observed to co-exist with methanogenic archaea possibly to increase energy extraction from indigestible plant materials 33 , and we hypothesize that they are the key drivers of the observed β-diversities since the relative abundance of this genus gradually decreased with the digestion process.
The digestion process of the AD was carried out by the integrated cross-kingdom interactions since both bacteria and archaea were simultaneously detected in this WMS-based study corroborating with several earlier reports 22,23,26 . The strain-level taxonomic profiling revealed that methanogenic archaeal strains were more prevalent than bacterial strains. The archaeal domain of the AD microbiomes was composed of different strains of methanogenic, hydrogenotrophic and thermophilic genera of Methanoculleus, Methanosarcina, Methanothrix, Methanobacterium and Methanobrevibacter genera. The current findings are corroborated with many of the earlier studies who reported that these genera to be predominantly abundant in the AD of manures 30 , and associated with biogas production under anaerobic conditions 49 . These methanogenic genera might reside in the microenvironments appropriate for anaerobic metabolism 26 and their presence has been reported in microbial communities producing biogas 50 . Members of Methanoculleus are hydrogenotrophic methanogens 51 , while Methanosarcina species or strains are mostly acetoclasic but also able to use H 2 24,52 . In addition, Methanosarcina spp. has been reported to have higher growth rates and tolerance to pH changes and could potentially lead to stable methanogenesis in the AD 10,53 . Methanobrevibacter was predominant in initial phase of digestion (Day 2 and 15) in the bioreactor, and are known to be hydrogenotrophic, by using CO 2 and H 2 as substrates to generate biomethane 24,50 . These archaeal genera are suggested to play vital role in hydrogenotrophic methanogenesis, and maintaining methanogenic community diversity 21 www.nature.com/scientificreports/ We found significant association across different methanogenic species with increased level of CH 4 production 54 , with particular emphasis in the highest CH 4 producing metagenome (Group-IV) of the AD. Mounting evidence suggested that, under anaerobic condition, the functional composition of the energy producing microbial communities are closely related to the physicochemical and environmental factors such as pH, temperature and pressure of the digester 10,30,38,55 . The ecological functions of microorganisms living in similar environments remained more similar, but the composition of microbial species (especially, the dominant species) performing the functions differed greatly corroborating with the results of several previous studies 10,24,56 . The methanogenic archaea are obligate anaerobes that can reduce simple substrates to methane (CH 4 ) to produce cellular energy. Most of the dominant archaeal species belonged to Methanosarcinales and Methanobacteriales orders which are syntrophic acetate oxidization digesters 39 . M. vacuolata can play roles in energy production by utilizing methylamines, and resistance to harsh conditions 40 . Methanoculleus sp. are acetate-oxidising bacteria, which convert acetate directly to methane and of considerable importance in ammonia-rich engineered biogas processes 41 . Several members of Firmicutes and Bacteroidetes phyla are well known as fermentative and syntrophic bacteria, and can ferment cellulose and cellobiose to produce acetic acid, ethanol, CO 2 /H 2 and CH 4 42 . Moreover, the enriched consortia of Fibrobacteres, Spirochaetes, Actinobacteria, Proteobacteria, and Thermotogae phyla had negative correlations with physicochemical parameters, indicating some members working on acetate and butyrate syntrophic oxidation might be dominant in the AD 42 . The Candidatus Cloacimonetes phylum is reported to play a role in hydrolysis of cellulose and/or fermentation of hydrolysis products in AD 43 .
The differences in ARGs profiles throughout different stages of anaerobic digestion were greatly affected by the taxonomic composition of microbiomes (archaea and bacteria), and the physicochemical parameters of the AD. The observed differences in microbiome composition and AD physicochemical properties are considered to be co-selection factors for antibiotic resistance 57 . Moreover, we also found that cross-resistance (where a single microbe is related to resistance to both antibiotics and metals) and coregulation (a shared regulatory system for antibiotic and metal resistance) may also lead to the increases in ARGs 55 . Therefore, to understand the mechanism that allows heavy metals to influence ARGs, the microbial community, heavy metal resistance genes, and their relationships with ARGs should be research forward. Overall, the diversity and abundance of the ARGs determined in this study were high, which clearly indicates that the unmonitored use of antibiotics and metals in dairy farms can subsequently increase the abundance of ARGs in AD systems. Our results also showed that common indicator bacteria such as E. coli, Salmonella and Staphylococcus species were not found in indicator species analysis, and these findings are supported by several previous reports of the absence of common indicator bacteria after 30 days digestion in the experimental AD at different temperatures (25-45 °C) 58 . Spearmen correlation analysis showed negative correlation between Group-II and Group-IV microbiomes. These findings are in line with the metabolic functional potentials of the AD microbiomes since fumarate hydratase (fumA/B), malate dehydrogenase (mdh) and bacterial secretion system associated genes were predominantly overexpressed in Group-I and Group-III microbiomes, which gradually decreased with advance of digestion process, and remained more than two-fold lower expressed in highest CH 4 production stage (Group-IV). The microbial communities present in the early phage of digestion process (Group-I) increased both composition and abundances in the second phase of digestion (Group-II metagenome) in ambient growth conditions of the AD. With the advance of digestion time, the increasing efficiency of anaerobic digestion creates a favorable environment for the methanogens, and thus found in higher composition and abundances in Group-III and Group-IV metagenomes 21,25 . In addition, the Group-IV-microbiomes showed higher genomic functional activities related to tetrapyrroles, one carbon and biotin metabolism, oxidative and osmotic stress, proteolytic pathways, MT1-MMP pericellular network, acetyl-CoA production, and motility and chemotaxis compared to the microbes of the other groups. The energy and one-carbon metabolism are associated with certain common central pathways which are universally present across in many bacterial species, despite the wide variety of growth substrates and conditions they utilize 59 . Evidence of environmental stress conditions prevailing within the bioreactors suggested that microorganisms underwent to oxidative stress corroborating our present findings 60 . Gut microbiota has considerable proteolytic power, converting ingested dietary protein. Identified bacterial species particularly Clostridia, Streptococci, Staphylococci, and Bacillus showed higher proteolytic activity which differed both in quantity and quality of protein degradation in different stage anaerobic digestion, and these findings are in line with many earlier studies 61 . Though, these findings also support the taxonomic dysbiosis of microbiomes in Group-IV metagenome, however further comprehensive study is needed to elucidate the modulation of microbiome shifts, their functional potentials and genomic expression using a larger dataset. It is therefore imperative to carryout future research on biogas-producing microbiomes which will further help in the development and improvement of anaerobic digestion models and enhance AD efficiency and stability. Such research will benefit greatly from the continued improvement of the omics technologies, particularly metaproteomic and metametabolomics, and standardized methods of bioinformatics analyses.

Conclusions
The level of biogas production increased gradually up to Day 35 (highest CH 4 concentration), and declined thereafter under controlled environment of the AD. The physicochemical analysis showed significant positive correlation with the level of energy production and digester pH, O 2 level, and environmental temperature. In the AD microbiomes, Firmicutes, Bacteroidetes, Fibrobacteres, Spirochaetes, Actinobacteria, Candidatus Cloacimonetes, Proteobacteria, and Thermotogae were the dominant bacterial phyla while Euryarchaeota was the only predominant archaeal phylum. The indicator species analysis revealed that microbiomes of Group-III and Group-IV, for instance, M. vacuolate, D. mccartyi, Methanosarcina sp. Kolksee and M. barkeri were highly specific species for energy production. The correlation network analysis of the indicator species showed that different strains of Euryarcheota and Firmicutes phyla were negatively correlated but associated with highest level of energy www.nature.com/scientificreports/ production. The AD microbiomes harbored different ARGs of which genes coding for tetracyclines (doxycycline and tetracycline), macrolides (erythromycin and streptogramin B) and beta-lactams resistance remained over expressed. Additionally, the physicochemical properties of the digester including pH, CH 4 concentration (%), pressure, temperature and environmental temperature were found to be positively correlated with these genomic functional potentials and distribution of ARGs and metal resistance pathways. Taken together, a high-level abundance of the dominant and indicator microbial communities, associated genomic functional potentials, and their correlation with the physicochemical parameters of the AD were found to be driving factors for high proportion of the energy production.

Methods
Digester setup and experiment design. The experiment was conducted using an AD plant prepared with provisions to measure temperature, slurry and gas sample collection and substrate charging. The biogas plant consisted of a digester, inlet-chamber, three slurry outlet pipes, gas outlet pipe and thermometer (Fig. S1). The volume of the AD was 3000 L and made of flexible polyvinyl chloride (PVC) fabric with thickness: 1.2 Mm. There was a specialized ball valve which ensures the anaerobic condition within the digester and control the flow of substrates. The temperature of the AD was monitored using a probe. The experiment was conducted for 45 days (Day 0 to Day 44). The AD was loaded 12 times throughout the experiment period with a total of 1,087 kg raw CD (highest input volume = 375 kg and lowest input volume = 35 kg) (Table S1). Initially, the digester was started with the highest loading dose of 375 kg feedstock where the ratio of raw CD and active sludge was 1:1. The raw semi-solid CD was mixed with seed sludge from previous biogas plant (slurry) before charged into the AD. The AD was portable, light in weight, low-priced and retained more heat inside.  (Table 1). Therefore, highest CH 4 production (74.1%) was recorded at Day 35 of the digestion process. Thereafter, the CH 4 concentration gradually declined (up to Day 44), and thus, no samples were collected during this stage (Day 36-Day 44) for metagenomic analysis. The TN content was measured by micro-Kjeldahl method 62 while phosphorus, potassium, heavy metals (Lead, Zinc, Nickel, Cadmium, Chromium), OC, Sulphur and moisture content were determined by spectrophotometric molybdovanadate 63 , flame photometric 64 , atomic absorption spectrophotometric 63 , wet oxidation, turbidimetric and gravimetric methods from the Department of Soil Science, University of Dhaka. The detection limit of metals was of the order of 0.1 μg L −1 30 . www.nature.com/scientificreports/ variance (PERMANOVA) with 999 permutations to estimate a p-value for differences among groups 72 . Phyloseq and Vegan packages were employed for these statistical analyses 73 . Dominant microbial community were also determined where ≥ 1% read in at least one samples were considered. After filtering there were 43 taxa remained in 16 samples. Spearman's correlation analyses between dominant microbiotas and physicochemical parameters were performed in rcorr function of Hmisc 74 and corrplot function of corrplot 75 R package. Indicator species specific to a given sample group (having ≥ 1000 reads assigned to a taxon) were identified based on the normalized abundances of species using the R package indicspecies 76 , and the significant indicator value (IV) index was calculated by the 999-permutation test. Larger IV indicates greater specificity of taxa and p < 0.05 was considered statistically significant 76 . Network analysis was used to explore the co-occurrence patterns of the energy producing bacterial and archaeal taxa across the metagenome groups. In addition, the Spearman's correlation coefficient and significance tests were performed using the R package Hmisc. A correlation network was constructed and visualized with Gephi (ver. 0.9.2). We detected the antibiotics resistance genes (ARGs) among the microbiomes of four metagenomes through the ResFinder 4.0 database (https:// bitbu cket. org/ genom icepi demio logy/ resfi nder/ src/ master/). The ResFinder database was integrated within AMR++ pipeline (https:// megar es. meglab. org/ amrpl usplus/ latest/ html/ v2/) 77 to identify the respective antimicrobial resistance genes and/ or protein families 26 . Functional profiling of the microbiomes. In addition to taxonomic annotations, the WMS reads were also mapped onto the Kyoto Encyclopedia of Genes and Genomes (KEGG) database 78 , and SEED subsystem identifiers, respectively, on the MR server 66 for metabolic functional profiling. The functional mapping was performed with the partially modified set parameters (e-value cutoff: 1 × 10 −30 , min. % identity cutoff: 60%, and min. alignment length cutoff: 20) of the MR server 26 .
Statistical analysis. To evaluate differences in the relative percent abundance of taxa in AD (Group-I, Group-II, Group-III and Group-IV) for PS data, we used the non-parametric test Kruskal-Wallis rank sum test. Pairwise Pearson's correlation was performed to reveal the association between each of the two physicochemical parameters of the AD and microbial alpha diversity characteristics. We normalized the gene counts by dividing the number of gene hits to individual taxa/function by total number of gene hits in each metagenome dataset to remove bias due to differences in sequencing efforts. The non-parametric test Kruskal-Wallis rank sum tests were also performed to identify the differentially abundant SEED or KEGG functions (at different levels), and antimicrobial resistance (ARGs) in four metagenomes. All the statistical tests were carried out using IBM SPSS (SPSS, Version 23.0, IBM Corp., NY USA). To calculate the significance of variability patterns of the microbiomes (generated between sample categories), we performed PERMANOVA (Wilcoxon rank sum test using vegan 2.5.1 package of R 3.4.2) on all four sample types at the same time and compared them pairwise. A significance level of alpha = 0.05 was used for all tests. Spearman's correlation of dominant microbial community, KEGG pathways and SEED subsystem functional data to physical parameters of the AD were performed in Hmisc 74 and corrplot 75 R packages as described above.

Data availability
The sequence data reported in this article have also been deposited in the National Center for Biotechnology Information (NCBI) under BioProject accession number PRJNA668799. Supplementary materials supporting the results of the study are available in this article as Data S1 and S2, Figs. S1-S8, and Table S1.