Comparison of Bacterial and Archaeal Microbiome in Two Bioreactors Fed with Cattle Sewage and Corn Biomass

The bacterial and archaeal communities of two full-scale biogas producing plants (P1 and P2), associated with a 999 kW cogeneration unit, both located in North Italy, were analyzed at start up and fully operating phases, by means of various molecular approaches: (i) Automated Ribosomal Intergenic Spacer Analysis; (ii) cloning and sequencing of PCR amplicons of archaeal genes 16Srrna and mcrA; (iii) 16S rDNA high throughput next generation sequencing. P1 and P2 use the same technology and both were fed with cattle manure and corn silage. During the study of P1 also the post digester (fed with pig manure) was analyzed. The aim of this research was to characterize the bacterial and archaeal communities in two very similar plants to profile the core microbiota. The results of this analysis highlighted that the two plants (producing comparable quantities of volatile fatty acids, biogas, and energy) differed in anerobic microbiota (Bacteria and Archaea). Notably the methanogenic community of P1 was dominated by the strict acetoclastic Methanosaeta (Methanothrix) (up to 23.05%) and the unculturable Candidatus Methanofastidiosum (up to 32.70%), while P2 was dominated by the acetoclastic, but more substrate-versatile, Methanosarcina archaeal genus (49.19%). The data demonstrated that the performances of plants with identical design, in similar operating conditions, yielding comparable amount of biogas (average of 7237 m3 day−1 and 7916 m3 day−1 respectively for P1 and P2), VFA (1643 mg L−1 and 1634 mg L−1) and energy recovery (23.90–24 MWh day−1), depend on the stabilization of effective and functionally optimized methanogenic communities, but these communities are taxonomically different in the two biodigesters.


