Variation of Metagenome From Feedstock to Digestate in Full-Scale Biogas Plants

Anaerobic digestion (AD) has been widely used to resolve the problem of organic wastes worldwide. Previous studies showed that the types of feedstock have a great influence on the AD microbiome, and a huge number of AD populations are migrated from upstream feedstocks. However, the changes of microbial compositions from feedstock to AD digestate are still less understood. We collected feedstock samples from 56 full-scale biogas plants, generated 1,716 Gb feedstock metagenomic data in total, and constructed the first comprehensive microbial gene catalog of feedstock containing 25.2 million genes. Our result indicated that the predominant phyla in feedstock are Firmicutes, Bacteroidetes, and Proteobacteria, which is similar to that in AD digestate, and the microbial diversity of feedstock samples is higher than that of AD digestate samples. In addition, the relative abundance of most genes involved in methanogenesis increase from feedstock to AD digestate. Besides, the amount of antibiotic resistance genes (ARGs) and pathogenic bacteria in AD are effectively reduced compared to feedstocks. This study provides a comprehensive microbial gene catalog of feedstock, and deepens the understanding of variation of microbial communities from feedstock to AD digestate of full-scale AD. The results also suggest the potential of AD to reduce the level of ARGs and pathogens in animal manure.

Introduction

The vast amount of global organic waste due to population increase, urbanization expansion, and agriculture intensification has become a huge environmental burden (Tyagi and Lo, 2013). Recently, anaerobic digestion (AD) has been widely used to resolve the problem of organic wastes like animal manure, crop residues, and wastewater sludge and can produce biogas as renewable energy (Luo et al., 2016; Stolze et al., 2016). AD includes four sequential metabolic steps: hydrolysis, acidogenesis, acetogenesis, and methanogenesis (Angenent et al., 2004; Hassa et al., 2018), which are performed by highly complex consortium of microorganisms. In our previous metagenomic study on digestate of full-scale biogas plants (BGPs) fed with diverse feedstocks, we constructed a microbial gene catalog of AD (22,840,185 genes), and showed that the type of feedstock (cattle, chicken, and pig manure) has a great influence on carbohydrate hydrolysis and methanogenesis (Ma et al., 2021).

Feedstock for biogas plants contains a huge number of organic wastes, such as crop residues, animal manure, wastewater sludge, and food waste, and they were temporary stored in a tank before pumped into anaerobic digesters. During storage period, many spontaneous fermentation reactions were performed by microbes. Previous studies had indicated that a huge number of AD populations were migrated from upstream feedstocks in BGPs (Mei et al., 2017). However, the changes of microbial compositions from feedstock to AD digestate were less understood. Besides, since animal manure contains high level of antibiotic resistance genes (ARGs) and human pathogens, the manure directly applied as fertilizer will increase the dissemination of ARGs and pathogens to natural environment, which poses high risk to human health (Youngquist et al., 2016; Wu et al., 2018). Existing studies have proposed AD as a promising method to reduce some ARGs and pathogens in manure. However, whether AD can reduce all types of ARGs and pathogens remain to be verified.

In this study, we performed metagenome sequencing for corresponding feedstock samples from the 56 full-scale BGPs that were used in our previous metagenome study for AD digestate (Ma et al., 2021). These BGPs were operated at different temperatures, fed with diverse feedstocks, and distributed widely in geographical regions. We constructed a microbial gene catalog of feedstock, and compared metagenome of AD digestate and feedstock, to explore the variation of microbial communities from feedstock to digestate, as well as the changes of antibiotic resistance genes and pathogens during AD process.

Materials and Methods

Sample Collection

Feedstock samples were collected from 56 different full-scale biogas plants (BGPs) across 15 provinces of China, with one sample from each BGP. The characteristics of the 56 full-scale biogas plants were already described in our previous study (Ma et al., 2021). Feedstock samples were collected before pumped into anaerobic digesters, and 300 ml of feedstock were sampled from each BGPs and stored in 6 sterile, gastight tubes (50 ml). The samples were frozen immediately in a cooler with dry ice, and then transported to the laboratory. Frozen samples were stored at −80°C before DNA extraction.

DNA Extraction, Library Preparation, and Sequencing

Genomic DNA of feedstock samples was extracted using the same protocol as for AD digestate samples (Ma et al., 2021), and all samples were processed together during the DNA extraction. Frozen feedstock samples were taken out from −80°C refrigerator and thawed at room temperature. Genomic DNA was extracted using PowerSoil DNA Isolation Kit (cat. no. 128880-100; MoBio Laboratories Inc., United States). The entire process was carried out according to the manufacturer’s instructions, with minor modification of adding an additional four freeze-thaw cycles (alternating between 65°C and liquid nitrogen for 5 min) before using the kit. The integrity of extracted DNA was checked on 0.7% (w/v) agarose gel with GelRed nucleic acid gel stain (cat. no. 41003; Biotium, United States). The quality and quantity of the extracted DNA were assessed using NanoDrop (Thermo Fisher Scientific, United States) and Qubit 2.0 Fluorometer (Thermo Fisher Scientific, United States). After DNA quality checks, the high-quality DNA (band length > 15 kb, A260/280 1.8–2.0, dsDNA concentration >20 ng/μL) of each sample were used for library construction. Sequencing libraries with 350 bp insert size were prepared using Illumina TruSeq DNA PCR-Free Library Preparation Kit (ref. 15037059; Illumina, United States) following the manufacturer’s instructions. Sequencing was performed on an Illumina Hiseq X10 platform (Illumina, United States) by Cloud Health Genomics Ltd (Shanghai, China) applying 150 bp paired-end strategy.

Metagenome Assembly and Construction of the Gene Catalog

