This meta-analysis includes sequences from published papers, all of which used bacterial 16S rRNA gene sequencing to survey the microbial community in full-scale anaerobic reactors. Seventeen studies were selected that comprised 62 methanogenic reactors: 27 reactors treating wastewater (WTAR) and 35 reactors treating solid waste (STAR). The configuration of the reactors, operational data, substrates, sequencing platforms, and primer sets are summarized in Table 1 (Supplementary Table 1). PCR-generated amplicons from bacterial 16S rRNA genes were sequenced using a diversity of primer sets targeting different regions of the 16S rRNA gene of bacteria and archaea (Table 1 and Supplementary Table 1). The data sets were retrieved from the Sequence Read Archive (SRA) of the National Center for Biotechnology Information (NCBI, United States) or were requested from the authors directly. The sample metadata for the processing of the raw data were extracted from published articles or provided by the authors. Some other operational data such as temperature, chemical oxygen demand removal (COD), pH, hydraulic retention time (HTR) and volatile fatty acid (VFA) were also collected in this study, but these data were only present for some reactors. Accordingly, we did not include those variables in this meta-analysis.

Raw reads (FASTQ) from each study were assessed for nucleotide quality and adaptors contamination with FastQC v0.11.5.¹ The raw reads were processed through Trimmomatic v0.38 (Bolger et al., 2014) for removal of ambiguous reads, adapters, and low-quality sequences to obtain high quality reads (slidingwindow:4:25, minlen:200). Depending on the sequencing platform (single-end or paired-end) the pre-processed raw data was imported in QIIME2 v.2019.10 (Bolyen et al., 2019). For demultiplexing the raw data, qiime cutadapt demux-paired or demux-single plugins were used with default settings. All chimeras and borderline chimeras were detected and discarded using the qiime vsearch uchime-denovo plugin (vsearch v2.7.0) with defaults settings (Rognes et al., 2016).

