Characterization of Core Microbiomes and Functional Profiles of Mesophilic Anaerobic Digesters Fed With Chlorella vulgaris Green Microalgae and Maize Silage

Microalgal biomass is an alternative feedstock for biogas production although its C/N ratio is usually lower than optimal, therefore co-fermentation is recommended. Identification of the core microbiome by metagenome analysis and prediction of functional characteristics are essential to make microalgal feedstock more sustainable and economically feasible. Biogas production from photoautotrophically grown Chlorella vulgaris (Ch. vulgaris) biomass (240 mL CH4 g oTS-1) and co-fermentation with maize silage (330 mL CH4 g oTS-1) has been studied in semi continuous laboratory biogas fermenters. Maize silage control yielded 310 mL CH4 g oTS-1. The microbial community and the read-based functional profiles, derived from these data, were examined during the process by using next-generation metagenome Ion Torrent sequencing technology. The read-based core microbiome consisted of 92 genera from which 60 abundant taxa were directly associated with the microbial methane producing food chain. The data-set was also analyzed in a genome-based approach. 65 bins were assembled, 52 of them belonged in the core biogas producing genera identified by the read-based metagenomes. The read-based and genome-based approaches complemented and verified each other. The functional profiles indicated a variety of glycoside hydrolases. Substantial rearrangements of the methanogen functions have also been observed. Co-fermentation of algal biomass and plant biomass can be carried out for an extended period of time without process failure. The microbial members of the inoculum are well conserved, feedstock composition cause relative abundance changes in the core microbiome.

Microalgal biomass is a promising feedstock for anaerobic digestion (AD), as it is usually rich in lipids, carbohydrates, and proteins, and does not contain recalcitrant lignin (Yen et al., 2013;Ward et al., 2014). The generated biogas yields are typically 10-20% lower, but the methane content is 7-13% higher than that from maize silage (Mussgnug et al., 2010). Successful AD of algal biomass depends on the efficient deconstruction of the cell wall (Córdova et al., 2018). Physical, mechanical, chemical, and thermal pre-treatment methods have been tested to enhance digestion efficiency (Alzate et al., 2012;Lam and Lee, 2012;Passos et al., 2013Passos et al., , 2014Passos et al., , 2018Lavrič et al., 2017). These methods could improve biogas yield, but the energy input is quite high (Carrere et al., 2016). Hydrolytic enzymes have also been proven effective in microalgal biomass pre-treatment (Demuez et al., 2015;Mahdy et al., 2015), although the economy of the process is frequently jeopardized (Vergara-Fernández et al., 2008). Another tool to make the microalgal biomass a suitable AD feedstock is to adjust the carbon to nitrogen ratio (C/N) to a range between 20 and 30 (Yadvika et al., 2004). Low C/N ratios, i.e., <10-15, which are characteristic of microalgal biomass, result in increased levels of free ammonia, i.e., >0.15 g NH 3 -N L −1 or >3-4 g NH + 4 L −1 , that may become inhibitory (Khalid et al., 2011;Kwietniewska and Tys, 2014). Co-digestion is a straightforward strategy to overcome this obstacle Wirth et al., 2015a,b;Rétfalvi et al., 2016;Wang W. et al., 2016).
The microbiological events taking place in the AD reactors fed with microalgae are not fully understood. Early studies used Denaturing Gradient Gel Electrophoresis (DGGE) and identified Methanosaeta sp. and acetoclastic methanogenesis as the main route to biogas formation (Lakaniemi et al., 2011;Zamalloa et al., 2012). The same conclusion was reached by analyzing the methyl-coenzyme M reductase α-subunit gene (mcrA) by high-throughput amplicon sequencing (Ellis et al., 2012). Amplicon sequencing has limitations arising from possible PCR bias, choice of primers, variable copy numbers in bacterial genomes, sequence variation in closely related taxa, and lack of information about the functions of microbes (Becker et al., 2000;Větrovský and Baldrian, 2013;Oulas et al., 2015;Campanaro et al., 2018b). Development of high-throughput molecular tools and bioinformatics allowed the investigation of total DNA and could provide a more detailed picture of the microbiome. Ion Torrent PGM TM was used to explore the microbial composition of Spirulina AD at alkaline conditions (Nolla-Ardèvol et al., 2015b). A broad range of active membrane proteins were identified, which were responsible for starch degradation, sugar uptake, metabolism of peptides and osmoprotectants. Bacterial community alterations were observed during the co-digestion of microalgae with used cooking oil, mill residue, and maize silage (Rétfalvi et al., 2016). Similar changes were reported in the AD of Scenedesmus obliquus and microalga-bacteria mix (Wirth et al., 2015a,b). Co-digestion of algal biomass with terrestrial plant biomass apparently facilitated the hydrolytic step. Search for core microbiota, i.e., the microbial community present throughout the time span of the process and independently of the feeding regime, has gained increasing attention in studying complex systems in microbial ecology (Nagai et al., 2010;Xia et al., 2011;Jami and Mizrahi, 2012;Schären et al., 2017). Core Operational Taxonomic Units (OTUs), based on 16S rRNA gene amplicon libraries have been reported as being present in diverse anaerobic communities (Rivière et al., 2009;Huse et al., 2012;Rui et al., 2015;Calusinska et al., 2018). Core microbial communities can be more precisely identified by the combination of read-based and binning metagenome techniques (Nolla-Ardèvol et al., 2015a;Campanaro et al., 2016Campanaro et al., , 2018aHeintz-Buschart et al., 2016;Stolze et al., 2016;Treu et al., 2016;Zhu et al., 2018). Metagenomics studies also allow for the extrapolation to ecological functional networks (Schlüter et al., 2008;Wirth et al., 2012;Campanaro et al., 2016;Stolze et al., 2016), although this is of limited information, because the abundances of genes coding for enzymes do not necessarily reflect biological activity (De Vrieze et al., 2016).
We have demonstrated earlier (Wirth et al., 2015a,b) that an algal-bacterial mixed biomass is a suitable substrate for biogas generation. In the present study we report on biogas production from photoautotrophically grown Chlorella vulgaris microalgal biomass. Although microbiologically pure algal biomass is certainly more expensive biogas substrate than an algal-bacterial mixture, it was selected for this study so that the bacterial components, other than the constituents of the microalgal phycosphere, would not distort the core AD microbiome. To our best knowledge, there are no data about the core microbiota digesting microalgal biomass in semi Continuously Stirred Tank Reactor (CSTR) type laboratory biogas fermenters. Therefore, our aim was to determine the response of the core biogas producing microbial community in time to the feedstock, i.e., Chlorella algal biomass. The microbial community and the metagenomics-based functional profile, derived from these data, were monitored during the process by using highthroughput sequencing technology (Ion Torrent PGM TM ). The AD parameters and microbial communities in the anaerobic reactors fed with C. vulgaris, maize silage (internal control in the experiment), or the mixture of the two feedstocks were compared.