The metagenome assembly and construction of gene catalog for feedstock samples were performed using the same pipeline for our previous study (Ma et al., 2021). In brief, (1) raw reads were filtered by trimming the adapter contamination and low-quality bases using clean_adapter and clean_lowqual with default parameters, resulting in the clean and high-quality reads with an average error rate <0.001; (2) de novo assembly of the clean reads was carried out using Megahit (v1.1.3) (Li et al., 2016) under paired-end mode, with the preset “meta-sensitive” and the minimum contig length was set to 1,000 bp; (3) the assembled contigs were subjected to a gene prediction using the Prodigal software (v2.6.3) (Hyatt et al., 2012) with parameter “-p meta,” and filter out the genes with codon sequence length less than 100 bp according to a previous study (Qin et al., 2010); (4) the pooled genes from all 56 feedstock samples were then clustered to get the final non-redundant gene catalog using CD-HIT (v4.6.6) (Fu et al., 2012), with parameter “-c 0.95 –n 10 –G 0 –aS 0.9,” adopts the criteria of identity ≥95% and alignment coverage ≥90% of the shorter genes.

Taxonomic and Functional Annotation of the Gene Catalog

Taxonomic and functional annotations of genes in the gene catalog of feedstock were performed using a previously established method (Ma et al., 2021). For taxonomic annotation, CARMA3 (Gerlach and Stoye, 2011) was used on the basis of DIAMOND (v0.8.28.90) (Buchfink et al., 2015) alignment against the NCBI-nr database. Functional annotations were made by DIAMOND searches against the KEGG database (release 79) (Kanehisa et al., 2004) by taking the best hit with an e-value < 1e-5. For the antibiotic resistance genes (ARGs) annotation, genes were predicted using Resistance Gene Identifier (v4.2.2)¹ against the Comprehensive Antibiotic Resistance Database (CARD, v3.0.0) (Jia et al., 2017). The species with their scientific names included in Virulence Factor DataBase (VFDB 2019) (Liu et al., 2018), Global Pathogenic Bacteria Database (GlobalRPH), or CARD v3.0.0 (Jia et al., 2017) were identified as human pathogenic bacteria.

The relative gene abundance was calculated according to a previous study (Huang et al., 2018; Ma et al., 2021). To calculate the relative gene abundance of each sample, the clean reads of each sample were mapped separately onto the gene catalog by BWA-MEM, and the reads with alignment length ≥50 bp and identity >95% were defined as qualified reads. For each sample, total number of reads mapped to all genes (TA) equal to the count of qualified reads, total number of reads mapped to one gene (TO) equal to the count of qualified reads mapped to the gene. Considering the great differences of assembled gene length (GL) and mapped reads number (TA), the value of GL and TA were centralized by dividing 1,000 and 10,000,000 times respectively, since the average GL was about 1,000 bp and the level of TA was 10 million. At last, the normalized gene abundance (NGA) for each sample was calculated as followings: NGA = TO/(GL/1,000)/(TA/10,000,000). The relative abundance of different taxonomic ranks (superkingdom, phylum, class, order, family, genus, and species), KEGG orthologous groups (KOs), and ARGs were calculated by summing the relative abundance of the respective genes belonging to each category.

Analysis of Metagenome-Assembled Genomes

Construction of metagenome-assembled genomes (MAGs) for feedstock samples was performed using the same protocol as for digestate samples (Ma et al., 2021). In brief, the reads of each sample were mapped back to the assembled contigs using BBmap v38.50², and then the genomes were independently recovered from each sample by using the software of MetaBAT2 v2.12.1 (Kang et al., 2019), setting minimum scaffold length of 2,000 bp. Completeness (Cp) and contamination (Ct) of the MAGs were determined using CheckM v1.0.7 (Parks et al., 2015). De-replication of MAGs was performed according to their genome-wide average nucleotide identity (ANI) among different MAGs, and MAGs with ANI value >95% and genome coverage >50% were considered to belong to the same species, then the one with the highest CC3 value (CC3 = Cp – Ct ^∗ 3) was selected as the representative. Finally, 2,101 MAGs were constructed for feedstock samples. Taxonomic annotation of these MAGs was performed with GTDB-Tk v1.3.0 (Chaumeil et al., 2019).

Statistical Analysis

The microbial communities of feedstock were analyzed using the same protocol as published previously for analyzing the AD digestate from the same 56 BGPs (Ma et al., 2021). The relative gene abundance data of AD digestate for these 56 BGPs were cited from our previous study (Ma et al., 2021), and the data of all feedstock and AD digestate were then compared together. For microbial diversity analysis, the Shannon index was used. Box plots show medians ± interquartile ranges (IQR) and 1.5 IQR ranges (whiskers), with outliers denoted by open black circles. The overall differences in the bacterial community structures were evaluated by principal coordinate analysis (PCoA) based on Bray-Curtis dissimilarity values and performed by the R package PHYLOSEQ. The significant differences of gene relative abundance between feedstock and digestate samples were determined by the Wilcoxon signed-rank test at P < 0.05 [P-values were adjusted by Benjamini-Hochberg (BH) method]. The relationships between BGP operating temperature and the difference of ARG or pathogen relative abundance in digestate minus that in feedstock were analyzed by Spearman correlations using the R package PYSCH.

Results

A Comprehensive Gene Catalog of Feedstock Metagenome

Feedstock samples were collected before they were pumped into anaerobic digesters (Supplementary Figure 1), and a total of 56 feedstock samples were collected from the 56 full-scale biogas plants (Ma et al., 2021), with one sample from each plant. Deep sequencing of the 56 metagenomes generated a total of 1,716 Gb of high-quality data, with an average of 31 Gb per sample. Based on the assembled contigs with an average N50 length of 3,358 bp, we got a non-redundant Microbial Gene Catalog of Feedstock (MGC-F) containing a total of 25,296,776 genes, with an average open reading frame length of 770 bp and a full-length gene percentage of 53.44% (Table 1), which were comparable to those of Microbial Gene Catalog of AD (MGC-A) (Ma et al., 2021). A rarefaction analysis including all samples revealed a curve approached saturation with the increase of sample number (Figure 1A), suggesting that MGC-F covered the vast majority of microbial genes in these samples.

