Comparative Analysis of Bacterial and Archaeal Population Structure by Illumina Sequencing of 16S rRNA Genes in Three Municipal Anaerobic Sludge Digesters

2.1 Sampling

The waste sludge samples were collected from three anaerobic digesters of the Jebel Ali Sewage Treatment Plant (JASTP), UAE. The JASTP is the largest tertiary wastewater treatment infrastructure in Dubai, which serves a population of approximately 3.37 million and process 375,000 m³ of wastewater per day. All the three digesters had a capacity of 7433 m³ and were operating at a mesophilic temperature ranging between 32–37 °C [19]. The digesters were fed with 60 and 40% of raw and activated sludge, respectively. Details on the configuration and characteristics of the digesters are listed in Table 1. The samples from the anaerobic digesters were collected in autoclaved plastic bottles, and the bottles were transported to the lab on ice at 4 °C within an hour. The samples collected from three anaerobic digesters over a period of six months were designated as AD1, AD3, and AD5. The collected samples were immediately mixed with 100% ethanol in a ratio of 1:1 (v/v) and stored at − 20 °C until DNA extraction. Temperature, pH, and electrical conductivity (EC) were measured in sludge samples at the time of collection using HORIBA U-50 Multi Water Quality Checker (HORIBA Instruments Incorporated, USA).

2.2 Genomic DNA extraction

Genomic DNA was extracted in triplicate for each sludge digester sample using Power Soil DNA Extraction Kit (MO BIO Labs. Inc., Solana Beach, CA). The samples were vortexed and were then centrifuged at 4000 rpm for 5 minutes. DNA was then extracted from 0.25 g of pellet obtained according to the manufacturer's protocol. The extracted DNA was stored at -20 °C until further use. DNA concentration and purity were checked using the Qubit fluorometer (Thermo Fisher Scientific, USA).

2.3 Illumina sequencing of samples