Introduction
In the last years, in order to reduce the consumption of fossil fuels and the environmental impact of methane emissions from zoogenic fluids, the biogas production from organic wastes has been taken into much more consideration as a solution and, consequently, it greatly increased [1]. The biogas produced by anaerobic digestion from such renewable sources is transformed in fuels that can be used for energy production and transport.
Despite the high operating costs of these plants and the necessity of significant financial incentives, the anaerobic digestion (AD) is a well-established technology, it remains a fundamental energy source in the emerging market for renewable energy and it is considered, in the circular economy perspective, a key technology on the path of independence from fossil fuels [2,3]. Italy is one of the leading countries in the biogas production from agricultural wastes, sewage and manure substrates [4].
As many commercial AD facilities suffer from several problems like long start up times or system instability, the optimization of the process is an important deal. The efforts to improve the anaerobic digestions process had formerly been focused on the evaluation of substrates and the implementation of plants characteristics, but it was soon hypothesized that the efficiency of the plants was more likely dictated by the successful stabilization of effective core microbiotas [5,6] composed by selected microorganisms derived from the incoming substrates. Some researchers have argued that it is difficult to define a single core microbiota among full-scale reactors [7]. However, what is clear is that the governing power of a core microbiota to AD performance is still poorly understood.
The overall process of biogas production is due to the complex, synergistic, interaction among four functional guilds performing the different reactions (hydrolysis, acidogenesis, acetogenesis, methanogenesis). Early in the process, microbial enzymes hydrolyse polysaccharides, proteins and fats into sugars, amino acids and fatty acids. These compounds are then converted into volatile fatty acids (VFAs) and alcohols that are further converted into acetic acid, CO 2 and hydrogen. Methanogenic microorganisms belonging to the Archaea domain carry out the final step, the production of methane [8,9]. Acetoclastic and hydrogenotrophic methanogens, transforming acetic acid, CO 2 and hydrogen into methane are particularly sensitive to process conditions, like temperature, pH, redox and inhibitors as ammonium and the sulfidrilic acid produced through the reduction of sulphates [10][11][12]. The methanogenic Archaea, mostly represented in biodigesters, belong to Methanomicrobiales, Methanosarcinales and Methanococcales orders [13][14][15][16][17][18].
The process of biogas formation has been thoroughly dissected, however a little number of prokaryotes involved in anaerobic digestion have been so far isolated, as most microbial species are not culturable, and this had slowed, in the past, the acquisition of information. Methods based on DNA analysis demonstrated to be the best way to study AD. The gene encoding ribosomal 16S RNA in Bacteria and Archaea is the most used target sequence for many molecular techniques, together with intergenic spacer (ITS) between the small (16S) and large (23S) rRNA genes in the rrna operon. In addition, methanogens can be identified by targeting the gene mcrA, encoding ɑ-subunit of methyl-coenzyme M reductase (MCR), the key enzyme involved in final steps of methanogenesis. Unlike other enzymes involved in methanogenic metabolism, MCR is highly specific for methanogens [19,20].
The molecular techniques used for microbiome analysis include ARISA (Automated Ribosomal Intergenic Spacer Analysis) analyzing ITS, clone library sequencing and NGS (Next Generation Sequencing). ARISA, based on the different length of ribosomal intergenic spacers, is a molecular technique for characterizing non-cultivable microbial communities directly from environmental DNA (the whole DNA extracted from samples). It has been widely applied for a rapid analysis of biodiversity in large samples, but suffers the limitation of not allowing the recognition, at the taxonomic level, of the species present in the sample [21][22][23].
The sequencing of 16Srrna or mcrA genes clone libraries offers a rather precise identification of microbial species, but it does not cover the whole microbiome diversity. 16S rDNA NGS is much more exhaustive about the microbiome diversity but is less precise in the identification at species level.
In this study, in order to describe the microbiome evolution in two full-scale biodigesters, managed in the same way, in the different steps of their operation, an integrated approach was followed consisting in (i) ARISA analysis, to preliminarily characterize the bacterial diversity among all the samples collected; (ii) amplification, cloning and sequencing of archaeal 16Srrna and mcrA genes in one sample for each digester; (iii) 16S profiling with the NGS approach to describe the relative abundance of phyla/genera of Bacteria and Archaea. The anaerobic biodigesters described (P1 and P2) are operating in two sites in Lombardia (Italy) separated by 200 km distance. Both plants were fed with cattle sewage and corn silages derived from the respective local farms and operated in mesophilic regimen. P1 was followed from the startup phase up to the phase of full operation (included), while P2 was analyzed when already operative at full capacity. As for P1, a sample of the post-digester was included in the analyses.
The results showed that bacterial and archaeal communities in P1 and P2 were different. Despite these differences, there was no difference in the efficiency of biogas production. It was concluded that an effective microbial succession developed, starting from the inoculum, and reached, in both cases, an optimized (although different) equilibrium of functions. The novelty of this work is to have emphasized the importance of core functions rather than core microbiota in digesters Moreover, results evidenced the role of a still poorly characterized methanogenic archaeum (Candidatus Methanofastidiosum) in the process, suggesting that other pathways than those so far identified and well studied are involved in biogas production.