TABLE 1

Table 1. Comparison of microbial gene catalog of feedstock (MGC-F) and AD digestate (MGC-A).

FIGURE 1

Figure 1. Comparison of microbial gene catalogs of feedstock and AD digestate. (A) Rarefaction curve based on the gene profile of all samples. The two curves approach saturation at the end of sampling. (B) Taxonomic annotation of the gene catalogs at the superkingdom and phylum levels. (C,D) Venn diagrams of shared KOs and genes between the two gene catalogs. The gene catalog of AD digestate was cited from a previous study (Ma et al., 2021).

Taxonomic annotation of MGC-F showed that 79.63% of the genes were classified at the superkingdom level, a little higher than that for MGC-A (76.73%) (Ma et al., 2021). Of these, more than 97% were assigned to Bacteria, and the remaining genes were assigned to Archaea (1.65%) and Eukaryote (0.78%). Similar to MGC-A, the predominant bacterial phyla were Firmicutes (31.52%), Bacteroidetes (13.89%), and Proteobacteria (9.18%) (Figure 1B), though higher proportion of genes were assigned to Firmicutes in MGC-F (Ma et al., 2021). At lower taxonomic levels, there were only 14.87 and 1.05% of the genes were assigned to specific genera and species, respectively.

Functional classification of MGC-F showed that 56.93% of the genes were annotated with KEGG orthologous groups (13,595 KOs), which are comparable to those of MGC-A (Ma et al., 2021). By comparing the pairwise overlap of KOs between the two gene catalogs, we found that ∼ 93% of KOs were shared between MGC-F (93.2%) and MGC-A (93.7%) (Figure 1C), which were higher than that for genes ∼ 46% shared between MGC-F (48.6%) and MGC-A (44.0%) (Figure 1D). Among these shared KOs, about 43.7% of KOs were unknown, while about 6.5, 4.9, 4.8, 4.6, and 3.0% of known KOs were assigned to top five of the KEGG pathways of signal transduction, amino acid metabolism, membrane transport, carbohydrate metabolism, metabolism of cofactors, and vitamins, respectively. However, for 927 unique KOs in feedstock samples, they were assigned to about 880 unique orthologous groups.

Microbes and Functional Pathways for Methane Production Enriched From Feedstock to AD Digestate

In general, feedstock for anaerobic digestion, consisting of a complex mixture of microorganisms, organic matters, and inorganic materials, is an essential factor that drives microbial community variation in anaerobic digesters (Zhang et al., 2014). Principal coordinate analysis (PCoA) based on Bray-Curtis dissimilarity at genus level showed a clear separation between the samples derived from feedstock and AD digestate (Figure 2A), displaying the obviously different microbial communities. Alpha-diversity analysis (Shannon-Wiener index) revealed that the microbial diversity of feedstock samples was higher than that of AD digestate samples at the gene (adjusted P < 0.05), genus (Wilcox rank sum test P < 0.05), and KO (adjusted P = 0.12) levels (Figure 2B).

FIGURE 2

Figure 2. Comparison of microbial diversity between the microbial gene catalogs of feedstock and AD digestate. (A) Principal coordinate analysis (PCoA) at the genus level based on Bray-Curtis dissimilarity. The plots were separated into two clusters according to the feedstock and digestate. MCA, cattle manure biogas plants (BGPs); MCH, chicken manure BGPs; MPI, pig manure BGPs; OTH, BGPs with other feedstocks. (B) Microbial diversity (Shannon-Wiener index) at gene, genus, and KO levels. Asterisks denote significant difference between feedstock and digestate (adjusted P < 0.05). The relative gene abundance data for the gene catalog of AD digestate was downloaded from a previous study (Ma et al., 2021).

To compare the differences of microbial communities between the feedstock and digestate of AD, we used the relatively abundant top 50 genera in digestate samples, which accounted for about 90% of the relative abundance of annotated genera. In brief, the relative abundance ratio of each genus in digestate to that in feedstock (digestate/feedstock) for each biogas plant was calculated, and then the median ratio value of all these 56 biogas plants was counted. As a result, the relative abundance of only 14 genera in AD digestate samples were higher than that of feedstock samples (with median ratio value > 1; Figure 3A and Supplementary Figure 2), indicating that only a small number of genera were enriched in AD. These enriched genera include several methanogens like Methanoculleus, Methanospirillum, Methanosaeta, and Methanobacterium (Figure 3A), and their relative abundance in AD digestates were significantly (adjusted P < 0.05) higher than that of feedstock samples, suggesting their importance to the AD process. In addition, two syntrophic bacteria of Smithella and Syntrophomonas were also significantly (adjusted P < 0.05) enriched (Figure 3A). They were reported to syntrophically oxidize propionate and butyrate in association with methanogens, respectively (Xia et al., 2019).

FIGURE 3

Figure 3. The changes of microbial communities from feedstock to AD digestate. (A) The relative abundance ratio of each selected genus in digestate to that in feedstock (upper part) and the relative abundance of these genera (lower part). The median value was showing. Asterisks denote significant difference between feedstock and digestate (adjusted P < 0.05). (B) Relative abundance of genes involved in the three methanogenesis pathways of acetoclastic, hydrogenotrophic, and methylotrophic methanogenesis. Genes involved in only one pathway were marked by gray background, while genes shared by two (mtr) or all three pathways (mcr and hdr) were not marked. The relative gene abundance data for the gene catalog of AD digestate was downloaded from a previous study (Ma et al., 2021).

