Evolution of bacterial consortia in an integrated tannery wastewater treatment process

Abstract

In the present study, the dynamics of microbial communities and their abundance associated with each stage of a tannery wastewater treatment process were investigated by polymerase chain reaction-denaturing gradient gel electrophoresis (PCR-DGGE) combined with high throughput sequencing. Both PCR-DGGE and high throughput sequencing results reflected the bacterial succession in the integrated treatment process and phyla of Proteobacteria, Bacteroidetes and Firmicutes dominated throughout the integrated treatment process. However, Actinobacteria, Planctomycetes and Chlofoflexi only predominant the anoxic/oxic (A/O) process. The bacterial richness during the A/O process was higher than that in the pre-treatment process. Moreover, quantitative polymerase chain reaction (qPCR) analysis indicated that the absolute abundance of 16S rRNA genes in the biological treatment stages were higher than in other stages. Finally, redundancy analysis suggested that [Thermi] should be involved in NH₄⁺–N removal and ammonia concentration had positive effects on the bacterial diversity. Overall, this study provided insight into the evolution of the bacterial community structure and diversity in integrated wastewater treatment processes and identified the correlations between the physicochemical characteristics of wastewater and bacterial community structures.

---

1. Introduction

As one of the most polluting sources in industry, the tanning industry generally produces plenty of wastewater which has serious effects on the environment. Tannery wastewater is complex due to a multiplicity of different chemicals used during leather production and is mainly characterized by the presence of high salinity, ammonia, organic loading and specific pollutants such as chromium and synthetic tanning agents.^1,2 Ammonia present in effluent can result in eutrophication in rivers, lakes and reservoirs, leading to the excessive growth of algae and aquatic plants and depletion of dissolved oxygen.³ Consequently, the uncontrolled released of tannery wastewater into the environment could have a harmful impact on the environment and there is an urgent requirement to efficiently treat tannery wastewater before its discharge.

In general, biological wastewater treatment which utilized the metabolism of microbiota to remove the pollutants is more widely used in wastewater treatment plants (WWTPs) than physical and chemical methods due to its low cost on energy and operation as well as high efficiency.^4,5 However, conventional biological treatment cannot completely remove the pollutants, especially tannins due to its low biodegradability, leading to the substandard effluent.⁶ The combination of a periodic submerged filter and ozone oxidation were applied for 97% and 98% removal on chemical oxygen demand (COD), NH₄⁺–N, respectively.⁷ Additionally, anoxic/oxic (A/O) process was considered as a desirable option with greater cost effectiveness and higher efficiency than conventional wastewater treatment.⁸ High removal rates of COD and nitrogen ammonia can be obtained by A/O system.^9,10 Noteworthy, the composition and activity of microbial community, which is dominate by prokaryotic microorganisms play essential roles in the efficiency and robustness of biological treatment reactors.¹¹ Therefore, better understanding of the structure and dynamics of microbial community in WWTPs will not only be conducive to clarify the mechanisms of nutrient and nitrogen removal, but also contribute to improve the performance and stability of biological wastewater treatment bioreactors.¹²

Currently, the microbial consortia in activated sludge from different industrial and domestic wastewater have been investigate by nested polymerase chain reaction-denaturing gradient gel electrophoresis (PCR-DGGE) approach and pyrosequencing approach.^13,14 Nevertheless, still comparatively little is known of microbial diversity in tannery wastewater treatment plants. Recently, the comprehensive structures of bacterial nitrifiers and denitrifiers in anaerobic and aerobic sludge from two tannery sewage treatment plants were investigated and the results showed that functional genes were more abundant in aerobic sludge.¹⁵ However, they only focused on the differences on microbial structures among the certain tanks in different WWTPs. Research on the microbial community structure in an integrated treatment system by 16S rRNA clone library showed that Firmicutes was the most represented phylum in anaerobic–aerobic reactors.¹⁶ Nevertheless, still comparatively little is known of how each step impacts on the bacterial composition. Therefore, it is worthwhile to investigate the evolution of microbial community in integrated process and their linkages to the characteristics of wastewater.

In the present study, 16S rRNA genes of microorganisms in different sections of the integrated process were systematically analyzed by nested PCR-DGGE fingerprinting and Illumina HiSeq sequencing to investigate the changes of bacterial diversity and abundance in the system of a full scale WWTP treating tannery sewage in Shandong Province. Principal components analysis (PCA)/principal coordinate analysis (PCoA) and cluster analysis were used to appraise α- and β-diversity among the samples. Subsequently, the absolute quantity of 16S rRNA gene copy numbers of total bacteria for each sample was quantified by quantitative polymerase chain reaction (qPCR). Then, redundancy analysis (RDA) and Spearman's rank correlation coefficient (SRCC) were conducted to investigate the correlations between microbial community composition and environmental factors.

2. Materials and methods

2.1 Mixed-liquid samples collection

Mixed-liquid samples were collected from all stages of an integrated tannery wastewater treatment process located in Binzhou City, Shandong Province, China (37°22′N 118°01′E). Details of the treatment process and sampling sites were shown in Fig. 1 and Table 1. Among these tanks, grit chamber (GC) was used to wipe off the large particles and grits in sewage. Subsequently, the settleability of suspended substances was improved with the addition of polymeric aluminum sulfate (PAS) and polyacrylamide (PAM) in homogenize tank (HT) and suspended substances deposited in primary sedimentation tank (PST) for primary removal of contaminants. Secondary sedimentation tank (SST), final sedimentation tank (FST) and central water tank (CWT) were used to separate the sludge and liquid. As the functional tanks, anoxic tank (AT) and oxic tank (OT) were mainly responsible for the removal of nitrogen and organic matter.