Feedstock Composition
Biomass of C. vulgaris was obtained from Lisi (Xi'an) Biotech INC (Shaanxi, Xi'an, China). C. vulgaris had been cultivated in photobioreactors, under natural light illumination at ambient temperature. The protein content was 58.50% and the bacterial count was <100,000 cfu g −1 , according to the specifications of the producer. We have confirmed the cfu value by plating the lyophilized material on complex medium (Luria-Bertani agar). The biomass was stored at 4 • C until utilization. Maize silage was obtained from the mesophilic biogas facility "Zöldforrás Ltd., " located near the city of Szeged, Hungary. The inoculum sludge was originated from the same place. The plant uses pig slurry and ensilaged plant biomass (maize silage and sweet sorghum mixture) as feedstock. The main parameters of the microalgal biomass and maize silage feedstocks are presented in Table 1.

Anaerobic Digestions
Anaerobic digestions were carried out in three of 5 L CSTRs  in fed-batch operational mode. The experimental design and time course followed the scheme described previously (Wirth et al., 2015a,b). In short, the apparatus can be fed continuously or intermittently via a pistontype delivery system, the fermentation effluent is removed through an air-tight overflow. The reactors are equipped with a spiral strip mixing device driven by an electronic engine. An electronically heated jacket surrounds the cylindrical stainless steel body, electrodes for the measurement of pH and redox potential are inserted through the reactor wall, in sealed sockets. The device can be drained at the bottom. The evolved gas leaves the fermentor through the top plate, where ports for gas sampling and the delivery of liquids by means of syringes through silicone rubber septa are also installed. Gas volumes are measured with thermal mass flow devices (DMFC SLA5860S, Brooks). Following the 1 month "start-up" incubation the three reactors were fed as follows: one fermentor received C. vulgaris biomass at a loading rate of 1 g oDM L −1 day −1 (oDM = organic dry matter), one reactor was supplied with a mixture of C. vulgaris + maize silage (0.5 g oDM L −1 day −1 each) and the third one with maize silage (1 g oDM L −1 day −1 ). After the "start-up" phase the accumulated gas volume was recorded in every 4 h. The composition of the evolved biogas was measured with a gas chromatograph (6890N Network GC System, Agilent Technologies) equipped with a 5Å molecular sieve column (length 30 m, I.D. 0.53 megabore, film 25 µm). Ultrapure Ar was used as carrier gas, the detection limit for CH 4 and CO 2 was about 30 ppm. The determination of digestion parameters, i.e., organic dry matter, C to N ratio, NH + 4 -N and volatile organic acids/total inorganic carbon was done as previously described (Wirth et al., 2015a,b).