It is commonly accepted that methanogenic microorganisms are essential for AD process, and methanogenesis can occur via three main pathways, namely acetoclastic, hydrogenotrophic, and methylotrophic methanogenesis (Wintsche et al., 2018). The three methanogenesis pathways share some genes such as methyl-CoM reductase (mcr) and H4MPT-methyltransferase (mtr), but most of genes were unique in each methanogenesis pathway (Supplementary Figure 3). The relative abundance of most genes involved in acetoclastic methanogenesis pathway were higher than that in hydrogenotrophic and methylotrophic methanogenesis pathways in both feedstock and AD digestate (Figure 3B). In addition, the relative abundance of most genes involved in methanogenesis increased from feedstock to AD digestate, consistent with that many methanogens were increased from feedstock to AD digestate (Figure 3A).

To compare the microbial communities between the feedstock and AD digestate at genome level, metagenome-assembled genomes (MAGs) were constructed for feedstock samples. As a result, a total of 11,198 MAGs were generated from all 56 samples, with 3,015 MAGs have Completeness ≥80% and contamination ≤10%. After filtering the redundant genomes of the same species, a final of 2,101 representative MAGs were obtained, including 1,001 (47.6%) MAGs with Completeness ≥90% and contamination ≤5%. Taxonomic annotation of the MAGs revealed that 97.3 and 2.7% of MAGs were assigned to Bacteria and Archaea, respectively, which were consistent to that of AD digestate samples (Ma et al., 2021). The predominant phyla were Firmicutes (45.7%), Bacteroidetes (24.8%), and Proteobacteria (5.5%), which were consistent with the result derived from the gene catalog. In addition, MAGs generated from AD digestate samples (Ma et al., 2021) were compared to those from feedstock samples, and about 23% of MAGs in AD digestate were identified as the same species to those in feedstock samples (with ANI value > 95% and genome coverage >50% between the two MAGs). These results were consistent with fact that a huge number of AD populations were migrated from upstream feedstocks in BGPs (Mei et al., 2017).

Anaerobic Digestion Reduces the Dissemination of ARGs From Manure Wastes

The dissemination of antibiotic resistance genes (ARGs) from manure wastes to natural environment has become a global health issue (Youngquist et al., 2016; Wu et al., 2018). Existing studies have proposed AD as a promising method to reduce ARGs in manure (Wallace et al., 2018; Wu et al., 2018). However, only a few researches have studied the occurrences and changes of diverse ARGs in several full-scale BGPs treating manure (Luo et al., 2017), and not all ARGs can be effectively reduced in AD process. Thus, we used the metagenomic data to investigate the prevalence of all known ARGs in both AD digestates and the corresponding feedstocks of 56 BGPs. According to their substrate types, the 56 BGPs were divided into four groups: MCA (13 cattle manure BGPs), MCH (6 chicken manure BGPs), MPI (27 pig manure BGPs), and OTH (10 BGPs with other substrates) (Supplementary Table 1). Overall, AD reduced approximately 40% of the ARGs present in different feedstocks, with medians of the total ARGs relative abundance decreased by 51, 62, 36, and 48% in MCA, MCH, MPI, and OTH, respectively Figure 4A, possibly due to the decreased selective pressure on ARGs caused by the degradation of antibiotics in AD (Cheng et al., 2016; Turker et al., 2018; Huang et al., 2019). Besides, the changes of feedstock dilution, temperature, pH, ammonia concentration, may also influence the prevalence of ARGs. In particular, the decrease of most ARGs by AD was more evident in BGPs operated at relatively high temperatures typically 30–40°C (Supplementary Figure 4 and Supplementary Table 2).

FIGURE 4

Figure 4. Reduction of antibiotic resistance genes (ARGs) from feedstock to AD digestate. (A) Total relative abundance of ARGs in digestate and the corresponding feedstock of each BGP. MCA, cattle manure BGPs; MCH, chicken manure BGPs; MPI, pig manure BGPs; OTH, BGPs with other feedstocks. Hollow circles represent feedstock samples, and solid circles represent digestate samples. (B) Relative abundance ratio of 10 major ARG types in digestate to that in the corresponding feedstock of each BGP. Asterisks denote significant difference between feedstock and digestate (adjusted P < 0.05). The relative gene abundance data for the gene catalog of AD digestate was downloaded from a previous study (Ma et al., 2021).

AD significantly reduced (P < 0.05) the relative abundances of 15, 21, 23, and 20 ARG types in, respectively MCA, MCH, MPI, and OTH (Supplementary Table 3), and the reduction rates varied greatly among different ARG types. The highly abundant (>0.001%) ARGs dominated in all BGPs were the widely reported tetracycline, multi-drug, aminoglycoside, lincosamide, and macrolide resistance genes (Luo et al., 2017; Huang et al., 2019), with reduction rates of 20–70% in AD (Figure 4B and Supplementary Figure 5). Similarly, the reduction of some genes conferring tetracycline, aminoglycoside, and macrolide resistance was also observed in AD reactors treating cattle manure, domestic wastewater, and pig manure (Resende et al., 2014; Christgen et al., 2015; Ben et al., 2017). The less abundant (< 0.0001%) ARGs were cephalosporin, elfamycin, nitrofuran, nitroimidazole, and streptogramin resistance genes, with reduction rates higher than 95% in AD (Figure 4B and Supplementary Table 3). Moreover, we have detected the rarely reported phenicol and oxazolidinone resistance genes and they were also reduced by about 70% in AD. However, nucleoside and glycopeptide resistance genes were not reduced in AD (Figure 4B and Supplementary Table 3), similar to the persistence of sulfonamide and chloramphenicol ARGs in AD (Christgen et al., 2015). In summary, our results show that most ARGs, in particular the less abundant ARGs, are effectively reduced by AD, and the persistence of nucleoside and glycopeptide resistance genes deserves attention in future researches.