3
International AG, Switzerland); a PT100 sensor was used to measure temperature and the portable pH 3110 (WTW-Xylem Analytics-Germany) for measuring pH. Chemical analyses were done by external laboratories (Accredia certificated) following the standard procedures: Water & Life Lab https:// www. water lifel ab. it/ and Studio Alfa-https:// www. studi oalfa. it/.
Plant P1 was monitored since the start-up phase in the autumn-winter. Four phases were identified ( In the phase (I) the plant was loaded with a total of 1120 t of bovine sewages. In addition, an inoculum of 540 t of a substrate coming from a similar plant, already operative, was added on 28-29 Oct to speed up the anaerobic digestion. In the phase (II) the manure amount added to biodigester was reduced and dispensed more regularly, although not daily (from 11.60 up to 84.20 t day −1 ). From 6 Dec the average daily amount of sewage added was 20.24 t.
The corn biomass was regularly loaded into biodigester from 7 Nov, starting with 1 t , then gradually increasing until 15 Dec when 44.50 t were loaded. The daily load of plant biomass was then regulated around 39-40 t day −1 . In addition, in the phase (III) a vertical axis mixer was used to integrate two smaller submerged mixers. Finally, in the phase (IV), the primary digester was integrated with a supporting post digester where pig manure was added. The whole process was daily monitored for temperature, pH, VFAs, biogas production and energy cogeneration.
The P2 biodigester was started on 15 March (day 0), and became fully operative by 13 April (day 31). Since then, the average load of cattle sewage was 37.84 t/day with corn silage varying from 19.90 to 49.20 t day −1 The collection of samples started 8 May (day 55).
Five aliquots of digestate were collected in both plants, at each specific point and stored at -20 °C until DNA extraction. In Table 1 are reported the data concerning VFA production (mg/L), T (°C), pH, Alkalinity (mgL −1 ), ratio VFA/ Alk, Ammonia nitrogen concentration (mgKg -1 ), the amount of corn silage (t day −1 ) and manure (t day −1 ) added to the biodigesters, biogas production (m 3 day −1 ), energy production (MWhday −1 ) and are highlighted the samples used for microbial analysis.
All the samples were analyzed with ARISA (see below) and only some representative samples were selected for clone libraries construction and 16S rDNA NGS profiling (see Table 1).
Samples were centrifuged at 1400×g for 10 min, supernatant was discarded and DNA was extracted from the pellets with FastDNA SPIN Kit for soil and FasPrep®Instrument (MP Biomedicals, Santa Ana, CA). DNA samples were checked with electrophoresis on agarose gel 0.9%.

PCR Conditions
For the amplifications of bacterial ITS, archaeal 16Srrna and mcrA and for checking inserts in recombinant plasmids, PCR reactions were conducted with primers pairs listed in Table 2, in a mix of 25 µl containing 5 µl of 1:1000 diluted DNA, 1X PCR Buffer, 200 µM dNTPs, 100nM primers and 0.8 U µl −1 of Taq. A thermal cycler TECHNE TC-512 (Bibby Scientific Ltd, Staffordshire, UK) was used as follows: hotstart 94 °C for 3 min, 35 cycles of denaturation (94 °C for 45 s), annealing (58 °C for 1 min) and extension (72 °C for 2 min) followed by a final extension (72 °C for 7 min).
The primer ReubR, for the amplification of bacterial the ITSReub, was 5′-6 FAM labelled to detect amplicons through capillary electrophoresis. The specific annealing temperature for each primer couple are reported in Table 2.

ARISA (Automated Ribosomal Intergenic Spacer Analysis)
Amplicons obtained from primers ITSF/ITSReub (Table 2) were purified with Wizard®SV Gel and PCR Clean-Up System (Promega, Madison WI, US). An aliquot of 4 µl of purified DNA was mixed with 9 µl of Hi-Di TM formammide (Applied Biosystems, Foster City, CA, US) and 1 µl of GeneScan-2500 ROX TM Size Standard (Applied Biosystems). After denaturation at 95 °C for 3 min and fast cooling, DNA was injected in an ABI Prism 310 (Applied Biosystems, CA, US) capillary filled with POP-4™ Performance Optimized Polymer for electrophoresis. ARISA electropherograms were analyzed using the GeneScan 3.1 software program and elaborated with Microsoft Excel. The raw data were processed excluding the peaks whose heights  and areas, proportional to the intensity of the fluorescence, were less than 1% of the sum of the same. Then the peaks were grouped according to differences of 4 bp for the fragments up to 700 bp and 5 bp for those up to a length of 1000 bp [24]. For each sample the presence (1) or absence (0) of a certain peak was determined obtaining a matrix utilised for calculation of Jaccard similarity coefficient and construction of dendrograms using the programme DendroUPGMA, available at the site http:// genom es. urv. cat/ UPGMA.

Amplification and Cloning of 16Srrna and mcrA (methyl coenzyme M reductase A) Genes from the Archaeal Community
The DNA extracted from selected biodigesters samples (see Table 1) was diluted 1:100 and amplified with the primers, specific for Archaea, reported in

ARDRA (Amplified Ribosomal DNA Restriction Analysis) and Sequencing
In order to select, for sequencing, the recombinant plasmids containing fragment of 16Srrna and mcrA genes from different Archaea, the ARDRA analysis was performed respectively on 100 and 50 clones for each clone library. Restriction analysis of the cloned inserts was conducted with the enzyme HaeIII for 16Srrn, and with enzymes RsaI and AluI for mcrA. 5 µl of amplified DNA was digested with 5 U of each enzyme and 1X of the respective buffer at 37 °C for 4 h. The restriction fragments were loaded on 2.5% NuSieve ® low melting agarose gel (Lonza Group Ltd, Basel, Switzerland) on TBE buffer 0.5% (45 mM Tris borate, 1 mM EDTA pH 8.0) and subjected to electrophoresis at 60 V for 4 h. Images were acquired with Biorad Gel Doc 2000 and the software Quantity One (Bio-Rad, Hercules, CA, US).
One representative clone for the most represented ARDRA profiles was selected for the partial sequencing of 16Srrna gene and mcrA at BMR Genomics srl in Padova (Italy). The obtained sequences were then compared with those stored in the GenBank database at the NCBI (National Center for Biotechnology Information) by using the BLAST (Basic Local Alignment Search Tool) program (http:// www. ncbi. nlm. nih. gov/ blast). The sequences were deposited in GenBank (https:// www. ncbi. nlm. nih. gov/ genba nk/) and their accession number are reported in Table 3 Biodigesters Microbiome Profiling by 16S rRNA Gene Amplification and NGS Sequencing The bacterial and archaeal community profiles of the samples were generated by NGS technologies at the Genprobio S.r.l. Laboratory in Parma (Italy). Partial 16S rRNA gene sequences of duplicate samples were obtained from the extracted DNA by Polymerase Chain Reaction (PCR), using the primer pair Probio_Uni and Probio_Rev, targeting the V3 region of the 16S rRNA gene sequence [25] for bacterial community and the primers pair ArchV56 [26] for archaeal community. Amplifications were carried out using a Verity Thermocycler (Applied Biosystems) and PCR products were purified by the magnetic Agencourt AMPure XP DNA beads (Beckman Coulter Genomics GmbH, Bernried, Germany) in order to remove primer dimers. Sequencing was performed using an Illumina MiSeq sequencer with MiSeq Reagent Kit v3 chemicals. To avoid biases caused by different bioinformatic analysis pipelines, the sequence read pools of each dataset were filtered and analyzed through the same custom script based on the QIIME software suite [27,28]. Paired-end read pairs were assembled to reconstruct the complete amplicons. Quality control retained sequences between 140 and 400 bp and mean sequence quality score of >20. The sequences with homopolymers of >7 bp and mismatched primers were removed to improve quality Table 2 Primers utilized in this study. "Y", "R" "M" and "W" mean degenerate bases *T a = annealing temperature of each primer couple Primer Sequence *T a (°C) Refs. 5′-GTC GTA ACA AGG TAG CCG TA-3′  5′-GCC AAG GCA TCC ACC-3′ 58 [23] 69F ARC934R 5′-YGA YTA AGC CAT GCR AGT -3′  5′-TGC TCC CCC GCC AAT TCC T-3′ 47 [15] SP6 T7 5′-TAT TTA GGT GAC ACT ATA G-3'  5′-TAA TAC GAC TCA CTA TAG GG-3′   50   mcrAF  mcrAR   5′-GGT GGT GTMGGA TTC ACA RTA TGCW-3′  5′-TTC ATT GCR TAG TTW GGR TAGTT-3′ 55 [19] 1 3

ITSF ITSReub
issues (even though Illumina sequencing is not known to be affected by homopolymer sequences). In order to calculate downstream diversity measures, 16S rDNA amplicon sequence variants (ASVs) were defined at 100% sequence homology using DADA2 [29] and singletons were removed. All reads were classified to the lowest possible taxonomic rank using QIIME2 and a reference dataset from the SILVA database v132 [30]. The biodiversity of the samples (alpha-diversity) was calculated with Chao1 and Shannon indices. The 16S rDNA reads were sent to Gen-Bank to obtain their under-accession numbers. They are available as bioproject PRJNA534196.

Results and Discussion
This study was aimed to the analysis of microbiome in two full-scale biogas producing plants (P1 and P2), located in two different sites in North Italy, associated with a 999 kW cogeneration unit, both fed with bovine manure and corn silage, and analyzed in two different seasons. The data discussed below, highlight that these two highly similar plants differed in the archaeal microbiome and, even at less extent, in bacterial microbiome, but did not differ in biogas and energy productivity and other performance and operating parameter.

Performance of the Reactors
In the plant P1, the biogas production started during phase II (Fig1), after the addition of corn silage and of an inoculum coming from a similar, fully operating, plant. As shown in Table 1 At the end of December, 2 month and 20 days from startup, P1 was supported with a secondary biodigester that operated in the same condition of pH and temperature, but was fed with pig slurry instead of cattle slurry. The VFA concentration detected in the secondary biodigester, in 29 Dec (day 10 of activity of the secondary plant), was lower than in primary digester (1152 mgL −1 ). This can be explained by the dilution effect or considering that cows and pigs have different diets and the dry matter content of pigs zootechnical effluent varies from 3 to 6% (versus 8-15% of cows), and has a faster biological degradability due to the lower amount of fiber. While the P1 digester was followed from the startup phase and monitored up to the stabilization phase, the sampling in P2 biodigester started when the production of biogas had already stabilized. In P2 the biogas production increased 40 days from the startup phase (similarly to P1) and reached 8662 m 3 by 22 May (day 68 from startup) with an average value of 7796 m 3 day −1 (comparable with P1 biodigester). The pH values were similar to those detected in P1 (in the range of 7.3 to 7.5) while the average temperature was higher (44.2°C vs 42.8°C). Finally the concentration of VFA was highly comparable to that of P1, with an average value of 1634 mgL −1 and, as in P1, the higher production of biogas (8662 m 3 in day 68) corresponded to the lower concentration of VFA (1379 mgL −1 ) detected.
In both plant the daily energy production (reported in Table 1) was the same in the full operational phase (23.90-24 MWhday −1 ). For a synoptic view of the VFA and biogas production in the two plants is possible to look at the graphical abstract.
In general, from the data reported in Table 1, it is possible to state that the performances of the two plants, identical in structure, but located in different sites therefore processing sewage from different local farms, are comparable.

Bacterial Biodiversity
A preliminary overview of the bacterial communities in digesters P1 and P2 was obtained with a method widely applied for a rapid analysis of biodiversity in large environmental samples, the Automated rRNA Intergenic Spacer Analysis (ARISA). All the samples collected in both biodigesters were analyzed with ARISA then, based on the results of this first screening, few samples were selected and a deeper bacterial community profiling was performed by NGS of 16S rDNA.

Bacterial Community Analysis by ARISA
With the DNA extracted from all the samples collected in P1 and P2 digesters and reported in Table 1, ARISA profiles were obtained and compared by using a pairwise similarity coefficient, the Jaccard index. The dendrogram resulting from analysis with the programme DendroUPGMA showed that bacterial microbiome of P1 and P2 clustered in two distinct groups and, for each plant, samples were similarly grouped (Fig 2). In P1 digester one cluster encompassed samples CM2 and CM3, collected when the VFA production was almost stabilized, but the biogas production was still low (see Table 1); in the second cluster the communities of CM4, CM5, CM6 and CM7, collected in the phase of temperature stabilization and steady biogas production, were more similar and separated from CM8 when the biogas production decreased and from CMP (the post digester). The separate clustering of CMP sample was expected as the post digester was fed with pig manure and had lower VFA concentration. The lowering of VFA concentration can be considered as the result of dilution, and the similarity of microbial communities with those of sample CM8 can be explained with the massive inoculum from the primary digester into the secondary digester (see Fig 1). Similarly, in P2, the samples clustered in two main groups, one encompassing VL1 to VL6, the other encompassing VL8, VL10 and VL12, characterized by a lower VFA concentration and higher temperatures. The most important result of this analysis was the clear difference found between the communities of the two biodigesters, despite being managed in the same way.
ARISA is a method widely applied for a rapid analysis of biodiversity in environmental samples and allows a general description of communities [21][22][23], but it doesn't describe the community in terms of bacterial genera or species. Therefore, in order to have more detailed taxonomic information, the bacterial community was profiled by NGS of 16S rDNA.

Bacterial Community 16S Profiling by NGS
The bacterial community profile by NGS of 16S rDNA was determined on selected samples, in duplicate. Samples CM3, CM8 and CMP were analyzed in P1, while, in P2, the profiling was performed on sample VL5. MiSeq runs produced an average of 75566 sequences for the samples collected in P1 and 83931 for those collected in P2. The rarefaction curves obtained with Chao 1 and Shannon indices (a measure used to estimate the alpha diversity in samples and evaluate whether sequencing efforts allowed to capture all the microbial diversity) highlighted the lowest index of diversity in VL5. In Fig 3 (panel a) are shown only results with Shannon index. The most represented phyla in both P1 and P2 biodigesters were the Firmicutes, with similar percentages in all samples (from 62.64 to 74.54%), followed by Bacteroidetes (14.55 to 19.43%) and Cloacimonetes (from 1.57 to 8.72%) (Fig 4). The members of these phyla have a wide metabolic capacity, including degradation of lignocellulose of corn silage and residual of bovine digestion [31]. Firmicutes are often present in fermentative plants, and in some cases, they are the most abundant bacteria [32][33][34][35][36]. Indeed, they possess the ability to metabolize a large variety of macromolecules into acetate and butyrate, VFA intermediate, and therefore they play an important role in the 1 3 hydrolysis of primary substrates in anaerobic digesters. The higher relative abundance of the phylum Firmicutes in both reactors is most likely related to the very similar inoculum and abundance of lignocellulosic substrates.
Among the Firmicutes phylum, the most represented class was that of Clostridia (between 50.50 and 61.98%), that includes the uncultured group MBA03 (10.79 to 21.27%), the genus Sedimentibacter (7.93 to 13.94%), and Clostridium sensu stricto, present at 12.29% in CM3 and in lower percentages (< 2%) in the other samples probably indicating a specific role of this bacterium in the first phases of the process when alkalinity is higher (14325 mgL −1 ) and biogas production is still low (Fig 5). The main differences between the plants P1 and P2 concern the very low presence of Clostridium sensu strictu, in sample VL5 of P2 (1.16%) where Hydrogenispora (a fermenting bacterium producing acetate, ethanol and hydrogen) represents 8.93% (vs < 1.5% in all samples of P1). More differences between P1 and P2 were detected for the Bacteroidetes phylum: the uncultured bacterium of Lentimicrobiaceae family (9.47% in P2 vs <2% in all samples of P1), Proteiniphilum (5.29% in P2 vs <1.5% in all samples of P1) and the uncultured bacterium belonging to the Bacteroidales family UCG-001 (3.10-3.58% in P1 samples vs < 0.28% in P2). Moreover, P1 was characterized by an increase in the CM8 and CMP samples of the phylum Cloacimonetes, very low in P2. Some authors suggest that members of this phylum are involved in lignocellulose degradation and play a role in the syntrophic oxidation of propionate [37,38]. The genome reconstruction of two candidates of Cloacimonetes, Candidatus Cloacimonas acidaminovorans [37] and Candidatus Syntrophosphaera thermopropionivorans [38] indicate that they possess the genes of the methylmalonylCoA pathway responsible of the oxidation of propionate. The oxidation of propionate and other short-chain fatty acids is a key step for methanogenesis to occur, as it produces acetate and H 2 utilizable by both acetoclastic and hydrogenotrophic methanogens.

Archaeal Biodiversity
Similarly to the analysis of bacterial communities, the first approach for archaeal community characterization was carried out with a more classic molecular method, the cloning and sequencing of amplified target genes (16Srrna and mcrA genes). Afterwards, on selected samples, 16S profiling by NGS was performed.

Archaeal Community Analysis by Cloning and Sequencing of 16S rRNA Genes and mcrA Genes
The sequences encoding partial 16S rRNA (about 870 bp) and part of the gene encoding methyl CoM reductase alpha subunit (mcrA) gene (about 500 bp) were amplified from the DNA extracted from samples CM8 (in P1) and VL10 (in P2) and cloned into the vector pGEM ® -T Easy. The cloned fragment (100 clones/gene library) were analyzed by ARDRA and generated eight reproducible profiles (A-E for P1 and F-I for P2) ( Table 3). Sequence analysis of the phylotype A, representing 36% of the analyzed clones, showed 98.20% similarity with Methanotrix soehngenii (Methanosaeta) an acetoclastic archaeum. Other methanogenic species belonging to genus Methanoculleus (hydrogenotrophic and similar to Methanobacterium) as well as uncultured bacterial species were identified in P1 (Table 3) while in P2 the dominant archaeal genera were the acetoclastic, but more versatile, Methanosarcina (98.80% similarity with the specie thermophila) and Methanoculleus. Within these samples, several Clostridiales and other uncultured bacteria were found. This may be due to the utilization of degenerate primers amplifying also bacterial sequences very abundant in biodigesters.

Archaeal Community 16S Profiling by NGS
The MiSeq runs produced an average of 10912 sequences for the samples collected in P1 and 11184 for those collected in P2. Rarefaction curves obtained with Chao1 and Shannon biodiversity indices highlighted, as for bacteria, the lowest index of diversity in sample VL5 (P2) and the highest in sample CM3 (P1). In Fig 3 (panel b) are shown only the rarefaction curves obtained with Shannon index, but identical results were obtained with Chao1 biodiversity index (data not shown).  This means that, although P1 and P2 had a similar trend in biogas and VFA production and operated in the same conditions of pH and temperature, they shaped in different way the metabolic fermentative pathways. It is worth mentioning that Methanosarcina and Methanosaeta are both acetoclastic methanogens, but operate, respectively, at high and low concentration of acetate [39]. In addition, differently from members of Methanosaeta genus, that are strictly acetoclastic, most Methanosarcina are mixotrophic, as they utilize not only acetate, but also hydrogen derived from acetate oxidation conducted by anaerobic bacteria (in general Clostridia) and are able to produce methane even from methanol and methylamines. Archaea belonging to Methanosarcina genus are able, themselves, of acetate oxidation, thus bypassing the bacterial action [40]. Differently from Methanosarcina, Methanosaeta are rather sensitive to VFA and ammonia but, in this case, the VFA concentration did not significantly differ in P1 and P2 and ammonium concentration, even thought different in P1 and P2, never exceeded 2000 mgL −1 . Therefore, the difference may be related to environmental condition of manure storage, considering that P1 started in autumnwinter, while P2 operates in spring-summer. Powell et al. (2008) observed that ammonia production and emission were significantly higher in spring than in fall-winter, not only because the temperature increases the evaporation, but also because the excreted N in urine is higher and urease activity is stimulated, producing ammonia [41,42]. In this case the concentration of ammonium was lower in the spring/summer biodigester (P2) were Methanosarcina was the most abundant methanogen.
Karakashev et al. [40] suggest that, in the absence of Methanosaeta, the hydrogenotrophic methanogenic pathway dominated, also in presence of Methanosarcina that shifted their metabolism from acetoclastic to hydrogenotrophic methanogenesis. The same authors showed that acetate oxidation is reduced in presence of Methanosaeta. The P1 microbiome showed an evolution characterized by a decrease of Methanosarcina (from 8.38 of CM3 to 0.73 of CM8) and a parallel increase of Candidatus Methanofastidiosum, (formerly called WSA2 group) from 10.30% in sample CM3 to 32.70% in CM8 (and 24.70% in CMP) while Methanosaeta remained substantially unchanged (16,12% in CM3, 23.05% in CM8 and 18.15%in CMP). In P2 the presence of Ca. Methanofastidiosum is low (0.44%) and that of Methanosaeta was negligible (0.04%).
The genome analysis of Ca. Methanofastidiosum by Nobu et al. [43] revealed that this group has the peculiar capacity of producing methane through methylated thiol reduction. Moreover, members of this group may utilize acetate (as well as malonate or propionate) with CO 2 as carbon source. It is possible that, in this way, they maintain low the acetate concentration favoring the growth of Methanosaeta. Our data showed a certain competition between Methanosarcina and Methanosaeta, but also between Methanosarcina and Ca. Methanofastidiosum. In general, Methanosarcina may prevail over Methanosaeta because of the greater tolerance to NH 4 + , acetate and the metabolic versatility; but it is also possible that Ca. Methanofastidiosum compete with Methanosarcina, favoring the development of Methanosaeta. The group of WAS2 (Ca. Methanofastidiosum) was detected in different environments, like freshwater and marine sediments, contaminated groundwater and sludge biodigester [43,44]. Recently this genus has been detected also in municipal sludge biodigester [45] and in bioreactors fed with pig manure [46] but, being uncultivable microorganisms, their metabolic features were deduced only by genome analysis. However, a possible competition of WAS2 with Methanosarcinales for acetate has been hypothesized by Rivière et al. [5] and this could explain its increase in P1 when Methanosarcina decreases.
While Methanosarcina, Methanosaeta and Ca. Methanofastidiosum are acetoclastic methanogenic Archaea, Methanobrevibacter is a typical microorganism of rumen and manure, operating a hydrogenotrophic metabolism [47]. It was detected in P1, in the initial phases of biodigester activity, in sample CM3 (at 10.80%), but it was strongly reduced in CM8 (1.45%), where it was replaced by Methanosaeta (23.05%). It was almost absent in P2 where Methanosarcina predominated in sample VL5. A similar shifting of Methanobrevibacter to Methanosarcina was observed in another research by Ciotola et al. [48]. Nevertheless, the global percentage of hydrogenotrophic Archaea remained substantially unchanged, in CM3 and CM8, as the reduction of Methanobrevibacter has been compensated by an increase of Methanobacterium and Methanoculleus (both hydrogenotrophic). The percentage of hydrogenotrophic Archaea (around 30%) was similar also in CMP and VL5 where, however, Methanobacterium was the prevalent hydrogenotrophic microorganism. The presence of acetoclastic Archaea remained constant in all samples, being the increase of Methanosaeta compensated by the decrease of Methanosarcina.
The NGS analysis of archaeal community indicated a trend similar to that of ARISA analysis for bacteria, with less differences between CM8 and CMP were a predominance of Candidatus Methanofastidiosum, Methanosaeta and about 30% of hydrogenotrophic Archaea were detected. Is rather surprising the increase of Methanosarcina to 4.62% in comparison to 0.72 % of CM8, and a slight reduction of Methanosaeta. These results are somewhat in contrast with data reported in the literature, but according to the observation that low concentrations of VFA favor the development of Methanosaeta and not of Methanosarcina [49]. The slight reduction of Candidatus Methanofastidiosum, seems to confirm its competition with Methanosarcina. Another peculiarity of CMP is the presence of Candidatus Methananoplasma, a group of obligated hydrogen-dependent methylotrophs, but being Candidatus Methanoplasma termitium the only deeply characterised species of this group [50] it is difficult to infer its role in CMP and in methanogenesis in general.

Conclusions
In this study an integrated molecular approach was applied to analyze and compare the bacterial and archaeal communities of two full-scale Anaerobic Digestors (7200 m 3 reaction volume). The two plants, identical in structure, are used for production of biogas and energy from corn silage and bovine manure and are located in Po Valley at 200 km distance. The plants are managed from the same company, in a comparable way. The molecular analyses were conducted from the phase of start-up to the phase of full operation in autumn-winter for P1, while P2 was analyzed when already in full activity during the spring.
It was possible to observe that, in the plant P1, there was, along four phases, a progressive stabilization of the biogas production, that reached values comparable to that of P2. This was accompanied by a reduction of biodiversity of bacterial and archaeal communities, thus reflecting a competition and selection of microbial populations and microbial functions. Although the plant technology, the operating characteristic, the biogas and energy production were comparable in both plants, the microbial communities, in particular those of Archaea, resulted to be notably different. This may be linked to the different seasons in which the two plants were monitored and the samples collected, but also highlights that biogas production efficiency does not necessarily depend on the phylogenetic structure of microbial community acting in AD, but rather on their optimal synergic activity. Indeed, despite the differences observed, the percentage of hydrogenotrophic methanogens remained substantially unchanged in the two plants, suggesting that this pathway is necessary. It was also interesting to observe the relevant presence of Candidatus Methanofastidiosum, an unculturable Archaeum found in extreme environment and, more recently, in some biodigesters, probably involved in methanogenesis through methylated thiol reduction. This underlines that different pathways are relevant in methane production, other than acetoclastic and hydrogenotrophic.
From the data obtained, it is possible to conclude that, in the two biogas plants, two different effective microbial succession developed independently, starting from the inoculum and reached, in both cases, an optimized (although different) equilibrium of the community, ensuring the necessary metabolic functions. Although the phylogenetic composition of the two communities showed important differences it appears that these microbial communities carry out similar functional processes, regardless of differences in their structure (functional similarity). Langer et al. (2015) already reached these conclusions describing the functional redundancy and structural changes of microbial communities in four lab-scale (12 L), continuously stirred tank reactors. The diverse microbial communities optimized their metabolism in a way that ensured efficient biogas production [51]. With the present study, conducted in full scale dimension, the same conclusion can be drawn.