Organic Acid Analysis
Samples for organic acid analysis were taken from the liquid phase of the reactors. The samples were centrifuged (16,000 g for 10 min,) and the supernatant was filtered through polyethersulfone (PES) centrifugal filter (PES 516-0228, VWR) at 16,000 g for 20 min. The concentrations of volatile organic acids were determined with HPLC (Hitachi LaChrome Elite) equipped with refractive index detector L2490. The separation was performed on an ICSep ICE-COREGEL−64H column. The temperature of the column and detector was 50 and 41 • C, respectively. The eluent was 0.01 M H 2 SO 4 (0.8 mL min −1 ). Acetate, propionate and butyrate were determined (detection range 0.01-10 g L −1 ), the latter two were present in traces relative to acetate and therefore they are not included in the results.

DNA Isolation for Metagenomics Study
The composition of the microbial community was investigated six times during the experimental period, i.e., at the starting point and at the end of weeks 1, 5, 9, 12, and 16 (n = 16). For total community DNA isolation 2 mL of samples were used. DNA extractions were carried out using a slightly modified version of the Zymo Research Fecal DNA kit (D6010, Zymo Research, Irvine, USA). The lysis mixture contained 100 µL of 10% CTAB (cetyltrimethylammonium bromide) to improve the efficiency (Wirth et al., 2015a,b). After lysis (bead beating was performed by Vortex Genie 2, bead size: 0.1 mm; beating time: 15 min, beating speed: max), the Zymo Research kit protocol was followed. The quantity of DNA was determined in a NanoDrop ND-1000 spectrophotometer (NanoDrop Technologies, Wilmington, USA) and a Qubit 2.0 Fluorometer (Life Technologies, Carlsbad, USA). DNA purity was tested by agarose gel electrophoresis and on an Agilent 2200 TapeStation instrument (Agilent Technologies, Santa Clara, USA).  Figure 1). The primary sequences produced by Ion Torrent PGM TM were further analyzed by the Diamond software, applying the default LCA (Lowest Common Ancestor) algorithm (Buchfink et al., 2014). Diamond filtering parameters were set as follows: Blast Mode: BlastX, Min Score: 50, Max Expected: 0,01, Top Percent: 10, Min Support: 1. NCBI nr and EMBL-EBI InterPro databases were used for taxonomic and functional alignment (Federhen, 2015). The data were analyzed by Megan6 (Huson et al., 2007(Huson et al., , 2011(Huson et al., , 2016.

Availability of Data and Material
Sequences are available on MG-RAST under the study name: "Anaerobic digestion of C. vulgaris microalgal biomass" (ID: mgp20184). Data are also available on NCBI Sequence Read Archive (SRA) under the submission number: PRJNA544905.

Statistical Analysis of Metagenomics Data
Megan6 (version 6.8.18) was used to investigate microbial community and export data for statistical calculations. UPGMA (Unweighted Pairwise Grouping Method with Arithmetic mean) with Bray-Curtis method was employed to cluster the samples (Bray and Curtis, 1957). For microbial and functional core calculation the interactive web tool MetaCoMET (Metagenomics Core Microbiome Exploration Tool) was used (Biom formatted data were exported by Megan6). Default parameter sets were fixed at absolute read count threshold: 100, with the persistence Venn diagram type (additional settings: between group relative abundance: 0.8, samples in group relative: 0.8) (Wang Y. et al., 2016). The functional profiles of the reactors, fed with the various feedstocks, were investigated using metagenomics data and EMBL-EBI InterPro database (Schlüter et al., 2008;Wirth et al., 2012;Campanaro et al., 2016Campanaro et al., , 2018aStolze et al., 2016;Finn et al., 2017;Maus et al., 2017).
Statistical Analysis of Metagenomic Profiles (STAMP; version 2.1.3 Parks and Beiko, 2010;Cai et al., 2016) was used to compute the abundance differences in the case of whole microbiome and core taxa (http://kiwi.cs.dal.ca/Software/STAMP). Dissimilar taxa were identified with two-sided White's non-parametric t-test (White et al., 2009) at 0.95 confidence intervals and the results with q-value (corrected p-value) of <0.05 were retained.

Process Parameters During AD
During the anaerobic degradation of the various substrates, temperature (37±0.5 • C), mixing speed (10 rpm min −1 ) and pH (7.5-8.3) were continuously monitored by a software . The redox potential was below −500 mV throughout the experiment. By the end of the "start-up" period the residual biogas production ceased completely. Gas production data were collected during weeks 1 through 16 in every 4 h and the cumulative CH 4 productions are plotted in Supplementary Figure 5. Figure 1A shows the CH 4 yields.
The CH 4 content in the evolved gas from C. vulgaris and from maize silage was 57-67 and 50-52%, respectively, these were comparable with previously reported data (Amon et al., 2007;Mussgnug et al., 2010), although in the case of C. vulgaris we measured somewhat lower CH 4 than in a previous study (Mussgnug et al., 2010). Co-digestion of C. vulgaris and maize silage yielded an intermediate CH 4 value of 57-59%.
The VOAs/TIC (Volatile Organic Acids/Total Inorganic Carbon) ratio is a reliable indicator of a stable digestion process (Lienen et al., 2013). The optimal ratio ranges between 0.1 and 0.5 (Leite et al., 2015). Most of our values fell within this range, due to the low organic loading rate (OLR), indicating a stable and balanced operation in all three reactors, although toward the end of the experimental period the VOAs/TIC in the reactor digesting C. vulgaris alone started to increase ( Figure 1B).
The amount of NH 3 -N, conveniently measured as NH + 4 -N, is also a critical indicator of stability of the biogas-forming process (Yenigün and Demirel, 2013). Levels above 4 g NH + 4 -N L −1 , corresponding to about 0.15 g NH 3 -N L −1 under the current experimental conditions, apparently have a negative effect on the methanogenic community (Chen et al., 2008;Nielsen and Angelidaki, 2008). The inhibitory concentrations reported in earlier studies varied, depending on feedstock, inoculum, environmental conditions, and acclimation (Chen et al., 2008;Yenigün and Demirel, 2013;Gonzales-Fernandez et al., 2018). In the reactor fed with C. vulgaris, the NH + 4 -N content tended to increase, and by week 9, approached, and later exceeded the limiting inhibitory concentration (week 12: 5.16 g NH + 4 -N L −1 or 0.26 g NH 3 -N L −1 ) ( Figure 1C and Supplementary Figure 3). It is noteworthy that co-digestion showed balanced digestion and approached the limiting NH + 4 -N values (3.99 g NH + 4 -N L −1 or 0.15 g NH 3 -N L −1 ) only in week 16. The acetate concentrations ( Figure 1D) exhibited a tendency similar to that of VOAs/TIC ( Figure 1B). The acetate concentration in the C. vulgaris reactor reached the inhibitory limit of 3 g L −1 (Siegert and Banks, 2005) by week 12. Ammonia and acetate accumulation was due to the low C/N ratio of the microalgal biomass, although the pH remained stable in the C. vulgaris reactor because of the high TIC.
The hierarchical clustering indicated good correlation with ammonia and acetate concentrations suggesting that these parameters prominently affected the development of the microbial communities (Supplementary Figure 3). Co-digestion efficiently prevented ammonia and acetate accumulation corroborating earlier observations (Yen and Brune, 2007;Wirth et al., 2015a,b).

The Core Microbial Community
Although we detected substantial alterations in the AD microbial communities fed with the various feedstocks, it is intriguing to see which members of the inoculum community persisted throughout the time span of the experiment in spite of the distinct feeding regimes. A mixed algal-bacterial biomass is the obvious choice of microalgal biomass production for largescale biogas generation due to trivial economic feasibility reasons. It is, however, difficult to separate the bacterial community introduced with a mixed algal-bacterial biomass from the community introduced with the inoculum and adapted to this particular substrate. That is why in this study we used microbiologcally pure C. vulgaris biomass substrate to properly establish the core community participating in the process. The core microbiome exploration tool was employed (MetaCoMET) in search for this core microbial community (Rivière et al., 2009;Rui et al., 2015). In the read-based metagenomics database the "genus" taxonomic rank was chosen to identify the key players because the depth of the sequencing allowed the taxonomic resolution at this level (Supplementary Figure 1). In addition to the read-based evaluation, we analyzed the same sequencing data-set by three binning programs to determine the genome-based taxa involved in the biogas forming process. Due to the high complexity of the AD microbial community and the drawbacks of shortread sequencing technology, only a limited number of taxa can be recovered by binning. The read-based and genome-centric approaches are complementary in providing a more reliable picture of the AD microbial community.
We identified first the sub-cores, i.e., the genera present above threshold abundance (0.8%) at all 5 time points (and in the inoculum) in the reactor fed with each feedstock (Figure 2A). This distinguished 125 genera in the maize silage AD (31% of all genera detected in this reactor), 106 genera in the C. vulgaris AD (26% of all microbes in this reactor), and 113 genera in the co-digestion reactor (23% of the total microbes in this reactor).
The core community overlapping in all reactors, i.e., the collection of genera found in all 16 metagenomes above threshold, comprised differently abundant 92 genera ( Figure 2B and Supplementary Table 2). The composition of the core microbiome is shown in Figure 3. The high number of genera preserved in the core community suggests that the inoculum has a remarkable impact on the evolution of the microbial community in the AD reactor (De Vrieze et al., 2014;Han et al., 2016;Liu et al., 2017;Mahdy et al., 2017;Gonzales-Fernandez et al., 2018).
In the following discussion we combine the read-based metagenomics data with binning results (Figure 4). This information was used to position the various core genera all along the AD process from polymer hydrolysis to biogas evolution. The assembly of the elements of the "core puzzle" is summarized in Figure 5. Sixty genera from the read-based core taxa could be associated with various stages of AD. From the 65 bins 52 could be related with the methane producing microbial food chain. It should be noted that 42 of these bins are listed in Figure 5 because bins identified under the same name but distinguishable as separate bins are listed only once.
It general, it is noteworthy that the read-based sequence analyses are complemented adequately by the genome-based approach. The former established a good picture of the participating microbial community at genus level, while the latter indicated several species identification, but also numerous "uncultured" or partly "unknown" bins. This confirmed that the AD community is only partly understood in details (Campanaro et al., 2018a,b).
The core genera involved in polymer hydrolysis were searched for first (green box in Figure 5). Diversity characterized the order Clostridiales. This is in accordance with previous observations (Wirth et al., 2015a,b) (Supplementary Figure 2). FIGURE 2 | The results of core microbiome calculation. (A) The Venn diagram presents the sub-core habitat associations within the anaerobic digesters. The sub-cores, i.e., genera present at the 6 sampling times are included in the diagram, and in the intersection of the sub-cores representing the core, i.e., genera shared among all samples. (B) Sub-cores: taxa shared in the communities of the specific feedstock digestion at all 6 time points. Other microbes: genera present in some but not in all 6 samples of the given feedstock.   (Figure 4). These bins contain a wide variety of putative glycoside hydrolases and genes coding for hydrogen production (Supplementary Table 4). The possible end products from polysaccharide degradation are volatile fatty acids (VFAs) (Lin et al., 2007;Wirth et al., 2012) (Figure 5).
The genome fragment of Herbinix luporum (Bin 25) is predicted to possess cellulose degrading activity, genes coding for endoglucanase and cellobiohydrolase (Koeck et al., 2015). The genera Clostridia and Herbinix showed the largest variances among Bacteria when the substrates maize silage and C. vulgaris biomasses were compared (Figure 6). These hydrolytic genera were considerably more abundant in the reactors supplied with maize silage and a similar trend was apparent in the codigestion reactor.
Ruminiclostridia are able to degrade cellulose and hemicellulose, and ferment various sugars to VFAs (Ravachol et al., 2015). Ruminiclostridium thermocellum (Bin 27) and an unknown Ruminococcaceae (Bin 2) were identified in the genome-based evaluation. Their membrane sugar transporter and phosphotransferase enzyme coding genes may place them among the hydrolyzing bacteria. Fermentimonas and Lentimicrobium genera are also among the abundant hydrolyzer core microbes (Figure 5). Fermentimonas caenicola (Bin 21) and Lentimicrobium saccharophilum (Bin 8) assembled genome fragments fit into this picture (Figure 4 and Supplementary Table 4).
The core genus, Candidatus Cloacamonas belongs in the WWE1 candidate division (Pelletier et al., 2008) and has been found primarily in AD reactors operating with pig slurry (Kröber et al., 2009;Wirth et al., 2012). This genus predominated by a large margin in our three digesters (Supplementary Table 2,  Figure 3). Automated binning and subsequent manual assembly resulted in 2 medium quality (Bins 11 and 65) and one low quality genome fragment (Bin 12) associated with this group (Bowers et al., 2017) (Figure 4). Tetranucleotide frequency and single copy gene content positively placed the two bins of CC. acidaminovorans separately on the clustering FIGURE 5 | The distribution of the core microbial genera along the biogas producing microbial chain of events. The identified core microbes are arranged according to their known physiological roles in the steps of the anaerobic degradation process. Binning results are applied to make more detailed food chain. For detailed explanation see text. SAOB, syntrophic acetate oxidizing bacteria; SPOB, syntrophic propionate oxidizing bacteria. The genera are listed as their mean abundance order in the boxes.
FIGURE 6 | Significant differences in core microbial communities. Significances were determined with two-sided White's non-parametric t-test at 0.95 confidence.
The core oligosaccharide utilizer and acid producer group was represented by Acholeplasma laidlawii (Bin 26), Lascolabacillus massiliensis (Bin 22), Sphaerochaeta globosa (Bin 37), Unknown Treponema (Bin 52), and Treponema brennaborense (Bin 51). Read-based metagenomics data showed that the genus Streptococcus had low abundance in the digesters fed with microalgae as a monosubstrate (Figures 5, 6). Unknown Streptococcus (Bin 34) was identified as corresponding genome fragment in genome-based analysis (Figure 4 and Supplementary Table 4). The minor saccharolytic genera Dysgonomonas and Alistipes were more abundant in the C. vulgaris fed reactor. From the members of the genus Alistipes (read-based), Alistipes sp. CAG:831 (Bin 53) was detected by the binning approach. Bins 48, 50, and 59 could be identified as unknown Candidatus zambryskibacteria, Acidaminococcus sp. CAG:917, and Anaerolinae bacterium 49_20 were additional oligosaccharide degraders. Unknown Porphyromonadaceae (Bin 20), unknown Clostridia (Bin 18), Peptococcaceae bacterium sp. 1109 (Bin 41) contained formyltetrahydrofolate synthetase (FTHFS), the indicator enzyme of anaerobic acetate oxidation, therefore these are possible Syntrophic Acetate Oxidizing Bacteria (SAOB) (gray box in Figure 5) (see section Metagenome Based Functional Changes). The genera Syntrophaceticus and Mesotoga are likely to have SAOB function, although bins representing them were not found, therefore this finding is not supported by the genome-based evaluation of data (Westerholm et al., 2010;Wang et al., 2013;Nobu et al., 2015;Müller et al., 2016). Syntrophic Propionate Oxidizing Bacteria (SPOB) were observed in all three reactors at all sampling time points by read-based analysis (Syntrophomonas, Pelotomaculum, and Smithella). These genera are able to convert propionate to acetate, H 2 and CO 2 (Lueders et al., 2004;Kosaka et al., 2008) in syntrophy with hydrogenotrophic methanogens (blue box in Figure 5). Indeed, Syntrophomonas wolfei (Bin 33) came up in the genome-based evaluation as a well-known SPOB bacterium (Figure 4) (McInerney et al., 1981;Popp et al., 2017). The low abundance of the genera Pelotomaculum and Smithella may explain why these taxa are missing from the binning results.
Similarly, low abundance characterized the sulfate reducing bacteria Desulfotomaculum, Desulfosporosinus, and Desulfitobacterium (Friedrich et al., 1996;Pester et al., 2012) in the read-based metagenomes. These microbes may represent an additional H 2 sink in AD. They were also present in all reactors examined and were part of the core microbiome as permanent rivals of reductive homoacetogens and hydrogenotrophic methanogens for H 2 (Liu and Lu, 2018) (Supplementary Table 2). Nevertheless, sulfate reducers were not found in the genome-based evaluation.
It is noteworthy that the following precautions have to be kept in mind when the above characterization of the core microbiome is interpreted. (1) The 92-member core genera seem sufficient to cover all basic metabolic steps necessary for biogas formation from the green alga and maize silage. This does not mean, however, that this diversity is the minimum requirement for successful degradation of organic matter and concomitant biogas production (Calusinska et al., 2018;Campanaro et al., 2018a). (2) Here we could reach the genus level resolution of the microbiota by the read-based metagenomes. Some genera have fairly diverse species composition leading to diverse properties of the taxon. The analysis of the community by genome-based evaluation may help to make more precise understanding of the abundant taxa. Nevertheless, depending on the sequence number and the observed microbiota complexity only the relatively abundant sequences can be compiled in MAGs (metagenome-assembled genomes) or quality genome fragments (Bowers et al., 2017). (3) The information content of the available sequence databases is limited, a portion of the sequenced but not classified reads may belong to "microbial dark matter, " e.g., "unknown" bins such as Bin 5, 15, 17, 29, 31, 36, 46, 57 (Hatfull, 2015;Nobu et al., 2015;Solden et al., 2016).

Metagenome Based Functional Changes
The metagenomics data were also used to extract potential functional information about the biological activity of the AD communities fed with the various feedstock combinations. Genes coding for characteristic and signature proteins have been detected and their relative copy number in the samples was determined by counting the reads related to the gene in question. This information was corroborated by the binning study, which helped to associate specific functions to taxa. Although the presence of a certain gene and the number of reads linked to it does not necessary mean that the gene is being transcribed and translated, it is safe to assume at least a limited correlation between these parameters (Wirth et al., 2018). The main functional profiles investigated were the carbohydrate-active enzymes (CAZy), protein degradation, syntrophic acetate oxidation, i.e., formate-tetrahydrofolate ligase, hydrogen metabolism and methanogenesis (Figure 7) (Hattori, 2008;Caspi et al., 2014;Lombard et al., 2014;Rawlings et al., 2018). The glycoside hydrolase families (GH) catalyze the breakdown of poly-and oligosaccharides. GH9 predominated the cellulase group. This enzyme family was found in our unknown Bacteroidetes (Bin 13), unknown Porphyromonadaceae (Bin 20), Herbinix luporum (Bin 25), unknown Bacteroides (Bin 43), and in unknown Bacteroidia (Bin 44) assembled genome fragments (Supplementary Table 4). GH9 has been identified in numerous biogas producing ecosystems (Wei et al., 2015;Güllert et al., 2016;Maus et al., 2016). The GH8, GH9, and GH39 families were detected among the core functions (Supplementary Table 3, green color). As a tendency, the GH genes coding for cellulases were more abundant in the reactors running on maize silage or co-digestion than in the C. vulgaris fed one.
The GH43 enzyme family emerged from the hemicellulases and was found in the unknown Porphyromonadaceae (Bin 20), Fermentimonas caenicola (Bin 21), Lascolabacillus massiliensis (Bin 22), Herbinix luporum (Bin 25), Peptococcaceae bacterium 1109 (Bins 41 and 42), unknown Bacteroides (Bin 43), and Unknown Prevotella (Bin 45). The GH43 family was identified as one of the 15 most abundant GH families in a thermophilic biogas plant . The richness of cellulases and hemicellulases were reduced in the reactors fed with C. vulgaris biomass. The large number of oligosaccharide-degrading GH enzymes covers a wide range of catalytic activity ( Figure 7A). GH2, GH4, GH29, GH31, GH42, and GH77 were present in all sequenced samples although none of them showed significant (p < 0.02) differences among the various reactor functional profiles. The GH2 and GH31 families were detected in the bins of the unknown Bacteroidetes  Table 4). It should be noted that not all the GH families presented in Figure 7A are part of the core functions, which are compiled in Supplementary Table 3, highlighted in green color. In general, maize silage and co-digestion appeared to develop similar CAZY profiles.
Based on InterPro data and the EMBL-EBI MERPOS peptidase database classification, a large number of endo-and exopeptidases were present in our digesters (Supplementary Table 3, highlighted in yellow color) (Rawlings et al., 2018). The great variety of peptidases identified by readbased metagenomics was confirmed by binning results. The genes involved in protein degradation were higher in representation in the reactors digesting C. vulgaris than those fed with maize silage (Figure 7B). The known functions of the peptidases in family S51 are nutritional; lack of these enzymes leads to incomplete degradation of intracellular proteins to amino acids (Lassy and Miller, 2000). This peptidase was found in the unknown Bacteroidetes (Bin 6) and unknown Porphyromonadaceae (Bin 20) assembled genomes. The peptidase family M20 was present in the largest abundance among functions (Supplementary Table 3 Formate-tetrahydrofolate ligase (i.e., formyltetrahydrofolate synthetase) (FTHFS) was also among the selected signature enzymes. FTHFS is a highly expressed key enzyme in both the Wood-Ljungdahl pathway of autotrophic CO 2 fixation (acetogenesis) and the glycine synthase/reductase pathway of purinolysis. The key physiological role of this enzyme in acetogens is to catalyze the formylation of tetrahydrofolate, an initial step in the reduction of CO 2 and other onecarbon precursors to acetate. In purinolytic organisms, the enzymatic reaction is reversed, i.e., it releases formate from 10-formyltetrahydrofolate with concurrent production of ATP (http://pfam.xfam.org/family/PF01268). The detailed biochemical and regulatory mechanisms behind the shift and the way in which the bacteria regain energy is uncertain in spite of intensive investigations (Pester and Brune, 2006;Ragsdale and Pierce, 2009;Xu et al., 2009). FTHFS was present among the putative core functions (Supplementary Table 3, gray color) and its relative abundance was somewhat elevated in the reactors digesting C. vulgaris relative to the ones fed with maize silage and/or co-digestion ( Figure 7C). The FTHFS gene was present in unknown Clostridia (Bin 18), unknown Porphyromonadaceae (Bin 20), Peptococcaceae bacterium 1109 (Bin 41), and in Candidatus Cloacimonas acetaminovorans (Bin 65) assemblies.
The equilibrium between H 2 producers and consumers is crucial in methanogenesis (Wirth et al., 2012;Szuhaj et al., 2016;Bagi et al., 2017). The relative numbers of identified hydrogenase genes were nearly equal in the three reactors ( Figure 7D) (Supplementary Table 3, purple color). The gene coding for [FeFe] hydrogenase maturation, HydE was among the predicted core functions. The HydE protein, together with HydF and HydG comprise the accessory gene products for the assembly of an active [FeFe] hydrogenase (King et al., 2006). Genes coding for enzymes involved in iron hydrogenase maturation were found in the metagenome assembled genomes of unknown Bacteroidetes (Bin 6 and 13), Herbinix luporum (Bin 25), Ruminiclostridium thermocellum (Bin 27), and in both Candidatus Cloacimonas acetaminovorans (Bins 11 and 65). The observation suggests an equally important role of H 2 metabolism in biogas formation, independent of the feedstock composition.
The genes coding for enzymes involved in the biogas evolution pathways were our final group of signature functions tested. A collection of 24 methanogenic enzymes were identified (Supplementary Table 3, red color); 14 of them were present in the core. Hydrogenotrophic methanogenesis predominated in maize silage and co-digestion. In contrast, acetoclastic methanogenesis was present in higher proportion when only C. vulgaris was the feedstock (Figure 7E). Among the core methanogenic enzymes, methyl coenzyme M reductase (subunit A, B, G, and D) and CoB-CoM heterodisulfide reductase (subunit B) were detected. Methyl coenzyme M reductase genes were observed in unknown Methanoculleus (Bin 39), Mathanosaeta concilii (Bin 47), unknown Methanosarcina (Bin 64) (Supplementary Table 4). The methyl coenzyme M reductase catalyzes the reduction of methyl-coenzyme M and coenzyme B to methane (Luton et al., 2002;Nunoura et al., 2008).

DATA AVAILABILITY STATEMENT
Sequences are available on MG-RAST under the study name: "Anaerobic digestion of C. vulgaris microalgal biomass" (ID: mgp20184). The data are also available in NCBI Sequence Read Archive (SRA) under the submission number: PRJNA544905.

AUTHOR CONTRIBUTIONS
RW with the help from TB and GL, developed the DNA extraction protocol, designed and performed the experiments, and contributed to the evaluation of metagenomic data. GM organized and performed the DNA sequencing work. KK and ZB conceived the project and participated in its design. RW and KK drafted the manuscript. GR critically evaluated the manuscript. All the authors have read and approved the final manuscript.

FUNDING
This study has been supported in part by the Hungarian National Research, Development and Innovation Fund projects GINOP-2.2.1-15-2017-00081, GINOP-2.2.1-15-2017-00033, and EFOP-3.6.2-16-2017-00010. RW, GL, and GM received support from the Hungarian NKFIH fund projects PD121085, PD123965, and FK123899 financed under the PD16 and FK16 funding schemes. This work was also supported by the János Bolyai Research Scholarship (for GM) of the Hungarian Academy of Sciences.

ACKNOWLEDGMENTS
We acknowledge the research funding received to carry out this study.