Elimination of Common Pathogens in Manure by Anaerobic Digestion

The outbreak of food-borne diseases threatens human health worldwide, partly due to the spreading of human pathogenic bacteria from animal manure (Altmann et al., 2011; Authority et al., 2017). Manure treatment by AD has been reported to reduce some culturable human pathogens (Massé et al., 2011; Fröschle et al., 2015; Orzi et al., 2015), while some other pathogens still persist in the digestate of AD (Owamah et al., 2014). Up to now, there lacks a comprehensive assessment on the risks of various human pathogens in AD, so we used the metagenomic data to investigate the change of all pathogens in the AD process of full-scale biogas plants. In total, we identified 79 human pathogenic bacteria (species) in the 56 BGPs, and digestates and corresponding feedstocks showed distinct difference in pathogen composition (Figure 5A). The medians of total pathogen relative abundance were reduced by 55, 85, 52, and 76% in MCA, MCH, MPI, and OTH, respectively (Figure 5B). Particularly, the BGP operating temperature showed negative correlations with the reduction of most pathogens in AD (Supplementary Figure 6 and Supplementary Table 4). On average, AD removed about 55% of the detected pathogens in different feedstocks, possibly because most pathogens are aerobic or facultative aerobic bacteria (Tang et al., 2015) whose growth is restricted by anaerobic environments and unfavorable high temperatures like 30–40°C. In addition, the abundance of pathogens may also be affected by the changes of feedstock dilution, pH, ammonia concentration, etc.

FIGURE 5

Figure 5. Reduction of pathogenic bacteria from feedstock to AD digestate. (A) Principal coordinate analysis (PCoA) based on Bray-Curtis dissimilarities of pathogens. MCA, cattle manure BGPs; MCH, chicken manure BGPs; MPI, pig manure BGPs; OTH, BGPs with other feedstocks. Hollow circles represent feedstock samples, and solid circles represent digestate samples. (B) Total relative abundance of 79 pathogens in the feedstock and digestate. Asterisks denote significant difference between feedstock and digestate (adjusted P < 0.05). The relative gene abundance data for the gene catalog of AD digestate was downloaded from a previous study (Ma et al., 2021).

For most of the highly abundant (0.001% ∼ 0.01%) pathogens in AD, their relative abundances were significantly (P < 0.05) reduced by 30–99% during AD process, and these pathogens included the widely reported Acinetobacter baumannii, A. lwoffii, Campylobacter coli, Clostridium perfringens, Enterococcus faecalis, E. faecium, Escherichia coli, Listeria monocytogenes, and Salmonella enterica (Eikmeyer et al., 2013; Fröschle et al., 2015), and the rarely reported Klebsiella pneumoniae, Pseudomonas aeruginosa, Staphylococcus aureus, and Streptococcus anginosus (Supplementary Table 5 and Supplementary Figure 7). Previously, the reduction of culturable species of Clostridium, Acinetobacter, Escherichia, Enterococcus, and Salmonella was also observed in AD reactors treating cattle manure (Resende et al., 2014; Qi et al., 2019). However, for other highly abundant pathogens of Bacteroides fragilis, Fusobacterium necrophorum, and Streptococcus pneumoniae, their relative abundances were not reduced during AD process (Supplementary Table 5 and Supplementary Figure 7). Therefore, the potential spreading of anaerobic B. fragilis, F. necrophorum and facultative anaerobic S. pneumoniae from manure or digestate deserves serious attention, for that they are the leading causes of intraabdominal infections, pharyngotonsillitis, pneumonia, and bacterial meningitis in children and the elderly (Sears et al., 2006; Janoff and Musher, 2015; Holm et al., 2016).

Discussion

Given the importance of anaerobic digestion (AD) in resource recycling and energy recovery, we constructed the first comprehensive microbial gene catalog of feedstock, and made a comprehensive comparison between the metagenome of feedstock and AD digestate. Our result showed that the predominant phyla in feedstock were Firmicutes, Bacteroidetes, and Proteobacteria, which were similar to that of AD digestate, and the microbial diversity of feedstock samples was higher than that of AD digestate samples. In addition, the relative abundance of most genes involved in methanogenesis increased from feedstock to AD digestate, consistent with that many methanogens were increased from feedstock to AD digestate. Besides, the obvious reduction of antibiotic resistance genes (ARGs) and pathogenic bacteria in AD was observed in almost all BGPs, suggesting the potential of AD to help control the environmental threat of ARGs and pathogens.

This study deepens the understanding of microbial compositions and functions in the feedstock and AD digestate, which has the potential to improve biogas production that can supply more clean energy, reduce the dependence on fossil fuels, and accelerate the development of the recycling economy. Furthermore, a good knowledge of the microbe-mediated process from organic matter to methane, a greenhouse gas (GHG) with a strong global warming effect, is crucial for the guidance of reducing methane emissions from various environments, such as livestock rumen, paddy field and peatland. Therefore, the results from this study will not only benefit the environment and socioeconomy but also provide some clues to resolving the challenge of global warming.

Data Availability Statement

The datasets presented in this study can be found in online repositories. The names of the repository/repositories and accession number(s) can be found below: https://www.ncbi.nlm.nih.gov/genbank/, PRJNA681092.

Author Contributions

FJ, YZ, SM, YH, HF, and QL collected the samples and performed the experiments. FJ, SW, YZ, and HW analyzed the data. AW and HL provide helpful suggestions. FJ, SW, and YZ wrote the raw manuscript. LC, YD, and WF conceived the study, designed the experiments, and revised the manuscript. All authors read and approved the final manuscript.

Funding