Fig. 1 The schematic diagram of sampling sites along the tannery wastewater treatment process. Sampling sites: GC: grit chamber. HT: homogenize tank. PST: primary sedimentation tank. AT: anoxic tank. SST: secondary sedimentation tank. OT: oxic tank. FST: final sedimentation. CWT: central water tank.

Table 1 Characteristics of wastewater in each tank of the integrated process^a

	GC	HT	PST	AT	SST	OT	FST	CWT
a NA: not available, RE: removal efficiency.
COD (mg L^−1)	2476.00	2899.33	2091.67	573.00	384.33	267.67	202.00	93.47
NH4^+–N (mg L^−1)	335.00	295.00	345.00	345.33	358.67	36.67	3.33	1.21
pH	7.44	7.73	7.05	7.29	6.28	7.37	7.83	7.26
COD RE (%)	NA	−17.10%	27.86%	72.61%	32.93%	30.36%	24.53%	53.73%
NH4^+–N RE (%)	NA	11.94%	−16.95%	−0.10%	−3.86%	89.78%	90.91%	63.70%

---

We trisected the diagonal and depth of each tank, then collected 1 liter mixed-liquid sample at each point. A total of 9 liters mixed-liquid were mixed as one sample for each tank. The mixed liquid were then taken in sterile polypropylene bottles and immediately transported to the laboratory on ice for further treatment. For further analysis of microbial diversity, each sample was dispensed into a sterile Eppendorf tube and centrifuged at 10 000g for 10 min, and then decanted the supernatant. The pellet was stored at −80 °C prior to further analysis.

2.2 Chemical analysis

The physicochemical characteristics of the mixed liquids were analyzed by the standard methods.¹⁷ The pH was measured by Sartorius PB-10 pH meter (Sartorius, Germany); COD was measured by the potassium dichromate method: successively add 3 mL of sample, 1 mL of masking agent (300 g L⁻¹ HgSO₄ in 10% m/v H₂SO₄), 3 mL of digestion solution (9.8 g L⁻¹ K₂Cr₂O₇, 50 g L⁻¹ KAl(SO₄)₂·12H₂O and 10 g L⁻¹ H₈MoN₂O₄) and 5 mL of catalyst solution (4.4 g L⁻¹ Ag₂SO₄ in 98% m/m H₂SO₄) in a digestion tube and immediately incubated with XJ-III digester device (Tomorrow Environmental Protection Instrument Co., LTD, China) at 160 °C for 25 min. After the digestion solution cooled to ambient temperature, the solution was determined by titration with standard ammonium iron(II) sulfate solution (9.88 g L⁻¹ (NH₄)₂Fe(SO₄)₂·6H₂O in 0.2% m/v H₂SO₄) exhibiting bronzing titration endpoints. The COD value of sample by equation:
where V₀ and V₁ were the volume of standard ammonium iron(II) sulfate solution for control and sample, respectively. V₂ was the volume of sample. C was the concentration of standard ammonium iron(II) sulfate solution. NH₄⁺–N was measured by using Nessler's reagent spectrophotometry: 1 mL zinc sulfate solution (100 g L⁻¹) was added into 100 mL sample and then adjusted its pH to 10.5 with NaOH (250 g L⁻¹). After the precipitation, 50 mL supernatants were mixed with 1 mL zinc sulfate solution (500 g L⁻¹) and 1 mL Nessler's reagent (160 g L⁻¹ NaOH, 70 g L⁻¹ KI and 100 g L⁻¹ HgI₂), the mixtures were inoculated at ambient temperature for 10 min and the absorbance was measured at 420 nm. Ammonium chloride standard solution with a theoretical NH₄⁺–N value of 0.01 g L⁻¹ was used to prepare the calibration curve. The concentration of NH₄⁺–N in sample by equation:
where ρ_N was the concentration of NH₄⁺–N in sample. A_s and A_b were the absorbance of sample and control, respectively. a and b were the intercept and slope of the calibration curve, respectively. V was the volume of sample. Controls were conducted with the re-distilled water.

2.3 DNA extraction

Microbial genomic DNA was extracted by the E.Z.N.A.® Soil DNA Kit (Omega Bio-Tek, Norcross, GA, USA), following the manufacturer protocol. DNA quality and quantity were measured by NanoDrop 1000 Spectrophotometer (Thermo Fisher Scientific, USA) and agarose gel electrophoresis.

2.4 PCR-DGGE analysis

PCR amplification of bacterial 16S rRNA genes was performed by a nested PCR approach. 16S rRNA gene was amplified using the primers 27F/1492R¹⁸ in the first PCR round and a secondary touchdown PCR round was carried out on the primary PCR product by using the universal bacterial primers GC338F/518R,¹⁹ targeting 16S rRNA V3 variable region.²⁰ The PCR reaction mixes was performed in a final volume of 50 μL, consisting of 25 μL of Premix ExTaq (Takara, Dalian, China), 2 μL of each primer (10 μM) and 2 μL of DNA templates. PCR amplification was carried out using an Eppendorf Mastercycler Gradient (Eppendorf, Hamburg, Germany) with the following touch-down temperature cycling program: initial denaturation at 98 °C for 10 min, 20 touchdown cycles of denaturation at 98 °C for 10 s, then annealing at 65 °C (with the temperature decreasing 0.5 °C per cycle) for 30 s, and extension at 72 °C for 45 s, followed by 10 cycles of 98 °C for 10 s, 55 °C for 30 s and 72 °C for 45 s and a final extension at 72 °C for 10 min before holding at 16 °C. Five microliters of the PCR product was analysed in 1% agarose gel at 120 V for 20 min. Gels were visualized with the Gel Doc™ XR+ Imaging System (Bio-Rad) and only strong bands of suitable size could be used for DGGE analysis.