Given the heterogeneity of the selected hypervariable regions targeted within the 16S rRNA gene, we used closed-reference OTU picking method. In this method, the input sequences are clustered against a reference database as a means of being able to combine non-overlapping tags (Rideout et al., 2014). Because input sequences are not compared directly to one another (de novo OTU picking method), but rather to an external reference, the input sequences do not need to overlap. This is essential for this analysis because we are performing a meta-analysis including sequences derived from different amplification products of the same marker gene, such as the V1–V2 and V3–V4 regions of the 16S rRNA (e.g., as in the meta-analysis performed in Caporaso et al. (2010). As a result, for each short read we obtained a representative reference sequence of approximately 1500 bp, which enabled the comparison of non-overlapping tags. For this reason, it was not possible to use the ASV method in which overlapping sequences are needed. OTUs were then picked at 97% sequence identity against the SILVA 138 reference database using the plugin of vsearch cluster-features-closed-reference plugin (–p-perc-identity 0.97) (Quast et al., 2012; Glöckner et al., 2017). Taxonomy was assigned with qiime feature-classifier-sklearn (–p-confidence 0.8) plugin against the MiDAS 3.6 database (Nierychlo et al., 2020). MiDAS 3.6 database is a comprehensive ecosystem-specific reference database for activated sludge and anaerobic digesters which provides a taxonomic classification at all ranks for all sequences (Dueholm et al., 2019; Nierychlo et al., 2020). The singletons were deleted and each count table was then normalized to the lowest number of reads (561 reads per sample). Even though this sequence number is low, we decided to prioritize the inclusion of all samples and all variable regions amplified by the most frequently used primers. After the normalization, 10 reactors presented 0 Chloroflexi sequences which were the reactors with lowest relative abundance (below 0.5%) of Chloroflexi in the microbial community.

As equal variances are a prerequisite of permutational multivariate analysis of variance (PERMANOVA) (Daniels et al., 2011), previous to this analysis we performed a variance homogeneity check (betadisper function of the “vegan” package) (Oksanen et al., 2013). In order to decipher the significant influence of explicative factors on the microbial diversity (e.g., reactor type, platform, target region, and substrate, Table 1), a PERMANOVA with adonis function, and principal coordinates analysis (PCoA) based on Bray–Curtis dissimilarity distance was calculated using R software version 3.5.1 (R Core Team, 2013) in R Studio environment Version 1.0.153. Correction for multiple hypothesis testing was performed with the Benjamini–Hochberg (BH) procedure using the p.adjust function in R (Benjamini and Hochberg, 1995). Two-way analysis of variance (ANOVA) performed with R was used to test the effect of reactor type, platform, target region, and substrate on the relative abundance of phylum Chloroflexi. The differences and correlations were considered significant at p < 0.05. In order to determine significant difference in the relative abundance of Chloroflexi between STAR and WTARs, the Shapiro–Wilk (to test normal distribution) and Kruskal–Wallis (to test equality of means) tests were applied in R. The correlation between the relative abundance of Euryarchaeota and Chloroflexi was performed with the Spearman’s rho correlation test in R.

The diversity indices (Shannon H, Evenness, and Chao-1), and the t-test and Kruskal–Wallis test (for normal and non-normal distribution, respectively) to test the equal means over diversity indexes was performed using past software (v4.02) (Hammer et al., 2001). Boxplots and PCoA plots were generated using R package ggplot2 (Wickham, 2017).

A phylogenetic tree was constructed to evaluate the phylogenetic diversity of the phylum Chloroflexi. It included the reference OTUs obtained with the closed-reference OTU picking method with a relative abundance greater than 5% in at least one sample, reference OTUs which were part of the core microbiome (defined as Chloroflexi members composed by OTUs shared by 50% of the STARs or WTARs) and 16S rRNA sequences from assembled genomes and cultured representatives from the phylum Chloroflexi. Thermotoga spp. sequence served as an outgroup for rooting the tree. The phylogenetic tree was performed using the maximum likelihood phylogenetic tree method with the generalized time reversible (GTR) (Nei and Kumar, 2000) substitution model and GAMMA distribution model using 1000 bootstraps in MEGA 10.0 (Kumar et al., 2018). The resultant phylogenetic trees were visualized in iTOL (Letunic and Bork, 2019).

The coverage rates for Chloroflexi populations of the different primer sets used can differ with the target location within the 16S rRNA gene. The Silva database Test Prime tool (Klindworth et al., 2013) was used to evaluate primer bias over Chloroflexi populations at different taxonomic levels with settings recommended in the Silva database Test Prime tutorial (1 mismatch and 5 bases).

As different primer pairs were used to generate sequences in the different investigations, we analyzed how this might affect the recovery of Chloroflexi populations, at different taxonomic levels, frequently found in methanogenic reactors. The in silico analysis performed showed different coverage of the phylum when different primer pairs are used. The primer sets directed to the region between V3 and V5 hypervariable regions (47 reactors) showed the highest coverage for the Chloroflexi phylum (between 78 and 91%) (Table 2). Meanwhile, the coverage for families commonly found in methanogenic reactor as Anaerolineaceae, Ardenticatenaceae and Caldilineaceae was between 71 and 100%. The primer sets that covered hypervariable regions V1–V2 and V5–V9 (15 reactors) were the least efficient covering the Chloroflexi phylum (58–70%). For Anaerolineaceae, Ardenticatenaceae, and Caldilineaceae the coverage was between 61 and 100%, with the exception of V1–V2 region which showed a coverage of 5% for Anaerolineaceae (Table 2).

In total, 1,082,904 high quality sequences from 62 reactors, were retrieved from public databases, ranging from 561 to 164,994 (median of 11,299). All samples were randomly re-sampled at the same depth of 561 sequences before analyzing the data. To be able to compare sequences from different regions, closed-reference OTUs picking method was used. In this method, input sequences are aligned against a reference database. Because input sequences are not compared directly to one another, but rather to an external reference, they do not need to overlap. The closed-reference method allows performing a meta-analysis including sequences derived from different amplification products of the same marker gene. Clustering the pruned and combined data sets from all studies resulted in 13,165 no-singletons OTUs at a similarity threshold of 97%. Bacterial community diversity at phylum level was evaluated for all reactors. Proteobacteria, Chloroflexi, Bacteroidetes, and Firmicutes were the predominant phyla in WTAR and STAR (Figure 1). For WTARs the most abundant phylum was Proteobacteria, while Firmicutes dominated STARs. According to the PCoA analysis, the type of reactor (solid waste vs liquid wastewater treatment) influences the total bacterial community composition. PERMANOVA tests revealed that 15% of the variation of the total community was attributable to the type of reactor (PERMANOVA, p < 0.001) (Figure 2 and Table 3).

Regarding the abundance of the Chloroflexi phylum, we detected higher abundance in WTARs than in STARs. Moreover, Chloroflexi was detected in all reactors studied except for four STARs. In order to determine which variable significantly affected the relative abundance of Chloroflexi, two-way ANOVA analysis was applied. The variables tested were reactor type, primer pair target region, sequencing platform and substrate on the relative abundance of Chloroflexi. Our results revealed that the reactor type had the strongest impact on the abundance of Chloroflexi at phylum level (R² = 9.944, p < 0.001) followed by the target region (R² = 6.045, p < 0.001) (Supplementary Table 2).

Interestingly, the hypervariable region V3–V5 recovered the highest relative abundance of the phylum Chloroflexi (median value 7%) followed by V1–V2 (median value 6%) and V5–V9 (median value 0.5%) (Supplementary Figure 1). This is in accordance with the results obtained with the Silva Test Primer analysis (Supplementary Table 3).

The classification of Chloroflexi sequences at the class level showed that Anaerolineae class was predominant in all reactors independently of the substrate, reactor type, 16S rRNA-targeted region and sequencing platform (Figure 3A). The relative abundance of the Anaerolineae in the total community ranged between 1.4 and 31.1% (average: 13.2%) for WTARs and between 0.2 and 24.9% (average: 4.8%) for STARs. The difference in the median relative abundance of Anaerolineae was significantly different when comparing STARs and WTARs (Kruskal–Wallis test p-value < 0.001). Other minor classes detected were Chloroflexia and Dehalococcoidia present in some STARs and WTARs.

The diversity of the Chloroflexi community (Shannon index) was significantly higher in WTARs than in STARs (t-test, t = 2.32, p = 0.024) (Supplementary Table 4). The Chao1 and Evenness indices did not show significant differences (Kruskal–Wallis test, p > 0.05).

In this meta-analysis, most of the genera detected within the phylum Chloroflexi had no cultured representatives. The recent ecosystem-specific MiDAS 16S rRNA gene amplicon sequencing-based survey, facilitates the understanding of wastewater treatment ecosystem diversity and function (Nierychlo et al., 2019). The advantage of MiDAS 3.6 database is that it provides provisional alphanumeric names at genus and species levels allowing comparisons between unclassified members from different investigations. We found that the four most abundant genera were the same in STARs and WTARs with changes in their relative abundances: midas_g_467, midas_g_156, Bellilinea and Leptolinea (Figure 3B). While midas_g_467 was abundant in both systems, midas_g_156 was more abundant in WTARs. Despite that both belong to the T78 cluster and are therefore closely related.

Principal coordinates analysis and PERMANOVA tests revealed that 33% of the variation within Chloroflexi populations was attributable to the primer target region (PERMANOVA, p < 0.001) (Figure 4 and Table 3), while 10% of the variation was attributable to reactor type (Supplementary Figure 2 and Table 3). When performing PERMANOVA test, only with reactors in which the V3–V5 hypervariable region was amplified to avoid the primers bias, 12% of the variation of Chloroflexi populations was attributable to the reactor type (Supplementary Table 2). STARs and WTARs shared 23 of the 25 most abundant OTUs. STARs were dominated by midas_s_6158 (midas_g_467), midas_s_1625 (midas_g_467), midas_s_1462 (midas_g_467), midas_s_3887 (Leptolinea), and midas_s_667 (midas_g_667) (Figure 5 and Supplementary Table 5), while WTARs were dominated by midas_s_6158 (midas_g_467), midas_s_956 (midas_g_156), midas_s_6727 (Bellilinea), midas_s_3887 (Leptolinea), and midas_s_1462 (midas_g_467).

To determine if Chloroflexi populations were shared between reactors or if each reactor type had specific Chloroflexi populations we determined the minimum Chloroflexi core populations. Interestingly, we found that WTARs and STARs showed different Chloroflexi core populations as no OTU was shared by 50% of the reactors (Table 4). Within each reactor type, cores were found. For WTARs the core was composed by 18 OTUs belonging to Anaerolineaceae (midas_g_156, midas_g_467, midas_g_2702, Leptolinea, and Bellilinea genus) and midas_f_7837 (not assigned at genus level) families, representing at least five genera (Table 4). In STARs the core microbiome was composed by four OTUs belonging to the Anaerolineaceae family (midas_g_467). Only one OTU belonging to the Anaerolineaceae family (midas_g_467) was present in the core of both STARs and WTARs. It should be taken into account that these would be the most abundant OTUs. The low sequencing depth in some samples does not allow to detect OTUs with low abundances.

On the other hand, the correlation between Chloroflexi and Euryarchaeota (methanogens) was evaluated using the relative abundances obtained from Li et al. (2015), which included 18 reactors (515F-806R). Interestingly, we found a significant correlation between these two groups (rho: 0.60, p < 0.01) (Supplementary Table 6).

In order to determine the phylogenetic position of the most abundant Chloroflexi microorganisms, OTUs with a relative abundance greater than 5% in at least one reactor and OTUs belonging to the minimum core microbiome (defined as Chloroflexi members composed by OTUs shared between 50% of the STARs or WTARs) were used to construct a phylogenetic tree. Most of the selected OTUs sequences were positioned within the order Anaerolineales (52 OTUs) followed by midas_o_1 (5 OTUs). Half of the sequences positioned within the Anaerolineales order clustered with sequences from the T78 genus (22/52 OTUs), including most OTUs sequences from the core microbiome (Figure 6). The genus T78 was recently splitted into midas_g_156 and midas_g_467 (Nierychlo et al., 2020). Within the Anaerolineales order, the 16S ribosomal RNA gene retrieved from the assembled genome of Brevefilum fermentans CAMBI-1 (LT859958) was phylogenetically closely related to the T78 cluster. A few OTUs sequences within Anaerolineales were closely related to 16S rRNA gene sequences from isolated strains or assembled genomes members (Figure 6). The sequences affiliated to the order midas_o_1 had no closely related sequences.

Despite the bias generated from using different hypervariable regions we found that the Anaerolineae class largely dominated in both STARs and WTARs which was consistent with previous DNA based surveys (Rivière et al., 2009; Nelson et al., 2011; Bovio et al., 2019). The dominance of Anaerolineae in both types of systems might be explained through several hypotheses. First, most of the isolated strains of Anaerolineae present the potential to degrade cellulose, carbohydrates, and/or proteins anaerobically playing an important role as primary and secondary fermenters, having relevance in the bottlenecking hydrolysis step of anaerobic digestion (Sekiguchi et al., 2001, 2003; Yamada et al., 2005, 2006, 2007a; Narihiro et al., 2015). Thus, the fact that these systems in general treat different fermentable compounds might be favorable for the growth of Anaerolineae.

A second explanation to the predominance of Anaerolineae could be that their growth is significantly stimulated when co-cultivated with methanogens (Sekiguchi et al., 2001, 2003; Yamada et al., 2005, 2006, 2007a; Narihiro et al., 2015) indicating some potential for synergistic relationships. If the Anaerolineae group was a major donor of acetate or hydrogen to methanogens (Narihiro et al., 2015; McIlroy et al., 2017; Zamorano-López et al., 2020), the association in a “spaghetti-like” structure could promote the metabolites transfer flux between both groups, resulting in high methane production rates. In addition, the adhesiveness properties of Methanosaeta (Wiegant, 1986; Angenent et al., 2004; Zheng et al., 2006; Li et al., 2015) and members of Anaerolineae (Xia et al., 2016) could facilitate the metabolites transfer flux improving its dominance in these systems.

Interestingly, Anaerolineae had higher relative abundance in WTARs, which generally have granular biomass, than in STARs, in which the microbial biomass typically does not exhibit granulation but has instead dispersed biomass growth (Hofman-Bang et al., 2003; Klocke et al., 2007). In a previous study, we also determined that UASB reactors with granular biomass presented higher abundance of Anaerolineae than start-up UASB reactors with disperse biomass (Bovio et al., 2019). This might be explained by the fact that the filamentous morphology of Anaerolineae contribute to the granule formation in WTARs (Sekiguchi et al., 2001; Yamada et al., 2005; Satoh et al., 2007) presenting more abundance in reactors with granules than in those with dispersed biomass.

Another explanation could be the fact that the higher flow velocity used in WTARs compared to STARs, may select organisms which can adhere to each other to form well-settling granular sludge, favoring the growth of Anaerolineae in these systems over other microorganisms (Xia et al., 2016).

We found that the predominant genera and species in STARs and WTARs were the same, but changes in their relative abundances occurred. Four genera were predominant, two of them were closely related to previously described genera (Bellilinea and Leptolinea) and two of them did not harbor any cultured member and were named by the MIDAS database as midas_g_467 and midas_g_156. As for the isolated members of Anaerolineae, species belonging to the genus Leptolinea and Bellilinea present filamentous morphology and ferment carbohydrates and proteins (Yamada et al., 2006, 2007a). However, midas_g_467 and midas_g_156 previously named genus T78 in MiDAS database version 2 (Nierychlo et al., 2020), does not have representative cultured members. Despite this, some possible physiological features and ecological roles have been suggested for these genera in anaerobic reactors, which could explain the predominance of genus T78 (Kirkegaard et al., 2017; Zhu et al., 2017; Petriglieri et al., 2018; Jiang et al., 2020). Some of these previous reports suggest that their filamentous morphology and particular position within flocs or granules, might indicate that they play a role in maintaining floc structure (Petriglieri et al., 2018) as well as the granular structure (Zhu et al., 2017). It has also been reported that genus T78 could play an important role in hydrolysis and fermentation steps in the anaerobic digestion process (Yuan et al., 2015; Jiang et al., 2019). T78 also could metabolize alcohols and carbohydrates through syntrophic interactions (Praveckova et al., 2016). In the literature we found that two species belonging to the midas_g_467 correlated positively with different operational parameters: being total ammonia nitrogen, biogas yield, temperature, and organic loading rate (OLR) (Jiang et al., 2020). B. fermentans CAMBI-1 (LT859958) was the most closely related metagenome-assembled genome of the T78 cluster. The annotation of this genome suggested that members of the phylotype are important fermenters in mesophilic anaerobic digestion systems (McIlroy et al., 2017). Also, they are co-localized with the filamentous Archaea Methanosaeta spp. which suggests a potential undetermined synergistic relationship (McIlroy et al., 2017). Both, the potential fermenting metabolism as well as the potential syntrophism with archaea could explain the dominance of the T78 genus in anaerobic reactors. Our correlation analysis between Euryarchaeota and Chloroflexi supports this hypothesis. Hence, our main hypothesis about the predominance of the same four genera in both systems, relies on their capacity to undergo fermentative and hydrolytic metabolism. Further studies are needed to elucidate the function of these uncultured bacteria, which seem to have different responses to different conditions.

On the other hand, STARs and WTARs showed different core microbiomes of Chloroflexi members. In STARs the Chloroflexi core was less diverse comprising only midas_g_467 genus (4 OTUs) while in WTARs it was composed by at least midas_g_156, midas_g_467, midas_g_2702, Leptolinea, and Bellilinea genus (18 OTUs). This result was in accordance with the differences observed for the Shannon diversity index. The existence of a core microbiome based on the type of reactor (STARs/WTARs) despite of the heterogeneity of wastewater/solids and environmental factors (e.g., pH, temperature) may suggest that this core could be functionally important in each system. Specific core of microorganism between STARS and WTARs has been detected previously (Calusinska et al., 2018) but regarding Chloroflexi, no previous reports were found. A possible explanation for the different Chloroflexi cores found could be related to the biomass structure found in each type of system. Further exploration of this hypothesis is needed which would contribute to a better understanding of the role of Chloroflexi in the structure of granules and flocs and their metabolic role. Systematic studies examining multiple anaerobic reactor designs with greater depth of coverage and using platforms that generate long-reads or primer-free alternatives should help further define the “core microbiome” of Chloroflexi populations in methanogenic reactors. However, it is important to keep in mind that further recovery of isolated members or interpretative genomes is required to disclose the ecological importance of Anaerolineae as a core population in the anaerobic digestion process.

One of the bottlenecks in comparing microbial profiles is that different studies use different 16S rRNA gene regions generating an important bias. In the present study, we found an underestimation of Chloroflexi populations using primers targeting regions V1–V2 and V5–V9. An example was the underestimation of Anaerolineae class when V1–V2 regions were used recovering only a 5% of the total population. Meanwhile, primers targeting the V3–V5 region in general better reflected the abundances at all taxonomic levels. Previous studies investigating the influence of different hypervariable regions of 16S rRNA gene over the total community in activated sludge identified no perfect region concerning abundant functional genera, except for some regions with less bias, i.e., V1 and V2 (Cai et al., 2013; Guo and Zhang, 2013). However, (Albertsen et al., 2015) recommended primers targeting the V1–V3 regions which better reflected the abundances of Chloroflexi in activated sludge. More studies are needed to identify the best taxonomic profiling effectiveness between different hypervariable regions of the 16S rRNA gene over Chloroflexi populations in anaerobic reactors. As a general rule, drawing conclusions based only on one sequencing region should be avoided due to the potential false negative results, additionally amplicon sequencing of the 16S rRNA gene often limits resolution at genus level. We need primer-free alternatives to get the entire picture of the microbial diversity in anaerobic reactors (Karst et al., 2018). Using the full-length sequence has the potential to become a tool for more precise microbial community profiling that better allows global comparisons of microbiome studies and should potentially increase the accuracy and the resolution of closely related taxa. We highlight the effort to build the MiDAS reference database for microbes in wastewater treatment systems since it allows comparing uncultured microorganisms in different studies (Nierychlo et al., 2020).

Regarding the bias of different sequencing platforms, alpha diversity is significantly affected by both sequence length and depth (Singh et al., 2015). Platforms such as pyrosequencing which were discontinued, had less sequencing depth. To be able to compare samples with different sequencing depths, it is necessary to perform normalization at the lowest value of reads per sample. This leads to a loss of information in the less dominant groups which is an important limitation that has to be taken into account. This problem becomes relevant when comparing data from studies performed in the past (using outdated sequencing platforms) with data obtained more recently where the technology allows a much deeper sequencing. Moreover, it has been reported that Ion Torrent as well as Pyrosequencing platforms had comparatively higher error rates than Illumina platform (Siqueira et al., 2012; Salipante et al., 2014). Although it has been reported that while primer choice considerably influences quantitative abundance estimations, the sequencing platform has relatively minor effects when matched primers are used (Singh et al., 2015). Singh et al. (2015) reported that beta diversity metrics are surprisingly robust to both primer and sequencing platform biases.

New developments in single-cell genomics and metagenomics have in recent years provided new insights into the ecology and evolution of many novel uncultured microorganisms (Albertsen et al., 2013; Sekiguchi et al., 2015; Dam et al., 2020). The genomes have enabled the construction of metabolic models that attempt to explain the physiology of these organisms in detail. The genome-based models form the basis for more extensive investigations, such as in situ single-cell characterization, metatranscriptomics, and proteomics (Koch et al., 2014).

A second issue on performing a meta-analysis is the lack of standardized information on key operational parameters. It is important to have more reactor performance data associated with each reactor, such as compositions of influent waste/wastewater stream, COD loading/removal rate, pH, and temperature. These parameters should be fully used to infer potential roles of lineages of interest and should also be included. In our analysis it was not possible to correlate the microbial community with operational parameters data because in several articles this metadata was incomplete. For this reason, it would be essential to supply the operational parameters together with the raw data in each new study. In this sense, the creation of metadata repositories specific for anaerobic reactors including the 16S rRNA gene datasets must be a great advance. As far as we know there are some initiatives in other ecosystems as the Microbial Antarctic Resource Systems (MARS)² in Antarctic ecosystems.