This project was supported by the grants from Agricultural Science and Technology Innovation Program (ASTIP), Chinese Academy of Agricultural Sciences (CAAS-ASTIP-2016-BIOMA), the Agricultural Science and Technology Innovation Program Cooperation and Innovation Mission (CAAS-XTCX2016), the Agricultural Science and Technology Innovation Program and The Elite Young Scientists Program of CAAS, the Fundamental Research Funds for Central Non-profit Scientific Institution, China (1610012016023), Fundamental Research Funds for Central Non-profit Scientific Institution (No. Y2017JC01), the Fund of Key Laboratory of Shenzhen (ZDSYS20141118170111640), the Infrastructure and Facility Development Program of Sichuan Province (2019JDPT0012), Science and Technology Program of Sichuan Province, China (2017JY0242), and Shenzhen Science and Technology Program (JCYJ20190814163805604). Central Public-Interest Scientific Institution Basal Research Fund (1610012017002_05103, Y2021PT02).

Conflict of Interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Acknowledgments

We thank Yanlai Liu, Xia Li, Bo Tu, Jing He, Yunfei Zhang, Lirong Dai, Lu Yang, Shouchao Lai, Nengmin Zhu, and Yinggang Zhang from Biogas Institute of Ministry of Agriculture for collecting samples. We also express our thanks to Jianjun Hu, Zuojun Liu, Liumeng Chen, Ling Qiu, and Xiaomei Ye for their assistance with sample collection, and Rui Chen, Yayue Deng, Yao Xiao, and Weidong Wu for assistance with physicochemical parameters measuring.

Supplementary Material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fmicb.2021.660225/full#supplementary-material

Abbreviations

AD, anaerobic digestion; ARG, antibiotic resistance gene; BGP, biogas plant; CARD, Comprehensive Antibiotic Resistance Database; GHG, greenhouse gas; Gb, gigabase; KO, KEGG orthologous group; MCA, cattle manure biogas plants; MCH, chicken manure biogas plants; MGC-A, Microbial Gene Catalog of anaerobic digestion; MGC-F, Microbial Gene Catalog of Feedstock; MPI, pig manure biogas plants; ORF, open reading frame; OTH, biogas plants with other feedstocks; PCoA, principal coordinate analysis.

Footnotes

1. ^ https://github.com/arpcard/rgi
2. ^ https://sourceforge.net/projects/bbmap/

References

Altmann, M., Wadl, M., Altmann, D., Benzler, J., Eckmanns, T., Krause, G., et al. (2011). Timeliness of surveillance during outbreak of shiga toxin-producing Escherichia coli Infection, Germany, 2011. Emerg. Infect. Dis. 17, 1906–1909.

Google Scholar

Angenent, L. T., Karim, K., Al-Dahhan, M. H., Wrenn, B. A., and Domiguez-Espinosa, R. (2004). Production of bioenergy and biochemicals from industrial and agricultural wastewater. Trends Biotechnol. 22, 477–485. doi: 10.1016/j.tibtech.2004.07.001

PubMed Abstract | CrossRef Full Text | Google Scholar

Authority, E. F. S., Authority, E. F. S., and Co, E. C. D. P. (2017). The European Union summary report on trends and sources of zoonoses, zoonotic agents and food-borne outbreaks in 2016. Efsa J. 15:5077.

Google Scholar

Ben, W. W., Wang, J., Pan, X., and Qiang, Z. M. (2017). Dissemination of antibiotic resistance genes and their potential removal by on-farm treatment processes in nine swine feedlots in Shandong Province. China Chem. 167, 262–268. doi: 10.1016/j.chemosphere.2016.10.013

PubMed Abstract | CrossRef Full Text | Google Scholar

Buchfink, B., Xie, C., and Huson, D. H. (2015). Fast and sensitive protein alignment using DIAMOND. Nat. Methods 12, 59–60. doi: 10.1038/nmeth.3176

PubMed Abstract | CrossRef Full Text | Google Scholar

Chaumeil, P. A., Mussig, A. J., Hugenholtz, P., and Parks, D. H. (2019). GTDB-Tk: a toolkit to classify genomes with the genome taxonomy database. Bioinformatics 36, 1925–1927.

Google Scholar

Cheng, W., Li, J., Wu, Y., Xu, L., Su, C., Qian, Y., et al. (2016). Behavior of antibiotics and antibiotic resistance genes in eco-agricultural system: a case study. J. Hazard. Mater. 304, 18–25. doi: 10.1016/j.jhazmat.2015.10.037

PubMed Abstract | CrossRef Full Text | Google Scholar

Christgen, B., Yang, Y., Ahammad, S. Z., Li, B., Rodriquez, D. C., Zhang, T., et al. (2015). Metagenomics shows that low-energy anaerobic-aerobic treatment reactors reduce antibiotic resistance gene levels from domestic wastewater. Environ. Sci. Technol. 49, 2577–2584. doi: 10.1021/es505521w

PubMed Abstract | CrossRef Full Text | Google Scholar

Eikmeyer, F. G., Rademacher, A., Hanreich, A., Hennig, M., Jaenicke, S., Maus, I., et al. (2013). Detailed analysis of metagenome datasets obtained from biogas-producing microbial communities residing in biogas reactors does not indicate the presence of putative pathogenic microorganisms. Biotechnol. Biofuels. 6, 49–63. doi: 10.1186/1754-6834-6-49

PubMed Abstract | CrossRef Full Text | Google Scholar

Fröschle, B., Messelhäusser, U., Höller, C., and Lebuhn, M. (2015). Fate of Clostridium botulinum and incidence of pathogenic clostridia in biogas processes. J. Appl. Microbiol. 119, 936–947. doi: 10.1111/jam.12909

PubMed Abstract | CrossRef Full Text | Google Scholar

Fu, L. M., Niu, B. F., Zhu, Z. W., Wu, S. T., and Li, W. Z. (2012). CD-HIT: accelerated for clustering the next-generation sequencing data. Bioinformatics 28, 3150–3152. doi: 10.1093/bioinformatics/bts565