The diversity of bacterial and archaeal communities in the samples was determined by amplifying the V3-V4 regions of bacterial and archaeal 16S ribosomal RNA (rRNA) genes. The following is a brief description of the main steps and protocol. The extracted genomic DNA samples quality check was performed by quantification using Qubit DNA BR Assay kit (Thermo Fisher Scientific, USA Cat#Q32853). For the generation of 16S amplicon, the extracted DNA samples were diluted to 10 ng and were amplified for 16S (~ 1500 bp) using 16S (5' AGAGTTTGATCCTGGCTCAG 3') and 16S reverse primers (5' GGTTACCTTGTTACGACTT 3') and positive control (internal metagenomic DNA sample) and no template control. These amplicons were checked on 1% agarose gel. To generate V3-V4 amplicon, 16S amplicon was used as a template with all the samples subjected to V3-V4 amplification (~ 460 bp) using V3-V4 forward and V3-V4 reverse primers (primer sequences V3-V4-Forward 5' CCTACGGGNGGCWGCAG 3' and V3-V4 Reverse 5' GACTACHVGGGTATCTAATCC 3') and a positive control (internal metagenomic DNA sample) and without template control [20]. The amplicons were checked on 1% agarose gel. The V3-V4 amplicons were then cleaned using AMPure XP beads (Beckman Coulter, CA, USA, Cat# A63882) to get rid of non-specific fragments. The V3-V4 products were used for DNA library preparation using NEBNext Ultra DNA Library Prep Kit for Illumina (New England Biolab, UK, Cat# E7370L). First, the amplicons were end-repaired and mono-adenylated at 3' end in a single enzymatic reaction. Next, NEB hairpin-loop adapters are ligated to the DNA fragments in a T4-DNA ligase-based reaction. Following ligation, the loop containing Uracil is linearized using USER Enzyme (a combination of UDG and Endo VIII), to make it available as a substrate for PCR based indexing in the next step. During PCR, barcodes were incorporated using unique primers for each of the samples, thereby enabling multiplexing. The prepared libraries were checked for fragment distribution using D1000 Screen Tapes (Cat# 5067–5582, Agilent, CA, USA) and reagents (Cat# 5067–5583, Agilent, CA, USA). The obtained libraries were pooled and diluted to final optimal loading concentration before cluster amplification on Illumina flow cell. Once the cluster generation is completed, the clustered flow cell is loaded on Illumina HiSeq2500 instrument (Illumina, Inc., San Diego, USA) for amplicon sequencing to generate 0.5M, 250 bp paired-end reads per sample using pair-end approach.

2.4 Bioinformatics analysis

The schematic bioinformatics analysis pipeline is represented in supplementary Fig. 1. Briefly, the following steps were involved in the bioinformatics analysis for further processing of the reads and quality filtering sequencing data. The quality checking of the raw fastq files was carried out using FASTQC to check for base quality, base composition and G.C. content. Based on the quality report of fastq data, the sequence reads were trimmed using fastq-mcf where necessary only to retain high-quality sequences for further analysis. Besides, the low-quality sequence reads were excluded from the analysis.

The 16S rRNA gene identification was conducted by identification of the amplicons from the raw fastq files based on the presence of conserved region forward and reverse. The paired-end sequences read contained a spacer region, some portion of the conserved region, and V3-V4 region. As a first step, the spacer and conserved regions were removed from paired-end reads. After trimming the unwanted sequences from original paired-end data, a consensus V4 region sequence was constructed using Fast Length Adjustment of SHort reads (FLASH V1.2.7) described by[21].

The pre-processing of reads was done by dereplication step. The dereplication is the identification of unique sequences so that only one copy of each sequence is reported. USEARCH is used for dereplication with the derep_full length. Chimeric sequences are formed as a result of experimental artifact, which is a fusion between two or more original sequences. UCHIME utility from usearch was used to remove the chimeras using the de novo approach[22].

The operational taxonomic unit (OTU) identification and singleton removal was done by performing the following steps. The sequences having a similarity of 97% were grouped using a closed reference method under a single operational taxonomic unit (OTU) for classification of the whole data. Any OTU that has a count of 1, i.e., only one sequence present once in a single sample (singleton) were filtered out as these may have been created as an experimental artifact rather than being some novel organism. For taxonomic classification, one representative sequence from each OTU was picked and classified using the RDP classifier against the green gene database at 97% similarity. The taxonomy classification was done at phyla, order, family, and genera level. The raw sequences have been deposited in National Center for Biotechnology Information (NCBI Bethesda MD, 20894 USA) under Bio-Project accession number PRJNA602372.

2.5 Statistical analysis

In the present study, two types of comparisons are carried out. Firstly, intergroup comparisons were carried out to find differences in genera across each of the groups as follows: AD1 v/s AD3; AD3 v/s AD5 and AD1 v/s AD5. The comparisons are performed using a t-test, and any genus having a p-value less than 0.05 in the t-test is considered as differentially present. Secondly, intragroup comparisons are carried out for each of the three groups, between samples collected at a different period. As these comparisons are between a pair of samples (single sample comparison) without any replicate, a different approach is taken to identify the differential genera. First, the unique genera for both the samples are listed out separately. Then, a ratio of relative abundance for both the samples are calculated for common genera between the samples. These ratios are then converted into log scale to get a normal distribution, and any genera having a value beyond mean ± 2SD are considered as differentially present genera.

3.1 Physicochemical conditions of Anaerobic digesters

The physicochemical conditions of the digester operation are given in Table 2. All three digesters were operated in almost similar physicochemical conditions, i.e., under mesophilic temperature (34 °C), neutral pH range of 7.13–7.55. At the same time, dry solids, volatile solids concentrations varied respectively, between 2.56 to 5.74% and 43.75 to 70.95 mg L^− 1. Among the three digesters, however, AD3 had a relatively low organic loading rate of 5.84 kg. ODS/m³. d compared to AD1 (6.84 kg. ODS/m³. d) and AD5 (6.61 kg. ODS/m³. d). The electrical conductivity (EC) of all sludge digestors found to be towards the higher range. The highest conductivity of sludge samples was found for AD5 (12.33–15.46 mS cm^− 1) followed by AD3 (11.71–14.44 mS cm^− 1), and AD1(9.77–13.45 mS cm^− 1). The EC value range found in this study was similar to a study conducted recently in Austria in which the effects of various co-substates were evaluated for the microbial community composition of seven full-scale anaerobic digesters fed with or without co-substrates [17]. Overall, three anaerobic digesters showed acceptable levels of main operational parameters and were performing stably during the sampling period.

3.2 Illumina sequencing of 16S rRNA genes

To check the quality of Illumina sequencing results, the number of raw reads, total data obtained, and the average base quality for each sample are shown in supplementary Fig. 2. It can be observed that all the samples, except two of them, have an average Phred score of more than 30. To retain only the high-quality data reads with a Phred score > Q30 is taken for down-stream analysis. The total number of raw reads of 16S rRNA genes identification, and CR (conserved region for V3-V4) trimmed reads are presented in supplementary Fig. 3. Overall, 3340 MB sequencing data were generated using Illumina sequencing, revealing 1.6 × 10⁶ raw reads. Using quality filtering, 60% of the total reads were removed, resulting in an average of 0.6 × 10⁶ reads per sample. A total of 419 bacterial community taxa found and classified at Phyla (33), Class (64), Order (99), Family (116), and Genera level (107). To find out the major taxa in each level, filtration is applied to retain taxa, which are present at least 1% in any of the samples, resulting in 84 OTUs based on 97% similarity. While 6.1 ± 1.24% reads were assigned to Archaea, which included five major phyla, 83 Genera, and 152 species. Bacteroidetes, Firmicutes, Synergistetes, and Theromotogate, OP8, and Chloroflexi were found to be dominant in all the digesters, accounting for 90.85% of all sequences (Fig. 1). In this study, Bacteroidetes were found to be the most abundant phylum (41.44%), followed by Firmicutes (23.42%), Synergistetes (7.22%), Thermotogae (6.70%), and Chloroflexi(6.69%) and OP8 (5.3%). Other phyla, detected at percentages > 1% were affiliated to Actinobacteria (3.1%), Proteobacteria (1.4%), and WWE1 (1.21%). The overall abundance of bacterial community phyla observed in this study is similar to closely related studies where a group of Phyla, namely Chloroflexi, Firmicutes, and Bacteroidetes were found to be dominant in full-scale anaerobic sludge digester samples [17, 23, 24]. However, unlike these studies, the percentages of Bacteroidetes (44%) detected in this study were significantly high and the population of proteobacteria was very low (1.4%). The higher population of Bacteroidetes phyla can be associated with high hydrolytic activity in full-scale anaerobic digesters[25].

3.3 Bacterial community diversity at genus level

The rarefaction curves show that among the three digesters, the highest diversity was present in AD5. At the same time, the least diverse microbial community was found in AD1 (Fig. 2A). Fifty-one genera were detected as core genera in all the digesters, including Acetobacteroides (Blvii28), Coprothermobacter, Fervidobacterium, Clostridium, Caldilinea, Allochromatium, Sediminibacter, and T78 (Fig. 2B). In addition to these genera, many of the bacteria were identified as unknown bacteria. Most of these genera are associated with an anaerobic digester, some have been reported from anaerobic digester earlier also[11]. Notably the highest population of genus Acetobacteroides (Blvii28) was present in all the digesters. The cultured representative member of the genus is known to produce acetate, H_2, and CO₂ as the product of fermentation [26]. Intra-group comparison of three AD3 samples collected between September 2016 and January 2017 showed that the concentration of Blvii28 increased significantly in samples from September 2016 to December 2016 and January 2017 (supplementary file 2). A similar trend was observed for intragroup comparison of AD5 samples collected between September 2016 and January 2017 (supplementary file 3). In contrast, T78 is possibly metabolizing carbohydrates and alcohol via syntrophic interactions [27]. T78 is the member of the phylum Chloroflexi often found to be abundant in full-scale mesophilic anaerobic digesters receiving excess sludge from wastewater treatment plants[28]. Intragroup comparison of AD5 (Jan-Feb 2017) samples found nine-fold increase in concentration of T78 in February 2017 sample compared to January 2017. Many genera found predominantly in the digesters were thermophilic or halophilic. Like, Coprothermobacter is a known proteolytic anaerobic thermophilic bacteria found in many thermophilic anaerobic digester [29]. This genus can also improve protein degradation by establishing a syntrophy with hydrogenotrophic archaea [29]. Other thermophilic bacteria found in the anaerobic digester include Fervidobacterium and Caldilinea, that ferment carbohydrates to lactate, acetate, hydrogen, and carbon dioxide [30, 31].

In comparison, genus Sediminibacter which was initially isolated from marine sediment, was also found as one of the dominant genera [32]. This can be due to the simple reason that most of the water used in Dubai is obtained from the sea and is being used after desalination [33, 34]. The incoming feed may contain elevated levels of salts, or it may be due to the fact that Sediminibacter has a unique light-driven sodium ion pump, which helps in their survival in marine habitats [35]. Several purple sulfur bacteria like Allochormatium and Thermotogales AUTHM297 were also part of the microbial community in these digesters. Many sulfate-reducing bacteria were also found, including Desulfomicrobium and Desulfobacter, as the two most dominant genera of SRB. In our previous work using FISH technique, we found that in all three digesters members of the genus Desulfobacter and Desulfobacterium were consistently present in large numbers[19]. It is interesting to note that Desulfomicrobium is also known to be associated with the marine habitat [36].

The process of anaerobic digestion involves the degradation of organic compounds to simple organic compounds, which is most likely carried out by the members of the phylum Bacteroidetes, Firmicutes, Proteobacteria, Synergistetes, and Thermotogae. Other dominant phyla in the fermenters are shown in Fig. 3. Previous studies also demonstrate the presence of these phyla in anaerobic digester as macromolecules degrading bacteria[11]. Mycobacterium, which is associated with the production of Lipases and Lipolytic activity, is also found as one of the predominant genera. Another critical process of anaerobic digestion is acidogenesis, wherein bacteria convert organic monomers into acids like acetic, propionic, and butyric acids. The population of Acidogenic bacteria like Acetobacteroides, Fervidobacterium, Clostridium, and Paludibacter producing acetate lactate or propionate was high in the digesters as reported above. The conversion of these organic molecules to CH₄ is carried out mainly by Archaea. However, some bacteria may influence the production of methane. For example, acetate-utilizing uncultured bacteria of Synergistes group 4 and its competition with aceticlastic methanogen [37].

The boxplots of intergroup comparisons for differentially present genera among three anaerobic digestor groups (AD1 v/s AD3; AD3 v/s AD5 and AD1 v/s AD5) are shown in Fig. 4 and supplementary file. Three genera Rhodobacter, Allochromatium, and Bifidobacterium, are found to be differentially present across the two groups AD1 and AD3, but among them, only Bifidobacterium showed differential abundance between AD3 and AD5. For the third Comparison (AD1 v/s AD5), genera Propionicimonas, Clostridium, and Paracoccus showed differential abundance between the two groups AD1 and AD5. To find out if the samples showed an overall group-specific trend, a principal component analysis based on all the genera detected is performed. The PC1 vs PC2 PCA plot (Fig. 5) showed no group-specific signature. Overall, from the taxonomic classification at the genera level, it can be observed that different digesters alter the microbiome profile differently. For intergroup comparisons, no abundant genera (present at least 1% in any of the samples) showed significant differentiation.

3.4 Archaeal community diversity at genus level

When analyzed at the genus level, the highest diversity was observed in digester AD5, followed by AD1 and AD3 (Fig. 6A). This was evident from the highest Shannon index obtained for AD1, followed by AD3 and AD5. The predominant Archaea in the three digesters were the members of the phylum Euryarchaeota followed by Crenarchaeota. The population of hydrogenotrophic archaea was higher in the digesters. Methanocorpusculum alone constitutes more than 50% of the total population of Archaea. Other predominant genera of archaea include Metallosphaera, Methanocella, Methanococus, Acidianus, Natronobacterium, and others shown in Fig. 6B. Methanocorpusculum is one of the hydrogenotrophic methanogens, which was isolated for the first time from the biodigester of the wastewater treatment plant [38]. In contrast, the second most predominant genus was Metallosphaera, which is an extreme thermoacidophile with optimal growth at 74 °C and pH 2.0 [39]. Methanocella, which was the third most dominant genus, is a mesophilic, hydrogenotrophic methanogen [40]. While Methanococcus is also a thermophilic hydrogenotrophic methanogen and at least one species of the genus is also known to fix nitrogen[41].

It is clear from the microbial community analysis that a core microbial community exists in the anaerobic digesters studied. However, the degradation is carried out by the members of various phylum reported earlier also from the anaerobic digesters. These include the members of the phylum Bacteroidetes, Firmicutes, Proteobacteria, Synergistetes, and Thermotogae [13]. While several genera mainly involved in the conversion of these organic substrates into acetate, butyric acid, propionate, and H₂ were found to be the dominant genera in the anaerobic digesters. However, Euryarchaeota, mainly the hydrogenotrophic methanogens, were found as predominant members in the digesters (Fig. 6). The predominant hydrogenotrophic methanogens activities compared to aceticlastic ones, is in agreement with another study[25]. Overall, abundant microbial genera observed in this study, and their predicted role in anaerobic digestion associated with a specific process is shown in Fig. 7. Although the digesters were operated under mesophilic conditions, many thermophilic genera of bacteria and archaea were found to predominate in the digester. It may be because of the presence of these bacteria in the incoming feed. Whether or not these bacteria are initially present in the feed is a matter of further investigation. Furthermore, the bacteria of marine origin like Sediminibacter and Haloterrigena as one of the predominant genera maybe because most of the water used in Dubai is desalinated seawater. The high salinity range suggested indirectly by higher electrical conductivity values observed in this study may be the reason for the occurrence of halophilic bacteria found in anaerobic digester systems.