DGGE was performed with the Dcode System apparatus (Bio-Rad). PCR amplicons were loaded onto 8% (w/v) polyacrylamide gels with 30% to 55% denaturing gradient (100% denaturant corresponds to 7 M urea and 40% (v/v) deionized formamide) in 0.5× TAE (20 mM Tris, 10 mM Na-acetate, 0.5 mM Na₂EDTA [pH 7.4]) and subjected to an initial electrophoresis at 50 V for 20 min before a 4.5 h continuous run at 165 V and 60 °C, followed by 15 min silver staining. Thereafter, the DGGE gel image was acquired by HP ScanJet G4050 scanner (Hewlett Packard, USA).

2.5 DGGE patterns and phylogenetic analysis

The software Quantity One (Version 4.6.2, BioRad, USA) was used to analyze the migration and intensity of DGGE bands. The bands with identical migration position were regarded as the same species. Sample patterns were normalized to ensure the accuracy of the results. Metrix contained bands intensity data was exported from Quantity One and used for PCA by Conoco for Windows version 4.5 software. Cluster analysis was performed using Unweighted Pair Group Method with Arithmetic Mean (UPGMA) method. The excision and identification of DGGE bands were conducted as previously descripted.¹⁹ Sequences were compared with the GenBank database by using BLAST program. Phylogenetic analysis was carried out with MEGA 6.0 through the neighbor joining tree with bootstrap value of 1000. The identified sequences were archived at GenBank under accession numbers (KX223781–KX223812, KX223814).

2.6 Illumina HiSeq sequencing of 16S rRNA and sequence analysis

The microbial community of the different sections was identified by amplifying and analyzing the V4 region of 16S rRNA from the genome DNA as previous study.¹⁹ The amplicons sequencing were determined on the Illumina HiSeq 2500 PE250 platform at Novogene Bioinformatics Technology (Beijing, China). Raw HiSeq sequencing data have been deposited to the NCBI sequence read archive (SRA) database with accession number of SRP075188.