PubMed Abstract | CrossRef Full Text | Google Scholar

Gerlach, W., and Stoye, J. (2011). Taxonomic classification of metagenomic shotgun sequences with CARMA3. Nucleic Acids Res. 39:e91. doi: 10.1093/nar/gkr225

PubMed Abstract | CrossRef Full Text | Google Scholar

Hassa, J., Maus, I., Off, S., Puhler, A., Scherer, P., Klocke, M., et al. (2018). Metagenome, metatranscriptome, and metaproteome approaches unraveled compositions and functional relationships of microbial communities residing in biogas plants. Appl. Microbiol. Biotechnol. 102, 5045–5063. doi: 10.1007/s00253-018-8976-7

PubMed Abstract | CrossRef Full Text | Google Scholar

Holm, K., Bank, S., Nielsen, H., Kristensen, L. H., Prag, J., and Jensen, A. (2016). The role of Fusobacterium necrophorum in pharyngotonsillitis – A review. Anaerobe 42, 89–97. doi: 10.1016/j.anaerobe.2016.09.006

PubMed Abstract | CrossRef Full Text | Google Scholar

Huang, P., Zhang, Y., Xiao, K. P., Jiang, F., Wang, H. C., Tang, D. Z., et al. (2018). The chicken gut metagenome and the modulatory effects of plant-derived benzylisoquinoline alkaloids. Microbiome 6:211.

Google Scholar

Huang, X., Zheng, J., Tian, S., Liu, C., Liu, L., Wei, L., et al. (2019). Higher temperatures do not always achieve better antibiotic resistance genes removal in anaerobic digestion of swine manure. Appl. Environ. Microbiol. 85, e02878–18.

Google Scholar

Hyatt, D., LoCascio, P. F., Hauser, L. J., and Uberbacher, E. C. (2012). Gene and translation initiation site prediction in metagenomic sequences. Bioinformatics 28, 2223–2230. doi: 10.1093/bioinformatics/bts429

PubMed Abstract | CrossRef Full Text | Google Scholar

Janoff, E. N., and Musher, D. M. (2015). “’201 – Streptococcus pneumoniae,” in Mandell, Douglas, and Bennett’s Principles and Practice of Infectious Diseases, 8th Edn, eds J. E. Bennett, R. Dolin, and M. J. Blaser (Philadelphia, PA: Content Repository Only!), 2310–2327.e5.

Google Scholar

Jia, B. F., Raphenya, A. R., Alcock, B., Waglechner, N., Guo, P. Y., Tsang, K. K., et al. (2017). CARD 2017: expansion and model-centric curation of the comprehensive antibiotic resistance database. Nucleic Acids Res. 45, D566–D573.

Google Scholar

Kanehisa, M., Goto, S., Kawashima, S., Okuno, Y., and Hattori, M. (2004). The KEGG resource for deciphering the genome. Nucleic Acids Res. 32, D277–D280.

Google Scholar

Kang, D. D., Li, F., Kirton, E., Thomas, A., Egan, R., An, H., et al. (2019). MetaBAT 2: an adaptive binning algorithm for robust and efficient genome reconstruction from metagenome assemblies. PeerJ 7:e7359. doi: 10.7717/peerj.7359

PubMed Abstract | CrossRef Full Text | Google Scholar

Li, D. H., Luo, R. B., Liu, C. M., Leung, C. M., Ting, H. F., Sadakane, K., et al. (2016). MEGAHIT v1.0: a fast and scalable metagenome assembler driven by advanced methodologies and community practices. Methods 102, 3–11. doi: 10.1016/j.ymeth.2016.02.020

PubMed Abstract | CrossRef Full Text | Google Scholar

Liu, B., Zheng, D., Jin, Q., Chen, L., and Yang, J. (2018). VFDB 2019: a comparative pathogenomic platform with an interactive web interface. Nucleic Acids Res. 47, D687–D692.

Google Scholar

Luo, G., Fotidis, I. A., and Angelidaki, I. (2016). Comparative analysis of taxonomic, functional, and metabolic patterns of microbiomes from 14 full-scale biogas reactors by metagenomic sequencing and radioisotopic analysis. Biotechnol. Biofuels. 9:51.

Google Scholar

Luo, G., Li, B., Lo, L.-G., Zhang, T., and Angelidaki, I. (2017). Antibiotic resistance genes and correlations with microbial community and metal resistance genes in full-scale biogas reactors as revealed by metagenomic analysis. Environ. Sci. Technol. 51, 4069–4080. doi: 10.1021/acs.est.6b05100

PubMed Abstract | CrossRef Full Text | Google Scholar

Ma, S., Jiang, F., Huang, Y., Zhang, Y., Wang, S., Fan, H., et al. (2021). A microbial gene catalog of anaerobic digestion from full-scale biogas plants. Gigascience 10:giaa164.

Google Scholar

Massé, D., Gilbert, Y., and Topp, E. (2011). Pathogen removal in farm-scale psychrophilic anaerobic digesters processing swine manure. Bioresour. Technol. 102, 641–646. doi: 10.1016/j.biortech.2010.08.020

PubMed Abstract | CrossRef Full Text | Google Scholar

Mei, R., Nobu, M. K., Narihiro, T., Kuroda, K., Munoz Sierra, J., Wu, Z., et al. (2017). Operation-driven heterogeneity and overlooked feed-associated populations in global anaerobic digester microbiome. Water Res. 124, 77–84. doi: 10.1016/j.watres.2017.07.050

PubMed Abstract | CrossRef Full Text | Google Scholar

Orzi, V., Scaglia, B., Lonati, S., Riva, C., Boccasile, G., Alborali, G. L., et al. (2015). The role of biological processes in reducing both odor impact and pathogen content during mesophilic anaerobic digestion. Sci. Total Environ. 526, 116–126. doi: 10.1016/j.scitotenv.2015.04.038

PubMed Abstract | CrossRef Full Text | Google Scholar

Owamah, H. I., Dahunsi, S. O., Oranusi, U. S., and Alfa, M. I. (2014). Fertilizer and sanitary quality of digestate biofertilizer from the co-digestion of food waste and human excreta. Waste Manag. 34, 747–752. doi: 10.1016/j.wasman.2014.01.017

PubMed Abstract | CrossRef Full Text | Google Scholar

Parks, D. H., Imelfort, M., Skennerton, C. T., Hugenholtz, P., and Tyson, G. W. (2015). CheckM: assessing the quality of microbial genomes recovered from isolates, single cells, and metagenomes. Genome Res. 25, 1043–1055. doi: 10.1101/gr.186072.114

PubMed Abstract | CrossRef Full Text | Google Scholar

Qi, G. D., Pan, Z. F., Yamamoto, Y., Andriamanohiarisoamanana, F. J., Yamashiro, T., Iwasaki, M., et al. (2019). The survival of pathogenic bacteria and plant growth promoting bacteria during mesophilic anaerobic digestion in full-scale biogas plants. Anim. Sci. J. 90, 297–303. doi: 10.1111/asj.13137

PubMed Abstract | CrossRef Full Text | Google Scholar

Qin, J., Li, R., Raes, J., Arumugam, M., Burgdorf, K. S., Manichanh, C., et al. (2010). A human gut microbial gene catalogue established by metagenomic sequencing. Nature 464, 59–65.

Google Scholar

Resende, J. A., Diniz, C. G., Silva, V. L., Otenio, M. H., Bonnafous, A., Arcuri, P. B., et al. (2014). Dynamics of antibiotic resistance genes and presence of putative pathogens during ambient temperature anaerobic digestion. J. Appl. Microbiol. 117, 1689–1699. doi: 10.1111/jam.12653

PubMed Abstract | CrossRef Full Text | Google Scholar

Sears, C. L., Franco, A. A., and Wu, S. (2006). “CHAPTER 28 – Bacteroides fragilis toxins,” in The Comprehensive Sourcebook of Bacterial Protein Toxins, 3rd Edn, eds J. E. Alouf and M. R. Popoff (London: Academic Press), 535–546. doi: 10.1016/b978-012088445-2/50033-0

CrossRef Full Text | Google Scholar

Stolze, Y., Bremges, A., Rumming, M., Henke, C., Maus, I., Puhler, A., et al. (2016). Identification and genome reconstruction of abundant distinct taxa in microbiomes from one thermophilic and three mesophilic production-scale biogas plants. Biotechnol. Biofuels 9:156.

Google Scholar

Tang, Y.-W., Sussman, M., Liu, D., Poxton, I., and Schwartzman, J. (2015). Molecular Medical Microbiology, 2nd Edn, Vol. 3. Amsterdam: Academic Press, 2097–2145.

Google Scholar

Turker, G., Akyol, C., Ince, O., Aydin, S., and Ince, B. (2018). Operating conditions influence microbial community structures, elimination of the antibiotic resistance genes and metabolites during anaerobic digestion of cow manure in the presence of oxytetracycline. Ecotox. Environ. Safe. 147, 349–356. doi: 10.1016/j.ecoenv.2017.08.044

PubMed Abstract | CrossRef Full Text | Google Scholar

Tyagi, V. K., and Lo, S. L. (2013). Sludge: a waste or renewable source for energy and resources recovery? Renew. Sust. Energ. Rev. 25, 708–728. doi: 10.1016/j.rser.2013.05.029

CrossRef Full Text | Google Scholar

Wallace, J. S., Garner, E., Pruden, A., and Aga, D. S. (2018). Occurrence and transformation of veterinary antibiotics and antibiotic resistance genes in dairy manure treated by advanced anaerobic digestion and conventional treatment methods. Environ. Pollut. 236, 764–772. doi: 10.1016/j.envpol.2018.02.024

PubMed Abstract | CrossRef Full Text | Google Scholar

Wintsche, B., Jehmlich, N., Popp, D., Harms, H., and Kleinsteuber, S. (2018). Metabolic adaptation of methanogens in anaerobic digesters upon trace element limitation. Front. Microbiol. 9:405. doi: 10.3389/fmicb.2018.00405

PubMed Abstract | CrossRef Full Text | Google Scholar

Wu, Y., Cui, E. P., Zuo, Y. R., Cheng, W. X., and Chen, H. (2018). Fate of antibiotic and metal resistance genes during two-phase anaerobic digestion of residue sludge revealed by metagenomic approach. Environ. Sci. Pollut. Res. 25, 13956–13963. doi: 10.1007/s11356-018-1598-x

PubMed Abstract | CrossRef Full Text | Google Scholar

Xia, X., Zhang, J., Song, T., and Lu, Y. (2019). Stimulation of Smithella-dominating propionate oxidation in a sediment enrichment by magnetite and carbon nanotubes. Environ. Microbiol. Rep. 11, 236–248. doi: 10.1111/1758-2229.12737

PubMed Abstract | CrossRef Full Text | Google Scholar

Youngquist, C. P., Mitchell, S. M., and Cogger, C. G. (2016). Fate of antibiotics and antibiotic resistance during digestion and composting: a review. J. Environ. Qual. 45, 537–545. doi: 10.2134/jeq2015.05.0256

PubMed Abstract | CrossRef Full Text | Google Scholar

Zhang, W., Werner, J. J., Agler, M. T., and Angenent, L. T. (2014). Substrate type drives variation in reactor microbiomes of anaerobic digesters. Bioresour. Technol. 151, 397–401. doi: 10.1016/j.biortech.2013.10.004

PubMed Abstract | CrossRef Full Text | Google Scholar