Samples in this study were individually barcoded to enable multiplex sequencing. Raw HiSeq sequencing data obtained from this study were filtered by appraising data quality and trimming off the primers and barcodes to obtain high-quality tags according to the QIIME (V1.7.0, http://qiime.org/scripts/split_libraries_fastq.html) pipeline. The chimera sequences in the tags were detected by (UCHIME Algorithm, http://www.drive5.com/usearch/manual/uchime_algo.html). After chimera removal, the effective tags finally obtained.

2.7 Data analyses

The effective tags from high throughput sequencing in this study were clustered into operational taxonomic units (OTUs) using Uparse (Uparse v7.0.1001, http://drive5.com/uparse/), with setting 0.03 distance limit (equivalent to 97% similarity). The taxonomic information of representative sequences for each OTU were annotated via GreenGene Database (http://greengenes.lbl.gov/cgi-bin/nph-index.cgi) and RDP Classifier (Version 2.2, http://sourceforge.net/projects/rdp-classifier/) with a set confidence threshold of 80%. From the cluster file, the alpha diversity indice containing Shannon–Wiener diversity index, Chao 1 richness estimator, the ACE estimator, Simpson diversity index and the Good's coverage were calculated in QIIME (Version 1.7.0) for each sample. Rarefaction curve and Rank abundance curve were generated by R software (Version 2.15.3). For beta diversity, UniFrac distance and unweighted pair-group method with arithmetic means (UPGMA) clustering were also calculated by QIIME (Version 1.7.0). Clustering analysis, PCoA analysis and weighted UniFrac distance metrics were conducted based on UniFrac metrics.

2.8 Quantitative real-time PCR

The abundance of bacterial 16S rRNA genes was appraised by qPCR with universal primers 338F/518R.¹⁹ qPCR was carried out in 25 μL volumes using SYBR® Premix Ex Taq™ II (Takara, Dalian, China) on an ABI 7500 Real Time PCR System (Applied Biosystems). The 25 μL reaction mixtures consisted of 12.5 μL SYBR® Premix Ex Taq™ II (Takara, Dalian, China), 0.5 μL of each primer (10 μM), 0.5 μL of Rox Dye II, 2 μL of genomic DNA and 9 μL ddH₂O. Cycling condition for 16S rRNA gene quantification was as follows: 95 °C for 30 s, followed by 40 cycles of 10 s at 95 °C, 34 s at 55 °C and 30 s at 72 °C. Fluorescence detection was performed at the annealing step. Standards were made from 10-fold dilutions of linearized clone QSD1 (KX223813) containing bacteria 16S rRNA genes in the range of 5.5 × 10⁵ to 5.5 × 10¹¹ copies. Three negative control was included in qPCR assay (template DNA replaced by double-distilled water), and all standard reactions and samples were carried out in triplicate. The efficiency of the qPCR assay was 100%, and the R² value for standard curve line exceeded 0.999.

2.9 Statistics analysis

The correlations between operational parameters and microbial abundance were descripted by SRCC,²¹ which was calculated by SPSS Statistics 22.0 (IBM, USA), and RDA was carried out using Conoco for Windows version 4.5 software to assess the relationship between the features of wastewater and microbial diversity.

3. Results and discussion

3.1 Analysis of DGGE band patterns

DGGE patterns of the total microbial community in each tanks of tannery wastewater treatment process were shown in Fig. 2a. A total of 35 prominent bands were excised from DGGE, amplified, cloned and sequenced which indicated that obviously microbial dynamics were observed during the integrated tannery wastewater treatment process. For example, band 30, which showed 99% similarity with Frigovirgula sp. (see ESI, Fig. S1†) appeared with high intensity in PST and GC while was dim in HT and AT, and subsequently disappeared in the remaining four tanks. On the contrary, band 24 was detected with high concentration in these samples except for AT. Furthermore, bands 18 and 19 showed higher intensity than other tanks and bands 31, 32 and 33 were only appeared in AT.

Fig. 2 Analysis of microbial diversity by DGGE method. (a) DGGE profiles of the microbial community in samples at different stages; (b) principal components analysis (PCA) of DGGE profiles; (c) clustering analysis of the DGGE patterns of the 8 samples.

A total of 4 phyla were detected in all tanks by DGGE analysis (see ESI, Fig. S1†). Proteobacteria was the dominant phylum on average, followed by Bacteroidetes, Firmicutes and Chlorobi, which were similar to the previous study on printing and dyeing wastewater treatment system.²² However, Firmicutes, following Proteobacteria, were more abundant than Bacteroidetes in PST and GC. In addition, Chlorobi showed higher abundance than Bacteroidetes in FST. Based on the DGGE fingerprints, PCA analysis of microbial community composition in whole treatment process were presented in Fig. 2b. PC1 and PC2 explained 30.9 and 26.9% of the variation, respectively. The samples collected from PST, GC and HT gathered into one group and CWT shared similarity on microbial composition with FST. This probably reflected the metabolic stability of bacteria in the pre-treatment stage. Moreover, AT, OT and SST were individually divided as a single branch. The results of cluster analysis by UPGMA method were in accordance with that of PCA analysis (Fig. 2b and c). Overall, the microbial community and composition of each tanks in the process were complex and distinct.

3.2 Diversity of microbial community

After filtering the low quality and trimming the primers, adapters and barcodes, a total of 33 770 to 64 655 effective sequence tags of 16S rRNA gene with an average length of about 250 bp for eight samples were obtained using Illumina HiSeq sequencing. Rarefaction curves of these samples tended to reach saturation, which indicated that microbial communities in samples were well-represented and the sequencing data were reasonable (see ESI, Fig. S2†). Furthermore, a 97% similarity cut-off was used to calculate the numbers of OTUs, Shannon, Chao 1, ACE, Simpson and Good's coverage for the eight samples in the downstream analyses at the same sequencing depth (Table S1†). Similar to the numbers of microbial OTUs (3% cutoff) in water purification plant,²³ the OTU numbers of these samples ranged from 995 to 2361, but much less than that in anaerobic sludge¹² and indicated that the sample from CWT had the richest diversity, followed by those from PST and OT. Additionally, the Good's coverage of these samples ranged from 97.8% to 99.2%, indicating that the sequence libraries constructed in this study could represent the diversity of microbial community. The OTU number, Chao 1 and ACE demonstrated that the richness values varied by 1.1–3.6 times among these samples from different tanks. Overall, the bacterial diversity among the integrated process were significantly different and samples in the A/O process had much higher diversity than that in pre-treatment stage (Table S1, ESI Fig. S2†).

Based on the weighted UniFrac distance metrics, PCoA was carried out to appraise the similarities of these samples from different tanks. PCoA analysis, with maximum variation of 44.22% (PC1) and 23.8% (PC2), demonstrated that samples from SST and AT, HT and GC tended to cluster together, respectively (Fig. 3a). Nevertheless, sample from PST shared slight similarity on microbial compositions with OT, CWT as well as FST, and obviously distinguished from other tanks. Phylogeny-based UniFrac analysis, based on abundances of phyla, revealed that bacterial communities in these samples could be clustered into three groups: (1) group I contains the samples from PST, OT, CWT and FST; (2) group II contains the two samples from AT and SST; (3) group III are the samples from GC and HT, which resulted in the same trend with PCoA analysis (Fig. 3b). Therefore, it was obvious that the bacterial diversity among pre-treatment stage, anoxic stage and oxic stage were entirely distinct and we hypothesized that the significant distinctions might be due to the characteristics of wastewater as well as the operational parameters.

Fig. 3 Multiplex analysis of microbial communities using weighted-UniFrac from high throughput sequencing. (a) Principal coordinate analysis of the samples; (b) phylogenetic tree.

3.3 Bacterial community analysis at high taxonomic levels

In present study, a total of 47 phyla across all samples were identified and the major phyla for each sample were shown in Fig. 4. The predominant phylum was Proteobacteria with 46.07% on average, followed by Bacteroidetes (17.28%), Firmicutes (11.13%), Actinobacteria (5.92%) and Chloroflexi (2.24%). These were consistent with the microbial composition in activated sludge treating tannery wastewater²⁴ and similar to the microbial structures in activated sludge taken from the aeration tanks of conventional activated sludge process in USA which mainly treated municipal wastewater.²⁵ Members of these phyla generally were dominant in industrial wastewater treatment plants, which could be explained by their spore-forming ability which could ensure survival under stressful conditions.^12,16 However, the results were different from the study on bacterial communities in upflow anaerobic filter sludge digesting the mixture of olive mill and abattoir wastewaters,²⁶ which demonstrated that Firmicutes was dominant phyla and represented 54.72% and 41.2% in mesophilic anaerobic sludge and thermophilic anaerobic sludge, respectively. In addition, tags affiliated to Planctomycetes, Gemmatimonadetes, Nitrospirae, Acidobacteria and TM7 accounted for a small proportion.

Fig. 4 Relative abundance of the major bacterial community at phylum level in each sample. Taxa represented occurred at higher than 1% abundance in at least one sample. Sequences with their maximum abundance less than 1% at any sample were summarized in the artificial group “others”.

Distinctions in microbial community composition among samples were discovered. Previous reports indicated that anaerobic and aerobic reactors significantly harbored different bacterial communities due to the variation of abiotic factors.¹⁶ In this study, sample of OT displayed a considerably higher proportion of [Thermi] (10.07%) and Planctomycetes (4.45%) than that detected in AT. Additionally, it has been reported that Planctomycetes were strongly involved in the degradation of dissolved organic matter (DOM) in nutrient-poor waters with supplementary N²⁷ and anaerobic ammonium oxidation (anammox) process.²⁸ Hence it could be concluded that Planctomycetes played a significant role in tannery wastewater treatment. In addition, the sample from AT were comprised of Proteobacteria (42.59%), Bacteroidetes (25.91%), Actinobacteria (17.72%) and Firmicutes (5.84%), while the consortia of OT were significantly multiplex, with a superiority of Proteobacteria (52.17%), Bacteroidetes (3.26%), Planctomycetes (4.45%), Gemmatimonadetes (2.03%) and Chloroflexi (5.83%). The proportion of Bacteroidetes and Actinobacteria in this study agreed with the results revealed by previous study,²⁹ which demonstrated that Bacteroidetes were more abundant at low dissolved oxygen level. However, Bacteroidetes (31.79%) and Firmicutes (26.43%) were much more abundant in HT and GC, respectively, which was consistent with that Bacteroidetes made major contribution to degrading protein and DOM.²⁹ Nevertheless, it should be noted that 6.29% of total sequences were not classified at phylum level.

Though Proteobacteria (mostly Alpha-, Beta-, Gamma- and Epsilon-proteobacteria) predominated in all samples, they had little in common in terms of composition at class level. Alpha-proteobacteria had been considered as one of the dominant group in tanneries effluent and industrial wastewater.³⁰ Analogously, Alpha-proteobacteria (16.67–27.5%, averaging 20.64%) was the predominant classes in samples from the tanks AT, OT, FST and CWT. However, this was distinct with the previous study using 16S clone library,¹⁶ which suggested that Beta- and Delta-proteobacteria within Proteobacteria were dominant in aerobic and anaerobic reactors, respectively. In addition, Gamma-proteobacteria (18.16–21.6%, averaging 19.88%) dominated in HT and PST. Moreover, compared with other tanks, Beta-proteobacteria (19.91%) and Epsilon-proteobacteria (29.37%) had a higher proportion in SST and GC, respectively. Nevertheless, Delta-proteobacteria (3.8%) appeared at low percentage in all samples on average. Except for the subdivisions of Proteobacteria, other major classes including Bacilli (3.19%), Deinococci (2.04%), Actinobacteria (4.04%), Clostridia (7.86%), Flavobacteriia (5.07%), TM7-3 (1.01%) and Bacteroidia (8.35%), were also detected on average.

The 35 most abundant genera in each sample were selected and compared with their abundances in other samples to assess the dynamic evolution in microbial community composition during the wastewater treatment process (Fig. S3 in ESI†). It can be seen that there was a distinct fluctuation on the microbial composition. For example, Arcobacter, which was generally considered as human pathogens due to their frequent isolation from ill animals, chicken carcasses, and humans with enteritis,³¹ was the most dominant genus in the GC (28.75%), HT (8.16%), PST (11.22%), whereas the samples from other tanks contained relatively less (0.028–1.37%). Therefore, we inferred that these pathogens might derive from the raw materials of leather production. Meanwhile, the genus Thioalkalimicrobium, which was known to oxidize sulfide and thiosulfate to sulfate at high rates in alkaline conditions,³² were only abundant in GC, whereas its abundances were <1% in other samples. Furthermore, some genera, such as Tissierella_Soehngenia, Lactococcus, Pseudomonas, Lysobacter, Proteiniclasticum, were the major (abundance > 1%) genera in at least two samples. Members of Pseudomonas were capable of degrading the azo dyes and hexavalent chromium (Cr VI)³³ and heterotrophic nitrification and aerobic denitrification.³⁴ Previous study also detected apr and npr genes encoding metallopeptidases in Pseudomonas fluorescens and Bacillus biotypes, respectively³⁵ so that the two biotypes could convert organic nitrogen into NH₄⁺–N which called ammonification for subsequent nitrification and denitrification resulting in the rise of NH₄⁺–N in PST which was in accordance with the high abundance of Pseudomonas and Bacillus in PST. In addition, members of Nitrosococcus, such as Nitrosococcus oceani, had previously been reported as ammonia-oxidizing bacteria³⁶ and Nitrosococcus appeared distinctly higher abundance in OT (1.973%) than that in AT (0.159%) and SST (0.013%) in this study which was consistent with 89.78% NH₄⁺–N removal in OT. Moreover, genera including Nitrospira and Rhodococcus were mainly detected in PST. Thus, most of above mentioned genera might make a significant contribution in reducing pollutants in tannery wastewater.

In this study, two archaeal phyla, Euryarchaeota and Thaumarchaeota were also detected. The species affiliated to vadinCA11 and Candidatus Nitrososphaera was found with differences in numbers for each samples. This suggested that they could pose a potential ability on contaminants removal. Previous studies reported that the species in the genus Nitrososphaera, such as Nitrososphaera viennensis and Nitrososphaera gargensis were known as an aerobic and mesophilic, ammonia-oxidizing archaeon from soil³⁷ and had the ability to transform mianserin and ranitidine.³⁸ Thus, it was worth to note that the potential correlation between bacteria and archaea in wastewater treatment process.

3.4 Comparison of high throughput sequencing and PCR-DGGE analysis

In present study, both PCR-DGGE and high throughput sequencing revealed the microbial diversity and composition in each step of the integrated wastewater treatment process, and the predominant consortia in the all samples were of similarity at phylum level. Proteobacteria were demonstrated as the dominant phylum in the integrated process by the two methods, while slight distinction appeared between the two methods. For instance, Petrimonas as well as Thauera were only detected to some extend in the process by PCR-DGGE method (see ESI Fig. S1†) but not appeared in high throughput sequencing results (data not shown). This could be attributed to the primer variability and PCR parameters.³⁹ Nevertheless, compared with PCR-DGGE analysis, high throughput sequencing explained accurate and detailed analysis of bacterial community due to the appearance of some minor species in HiSeq sequencing results. For example, consortia affiliated to Azoarcus, Rhodococcus, Pseudomonas, Nitrospira and Nitrosococcus were only detected by high throughput sequencing (see ESI Fig. S3†).

The propinquity and differentiation of DGGE fingerprints in different steps were graphically displayed by PCA and cluster analysis (Fig. 2b and c). Cluster analysis based on PCR-DGGE and HiSeq sequencing data indicated that samples from CWT and HT clustered together with FST and GC, respectively (Fig. 2c and 3b). However, microbial diversity of PST shared similarity with those on the group consisted of HT and GC and the group composed of OT, CWT and FST as shown by PCR-DGGE (Fig. 2c) and high throughput sequencing (Fig. 3b). In addition, the results of PCA (PCoA) analysis described the same trends with cluster analysis (Fig. 2b and 3a).

Overall, both the two methods complemented each other and revealed the microbial diversity in tannery wastewater treatment process. What is noteworthy is that PCR-DGGE method was sufficient to detect the abundant species but deficient in identifying low-abundance communities, whereas high throughput sequencing was more accurate in analyzing integrated community composition and investigating the slight composition variation.

3.5 Quantitative analysis of total bacteria

qPCR was further applied to validate the abundance of total bacteria 16S rRNA gene in different tanks. Results of this indicated the variation of 16S rRNA gene abundances of the total microbes in different stages of wastewater treatment process (Fig. 5). The copy numbers of the total bacterial 16S rRNA gene varied from 4.22 × 10⁸ copies per mL to 4.93 × 10¹¹ copies per mL, which was in agreement with previous study⁴⁰ but higher than that in effluent of treated municipal wastewater.⁴¹ The copy numbers of 16S rRNA in HT, AT and OT were in the same order of magnitude, with nearly two to three order of magnitude higher than that of PST, GC and SST, CWT, respectively. As shown in Fig. 5, the bacterial copy numbers noticeably increased in HT, AT and OT, which were mainly steps in biological wastewater treatment process and these were in accordance with that bacteria played an important role in wastewater treatment.⁴² Consequently, high abundance of bacteria in HT might guarantee the degradation of pollutants with complicated structure resulting in the temporary rise on COD concentration.

Fig. 5 Abundance of total bacterial 16S rRNA gene in each tank of the integrated process. Error bars represent standard deviation calculated from three independent assays.

3.6 Correlations of environmental parameters and bacterial communities

Environmental conditions generally had significant impacts on microbial ecology.⁴³ Multivariate analyses were thoroughly conducted to assess the potential relationships between the microbial community composition and the environmental factors. SRCC and redundancy analysis plots were calculated from the bacterial profiles at phylum level and the results were shown in Table S2† and Fig. 6, respectively. Results of SRCC indicated that the abundance of [Thermi] was negatively correlated to the concentration of COD (SRCC = −0.714, p = 0.047), and positively correlation with NH₄⁺–N removal efficiency (SRCC = 0.786, p = 0.036) (Table S2†). Although the presence of [Thermi] had been well described,¹² the links between the appearance of [Thermi] and the concentration of COD and NH₄⁺–N removal need further study.

Fig. 6 RDA correlation triplot depicting the abundance and distribution of microorganisms in relation to physicochemical characteristics of wastewater (RE: removal efficiency).

RDA analysis, with constrained ordination, was used to examine how much of the variance in environmental parameters could explain the variance in microbial composition.⁴⁴ As shown in Fig. 6, the first two canonical axes explained 37.7% and 22.6% of the variance in data, respectively. Moreover, it could be found that microbial structure in GC had similarity with this in HT, which was identical with the cluster analysis (Fig. 2c and 3b). In addition, the results demonstrated that samples from PST, OT and CWT had higher microbial diversity than these of other tanks, which were in accordance with the taxonomic distribution of microbial community composition (Fig. 4). Meanwhile, the RDA analysis also showed that ammonia concentration was one of the main variables affecting the bacterial diversity, which was in accordance with the study on bacterial diversity in Dongjiang River.⁴⁵

Triplot analysis revealed that Bacteroidetes, Firmicutes and Actinobacteria were sensitive to the COD concentration, pH and COD removal efficiency respectively. Therefore, we deduced that the COD removal efficiency could be affected by the abundance of Actinobacteria. Moreover, previous study also demonstrated that the dominant bacterial groups changed from Beta-proteobacteria and Firmicutes to Acidobacteria and Alpha-proteobacteria due to pH variation.⁴⁶ So we deduced that the abundance of Firmicutes had positive correlation with pH and this result was consistent with the results of SRCC, which suggested that the abundance of Firmicutes might be influenced by pH variation. It was interesting that the abundance of Nitrospirae presented an opposite trend with COD value (Fig. 6, Table S2†). This was in accordance with the previous study and it was postulated that this might be explained by the stronger competition between nitrifying bacteria and strictly heterotrophic bacteria.²⁴ In addition, Planctomycetes had positive correlation with Chloroflexi. In the study on syntrophic butyrate oxidation, detection of Chloroflexi and Planctomycetes by ¹³C isotope probing revealed that they possibly assimilated [¹³C] acetate produced intermediately during the butyrate degradation.⁴⁷ Therefore, we inferred that the two phyla made contribution to removing COD and the coexistence and interaction between the two phyla for treating wastewater containing organic matter in anoxic environments need further research.

Overall, the correlation between the abundance and diversity of microorganisms and environmental parameters were systematically analyzed, but the coexistence and interaction among the microbes need further study.

4. Conclusion

Combination of PCR-DGGE and high throughput sequencing demonstrated that Proteobacteria, Bacteroidetes, Firmicutes and Actinobacteria composed the major microbial consortia in the integrated tannery wastewater treatment process. Quantitative PCR revealed that the abundance of 16S rRNA gene in AT and OT were significantly higher than that in other tanks. Multivariate analyses suggested that the bacterial diversity could be affected by physiochemical characteristics of wastewater and pre-treatment, anoxic stage as well as oxic stage shared different bacterial compositions. In addition, the relationship between NH₄⁺–N removal efficiency and [Thermi] should be investigated in the further study.

Acknowledgements

This study is financially supported by the National Natural Foundation of China (grant 31271924) and the National Training Program of Innovation and Entrepreneurship For Undergraduates of China (201410561102). We thank our institutions and lab colleagues for support.

References

1. K. Cooman, M. Gajardo, J. Nieto, C. Bornhardt and G. Vidal, Environ. Toxicol., 2003, 18, 45–51.
2. R. Ganesh, G. Balaji and R. Ramanujam, Bioresour. Technol., 2006, 97, 1815–1821.
3. M. Zamparas, M. Drosos, Y. Georgiou, Y. Deligiannakis and I. Zacharias, Chem. Eng. J., 2013, 225, 43–51.
4. L. Ye and T. Zhang, Appl. Microbiol. Biotechnol., 2013, 97, 2681–2690.
5. G. Vidal, J. Nieto, K. Cooman, M. Gajardo and C. Bornhardt, J. Hazard. Mater., 2004, 112, 143–149.
6. L. Szpyrkowicz, S. N. Kaul, R. N. Neti and S. Satyanarayan, Water Res., 2005, 39, 1601–1613.
7. C. Di Iaconi, A. Lopez, R. Ramadori, A. Di Pinto and R. Passino, Water Res., 2002, 36, 2205–2214.
8. M. Zhang, J. H. Tay, Y. Qian and X. S. Gu, J. Environ. Eng., 1997, 123, 876–883.
9. O. Canals, H. Salvadó, M. Auset, C. Hernández and J. Malfeito, Water Res., 2013, 47, 3141–3150.
10. K. Rasool, D. H. Ahn and D. S. Lee, Bioresour. Technol., 2014, 162, 373–378.
11. M. Wagner, A. Loy, R. Nogueira, U. Purkhold, N. Lee and H. Daims, Antonie van Leeuwenhoek, 2002, 81, 665–680.
12. D. Shu, Y. He, H. Yue and Q. Wang, Bioresour. Technol., 2015, 186, 163–172.
13. N. Boon, W. De Windt, W. Verstraete and E. M. Top, FEMS Microbiol. Ecol., 2002, 39, 101–112.
14. F. M. Ibarbalz, E. L. Figuerola and L. Erijman, Water Res., 2013, 47, 3854–3864.
15. Z. Wang, X.-X. Zhang, X. Lu, B. Liu, Y. Li, C. Long and A. Li, PLoS One, 2014, 9, e113603.
16. A. F. Desta, F. Assefa, S. Leta, F. Stomeo, M. Wamalwa, M. Njahira and D. Appolinaire, PLoS One, 2015, 10, 443–461.
17. N. Chinese, Water and wastewater monitoring methods, Chinese Environmental Science Publishing House, Beijing, 2002.
18. D. Lane, in Nucleic acid techniques in bacterial systematics, ed. E. Stackebrandt and M. Goodfellow, John Wiley & Sons, Chichester, United Kingdom, 1991, pp. 125–175.
19. P. Li, H. Liang, W.-T. Lin, F. Feng and L. Luo, Appl. Environ. Microbiol., 2015, 81, 5144–5156.
20. G. Muyzer, E. C. De Waal and A. G. Uitterlinden, Appl. Environ. Microbiol., 1993, 59, 695–700.
21. M. Mußmann, I. Brito, A. Pitcher, J. S. S. Damsté, R. Hatzenpichler, A. Richter, J. L. Nielsen, P. H. Nielsen, A. Müller and H. Daims, Proc. Natl. Acad. Sci. U. S. A., 2011, 108, 16771–16776.
22. Q. Yang, J. Wang, H. Wang, X. Chen, S. Ren, X. Li, Y. Xu, H. Zhang and X. Li, Bioresour. Technol., 2012, 117, 155–163.
23. W. Lin, Z. Yu, H. Zhang and I. P. Thompson, Water Res., 2014, 52, 218–230.
24. C. Giordano, V. Boscaro, A. Mannucci, G. Munz, G. Mori and C. Vannini, SIMGBM, Ravenna, September 2015.
25. T. Zhang, M.-F. Shao and L. Ye, ISME J., 2012, 6, 1137–1147.
26. H. Gannoun, I. Omri, R. Chouari, E. Khelifi, J.-J. Godon, M. Hamdi, A. Sghir and H. Bouallagui, Bioresour. Technol., 2016, 201, 337–346.
27. R. D. Tadonléké, FEMS Microbiol. Ecol., 2007, 59, 543–555.
28. C. J. Schubert, E. Durisch-Kaiser, B. Wehrli, B. Thamdrup, P. Lam and M. M. Kuypers, Environ. Microbiol., 2006, 8, 1857–1863.
29. S. J. Ma, L. L. Ding, H. Huang, J. J. Geng, K. Xu, Y. Zhang and H. Q. Ren, Bioresour. Technol., 2016, 214, 645–652.
30. O. Lefebvre, N. Vasudevan, K. Thanasekaran, R. Moletta and J. J. Godon, Extremophiles, 2006, 10, 505–513.
31. C. Phillips, Food Control, 2001, 12, 1–6.
32. D. Y. Sorokin, A. M. Lysenko, L. L. Mityushina, T. P. Tourova, B. E. Jones, F. A. Rainey, L. A. Robertson and G. J. Kuenen, Int. J. Syst. Evol. Microbiol., 2001, 51, 565–580.
33. S. Mahmood, A. Khalid, T. Mahmood, M. Arshad and R. Ahmad, J. Chem. Technol. Biotechnol., 2013, 88, 1506–1513.
34. S. Yao, J. Ni, T. Ma and C. Li, Bioresour. Technol., 2013, 139, 80–86.
35. H.-J. Bach, A. Hartmann, M. Schloter and J. Munch, J. Microbiol. Methods, 2001, 44, 173–182.
36. M. G. Klotz, D. J. Arp, P. S. Chain, A. F. El-Sheikh, L. J. Hauser, N. G. Hommes, F. W. Larimer, S. A. Malfatti, J. M. Norton and A. T. Poret-Peterson, Appl. Environ. Microbiol., 2006, 72, 6299–6315.
37. M. Stieglmeier, A. Klingl, R. J. Alves, K.-M. R. Simon, M. Melcher, N. Leisch and C. Schleper, Int. J. Syst. Evol. Microbiol., 2014, 64, 2738–2752.
38. Y. Men, P. Han, D. E. Helbling, N. Jehmlich, C. Herbold, R. Gulde, A. Onnis-Hayden, A. Z. Gu, D. R. Johnson and M. Wagner, Environ. Sci. Technol., 2016, 50, 4682–4692.
39. R. Sipos, A. J. Székely, M. Palatinszky, S. Révész, K. Márialigeti and M. Nikolausz, FEMS Microbiol. Ecol., 2007, 60, 341–350.
40. H. M. Jang, J. H. Ha, M.-S. Kim, J.-O. Kim, Y. M. Kim and J. M. Park, Water Res., 2016, 99, 140–148.
41. S. Liu, Y. Luo and L. Huang, Appl. Microbiol. Biotechnol., 2016, 1–10.
42. T. Osaka, K. Shirotani, S. Yoshie and S. Tsuneda, Water Res., 2008, 42, 3709–3718.
43. C. Yang, W. Zhang, R. Liu, Q. Li, B. Li, S. Wang, C. Song, C. Qiao and A. Mulchandani, Environ. Sci. Technol., 2011, 45, 7408–7415.
44. L. H. Hagen, V. Vivekanand, R. Linjordet, P. B. Pope, V. G. Eijsink and S. J. Horn, Bioresour. Technol., 2014, 171, 350–359.
45. Z. Liu, S. Huang, G. Sun, Z. Xu and M. Xu, FEMS Microbiol. Ecol., 2012, 80, 30–44.
46. S. K. Goffredi, A. H. Kantor and W. T. Woodside, Microb. Ecol., 2011, 61, 529–542.
47. P. Liu, Q. Qiu and Y. Lu, Appl. Environ. Microbiol., 2011, 77, 3884–3887.

---

Footnote

† Electronic supplementary information (ESI) available. See DOI: 10.1039/c6ra19603a

---