Enhanced Acetogenesis of Waste Activated Sludge by Conditioning with Processed Organic Wastes in Co-Fermentation: Kinetics, Performance and Microbial Response

: Aimed at the low ratio of carbon and nitrogen (C / N, approximately 7 / 1) of waste activated sludge (WAS), which would inhibit the acetogenesis process during anaerobic fermentation, this study introduced three brewing wastes, including vinegar (VR), stillage (SR) and soy sauce (SSR) residues, to promote acetogenesis by co-fermenting with WAS. Results showed that di ﬀ erent brewing wastes contributed di ﬀ erently to the volatile fatty acids (VFAs) yield. The best performance was observed with SSR (4517 ± 367 mg COD / L), particularly rich in C 2 –C 3 VFAs, corresponding to 40% and 52% higher concentrations than with SR and VR, respectively. Meanwhile, the hydrolysis rate constant peaked at 0.0059 h − 1 in the SSR test, compared to the sole WAS test (0.0018 h − 1 ). Furthermore, canonical correlation analysis reﬂected that the functional consortia, known to ferment saccharides / amino acids into C 2 –C 3 VFAs (i.e., Proteiniclasticum, Petrimonas, Cloacibacillus and Gemmobacter ), was related to the characteristics of the feedstock.


Introduction
Huge quantities of sludge, i.e., primary sludge and waste activated sludge (WAS), are produced from municipal wastewater treatment plants (WWTPs), thus displaying a pressing challenge inherent to traditional activated sludge management. In China, 6.25 million tons of sludge were produced in 2013, while the average annual growth rate increased by 13%, from 2007 to 2013 [1]. Anaerobic digestion (AD), known to be a promising method for renewable energy/bio-metabolites recovery from biomass, has gained worldwide attentions and has been strongly recommended. On the other hand, volatile fatty acids (VFAs), as an important intermediate product of AD, are a promising substrate for many different bioprocesses, e.g., biopolymer production [2], bioenergy generation [3] as well as biological nutrient removal [4]. The best way to utilize VFAs is without any doubt related to their composition, especially for the low-molecular-weight (C 2 and C 3 ) VFAs. Targeting the promotion of the acetogenesis step during WAS digestion is relatively scarce. Therefore, in order to enhance the VFAs yield and/or enrich for a VFA, current studies have been focusing on the optimization of operating

Materials
WAS was collected from a secondary sedimentation tank in the Jinzhong municipal wastewater treatment plant (WWTPs) (Jinzhong City, China). Raw WAS was concentrated by settling at 4 • C for 24 h, prior to testing. The main characteristics are displayed in Table 1. VR, SR and SSR were obtained respectively from Xiaoqi Vinegar Workshop (Qingxu City), Fenjiu Group CO., LTD (Fenyang City) and Jinjia Brewing Factory (Taiyuan City). The dried PWs were then immersed in 2% NaOH solution at 85 • C in a thermostatic water bath. The solid to liquid ratio was 1:10 (g dry weight to mL) and the residence time was 1 h. The main composition of pretreated PWs on a dry basis is reported in Table 2.

Fermentation Setup
Batch experiments were carried out in the 500 mL serum bottles with a sealed lid, in which one port was used to collect gas and another port was equipped with one tube for sampling. Four sets of experiments were conducted, each containing 300 mL of the substrates. Three sets were fed with the mixtures of WAS and one of the PWs, namely VR, SR and SSR (hereinafter referred to as VR, SR and SSR tests). The control test was performed using only concentrated WAS. The PWs was added with 35% PWs to 65% WAS based on VSS concentration. The fermenters were flushed with nitrogen 10-15 min and sealed immediately to guarantee strict anaerobic atmosphere. In order to promote fermentation, all bottles were incubated at 35 ± 1 • C and stirred with 100 rpm. The samples were withdrawn from each reactor at an interval of about 24 h to study the influence of different types of PWs on WAS fermentation.

DNA Extraction and Illumina MiSeq Sequencing
DNA extraction and sequencing were conducted according to a previous study (Wen et al. 2017). The V3-V4 region of the 16S rRNA gene were amplified using 341F and 805R primers. Raw data was purified and quantified, then analyzed in the Illumina MiSeq platform. Raw sequences were deposited in the NCBI Short Read Archive database with the accession No. SRR6179062. Sequences were clustered into operational taxonomic units (OTUs) at identity thresholds of 0.97. A venn diagram with OTUs was used to visualize the similarity and difference between the four fermentation sets. Canoco 4.5 was used to explore the potential interrelation between process performance and key microflora.

Analytical Methods and Statistical Analysis
Samples from each fermenter were centrifuged at 10,000 rpm for 10 min, then filtered through a 0.45 µm cellulose nitrate membrane filter and finally stored at 4 • C. The determination of total suspended solids (TSS), volatile suspended solids (VSS), soluble chemical oxygen demand (SCOD), total chemical oxygen demand (TCOD), NH 4 + -N, PO 4 3--P was conducted according to standard methods for examination of water and wastewater [17]. The lignocellulose contents of PWs was analyzed by a Fiber Analyzer (Foss, DK). Lipid and oil were measured as the previous study [15]. Proteins and carbohydrates were determined as our previous study [16]. Proteins were determined by the bicinchoninine acid assay (BCA)protein assay kit (Sangon, Shanghai, China) and carbohydrates were measured by the phenol-sulfuric method. The activities of α-glucosidase and protease were measured according to Goel et al. [18]. Biogas and VFAs were detected using gas chromatographs (GC) (4890D and 7890, Agilent, USA, respectively). 3D-excitation-emission matrix (EEM) of fermentation filtrates (48 h) was monitored on an F-7000 spectrofluorometer (Hitachi, Tokyo, Japan). The spectra were analyzed using the parallel factor analysis (PARAFAC) algorithm to quantitatively determine the components of fluorophores. For the kinetic characterization, a Monod-type model was used to simulate the hydrolysis constant, according to the following equation [19]: where Y CH4 is the methane yield (mLCH 4 /gVSS), Y CH4 ∞ is the ultimate methane yield (mLCH 4 /gVSS), k h is the hydrolysis constant (h −1 ), and t is the treatment time (h). The infrared spectra were recorded by a FTIR spectrometer (PerkinElmer, USA). The crystallinity indexes (CI) of the biomasses were determined using Nelson and O'connor's method [20] as illustrated in Equation (2): where CI represents the crystallinity indexes; α 1372cm −1 and α 2900cm −1 represents the infrared spectral intensities determined at the specified wavenumber of 1372 cm −1 and 2900 cm −1 .

The Compositions of Raw and Pretreated PWs
As can be seen in Table 2, higher amounts of proteins, lipid and oil were found in the dried SSR (22.6% and 8.4%) than that in the VR (11.2% and 2.4%) and SR (15.5% and 7.0%). By contrast, the proportions of fibers were higher in VR and SR. Compared to cellulose in the three PWs, TA pretreatment removed higher amounts of hemicellulose and lignin. The lignocellulosic structures could be opened under this pretreatment, thereby increasing the biocompatibility of cellulose in the subsequent AD process. Thus, a part of hemicellulose could be decomposed to monosaccharides, used as the easily degradable carbon source for the subsequent AD process. The observed reduction of VS also partially supported this result ( Table 2).
Infrared spectroscopy was conducted to observe the structure difference between PWs. As shown in Figure 1A,B, the transmissivity of absorption peak was clearly distinct. 3300-3500 cm −1 , which stands for stretching vibration of O-H in saccharides or -NH in amide, which showed big differences between the three kinds of PWs. After TA pretreatment, the strong absorption band corresponding to the CH 2 -stretching at -2920 cm −1 became flat ( Figure 1B), suggesting a decrease of the aliphatic constituents. -1720 cm −1 peak ascribed to the acetyl, uronic esters and carboxylic groups of the lignin and/or hemicelluloses [21]. The bands near 1635 cm −1 (aromatic skeletal vibrations) also indicated the varied swelling of different PWs. 1540 cm −1 stretching vibration of the benzene ring, (one of the lignin characteristic peaks) [22], indicated the delignification difference of PWs. Evident intensity reduction after TA pretreatment revealed the degradability of hemicellulose and lignin. Additionally, 1163 cm −1 (C-O-C asymmetric) and 1075 cm −1 (C-O stretch) vibrations of cellulose and hemicellulose, showed much more difference than other peaks, which changed slightly for VR and SSR, but decreased for SR after TA pretreatments. cm −1 , which stands for stretching vibration of O-H in saccharides or -NH in amide, which showed 160 big differences between the three kinds of PWs. After TA pretreatment, the strong absorption band 161 corresponding to the CH2-stretching at -2920 cm -1 became flat ( Figure 1B), suggesting a decrease of 162 the aliphatic constituents. -1720 cm -1 peak ascribed to the acetyl, uronic esters and carboxylic groups 163 of the lignin and/or hemicelluloses [21]. The bands near 1635 cm -1 (aromatic skeletal vibrations) also 164 indicated the varied swelling of different PWs. 1540 cm -1 stretching vibration of the benzene ring,

165
(one of the lignin characteristic peaks) [22], indicated the delignification difference of PWs. Evident 166 intensity reduction after TA pretreatment revealed the degradability of hemicellulose and lignin.

167
Additionally, 1163 cm -1 (C-O-C asymmetric) and 1075 cm -1 (C-O stretch) vibrations of cellulose and 168 hemicellulose, showed much more difference than other peaks, which changed slightly for VR and 169 SSR, but decreased for SR after TA pretreatments.

Effect of Different Feedstock on VFAs Production and Composition
As shown in Figure 2A, the VFAs concentration maximized at 96 h and then decreased over time in all tests. The maximum VFAs concentrations were 3105 ± 35, 3144 ± 69, 4517 ± 367 mg COD/L in VR, SR and SSR test (96 h), respectively, while it was 1083 ± 40 mg COD/L in the sole WAS. As can be seen in Figure 2B, representing the VFAs composition, most reactors reached a plateau of total VFAs in 96 h. The top three individual VFAs produced in the control test were acetate acid, propionic acid and iso-valeric acid (18%, 31% and 27%, respectively), while, acetate acid, propionic acid and iso-butyric acid represented the top three VFAs in the co-fermentation tests (with different percentages). More importantly, acetate acid content was clearly promoted, going from 18% in the raw WAS test to 34%, 30% and 48% in the VR, SR and SSR tests, respectively. Likewise, the corresponding propionic acid concentrations were also increased with PWs conditioning (compared with the control). The production of VFAs was mainly derived from the consumption of carbohydrates and proteins by acidogens, the increase of soluble organics by PWs addition especially SSR further enhanced both hydrolysis and acidification processes.

Time-Course Profiles of Soluble Organics in Different Fermenters
Typically, the concentration of soluble organics in the fermentation system depends on the rates of hydrolysis and acidification. By virtue of pretreatment, soluble proteins and carbohydrates all sharply peaked at 24 h ( Figure 2C,D). In conjunction with the subsequent acidification process, these released soluble organics were rapidly consumed by acidogenic bacteria, and then stabilized at different levels. We assumed that the rates of soluble organics production (hydrolysis step) equaled that of consumption (acidification step). Taking the consumption of proteins into consideration, 399 ± 45, 398 ± 9, 619 ± 18 mg COD/L proteins were consumed in the VR, SR and SSR tests during the fermentation time of 24 h to 120 h fermentation. In comparison, a notable increase of soluble proteins was observed during the sole WAS fermentation, thus implying that the production rate was bigger than the consumption. Clearly, PWs conditioning improved the consumption of organics during co-fermenting with WAS.   (Table S1). Based on the Shannon index, bacterial communities in

Key Functional Microbiome Analysis
Flat rarefaction curves ( Figure S1) for all libraries indicated that sampling of the microbial community was effective. α-diversity was assessed to show the microbial diversity of the communities in different samples (Table S1). Based on the Shannon index, bacterial communities in the Control (5.36) were much more diverse than those in PWs conditioning tests (5.11, 5.25 and 4.75 for VR, SR and SSR, respectively). Simpson index showed the same phenomenon, with SSR having the highest value (0.036). In other words, a reduction in microbial diversity happened in the WAS treatment after co-digesting with PWs. The nonproportional venn diagram showed that the shared OTUs (by the four groups) were 498, only 10.7% of total OTUs ( Figure S2A). They mainly belonged to the phyla Protebacteria (32.5%), Firmicutes (15.5%) and Chloroflex (9.8%) ( Figure S2B). Principal component analysis (PCoA) of β-diversity, based on unweighted UniFrac, was calculated to illustrate the similarity of the microbial community between the four samples ( Figure S2C). Results showed a clear diversity between the three co-digesting groups and the control test. This could be also demonstrated by the hierarchical clustering analysis ( Figure S2D). This indicated that characteristic conditioning, by the external PWs addition, clearly shifted the microbial community composition of raw WAS, with specific bacteria predominating.
Communities distribution among the four groups was further illustrated at phylum, class and genus level. Out of forty identified phyla, nine (having relatively high abundance) were analyzed and represented in Figure S3A. Bacteroidetes, Chloroflexi, Firmicutes and Proteobacteria were dominant in the four groups, with total relative abundance between 69.8% and 72.9%. They were also reported to be the main phyla accumulated in fermenters, taking WAS, WAS-cornstover and WAS-food wastes as feedstocks [23][24][25]. At the class level, the distribution of 12 taxonomic groups is shown in Figure S3B. The relative abundance of Clostridia and Bacilli substantially increased in the co-digestion tests, which explained the increase of Firmicutes, especially in the case of SSR. Clostridia become the dominant class in PWs-conditioning tests with 30.5%, 21.6% and 33.9% in VR, SR and SSR, respectively (while being only 8.6% in Control). They most likely had a significant impact on the studied process, being able to release hydrolases and utilize various organic matter types to produce VFAs [26]. Bacilli also increased from 0.8% in the Control to 3.8%, 2.7% and 8.7% in VR, SR and SSR, respectively. Alphaproteobacteria, closely related with saccharides utilization and propionic acid production [27], increased in VR, but decreased slightly in SR and SSR.
The relative abundances of genera in different samples is shown in Figure S3C and Table S2 (relative abundance > 1%). Gemmobacter, represented by facultative fermentative organisms, maximized in VR (2.7%), thus contributing to the conversion of large molecular carbohydrates into small molecular ones. Acetoanaerobium, one of the anaerobic H 2 -oxidizing acetogenic bacteria that can produce acetate acid from H 2 and CO 2 [28], decreased in SR (2.8%), but increased in VR (4.9%) and SSR (3.1%). Proteiniclasticum is known for its ability to metabolize soya peptone, tryptone, amino acids as carbon and nitrogen sources to acetate acid, propionic acid and iso-butyric acid, without consuming any of the saccharides [29]. Its proportion was less than 1% in the Control, while it increased in PWs-conditioning samples, reaching the highest proportion in VR (7.4%). Levilinea, isolated by Yamada and colleagues from an anaerobic reactor, is an anaerobic bacterium, which can ferment saccharides and amino acids to produce hydrogen, acetate acid and lactic acids [30]. It clearly increased in PWs-conditioning samples and even reached 9.0% in SR (having a relative abundance of only 3.1% in the Control). Some species of Cloacibacillus are known to ferment amino acids into acetate acid, propionic acid and formate, while producing only butyric acid from serine [31]. In the present study, they represented only 0.2% in the Control, but were enriched to 2.5% in VR and SR. Petrimonas, which can make use of various saccharides to produce acetate acid and propionic acid, increased in PWs-conditioning samples and reached 2.0% in SSR. However, there was less than 0.1% in the Control.

Assessment of PWs Composition
Biological accessibility of different PWs clearly varied because of differences in the processing feedstocks and technologies. For instance, soybean and wheat beans are commonly used for the production of soy sauce, which leads only to a partial utilization of proteins and easily degradable carbohydrates, like starch. This implies that a large proportion of lipids, lignocellulosic substrates and proteins are still available and present in the SSR [32]. In the case of vinegar and alcohol production, crops such as sorghum, barley and peas are widely used substrates in the Shanxi Province. Due to the brewing processes and microbial consortia used, small proportions of lipids and proteins remain in the VR and SR, but with plenty of crude fibers.
The results obtained by FTIR spectroscopic analysis outlined in Section 3.1 confirmed that the employed TA method could efficiently open up the lignocellulosic structures. This could also be further confirmed by the results of the crystallinity of the raw and pretreated PWs ( Figure 1C). The CIs were 1.00, 1.01 and 0.97 for the raw VR, SR and SSR. The raw SR possessed higher crystallinity (CI, 1.01), i.e., higher tensile strength, larger hardness and relative densities [33][34][35], which revealed the relatively complex aggregation state of celluloses. On the contrary, the SSR had the lowest CI of 0.97. After TA pretreatment, 1.5%-3.0% decrease of crystallinity was obtained, which may increase the tenderness, hygroscopicity and chemical reactivity of the feedstock to some extent.

Comprehensive Process Assessment Over Hydrolysis and Acidification Steps
The hydrolysis step is commonly considered as the rate-limiting step in AD processes [36]. Consequently, the efficiency of anaerobic digestibility should be improved if the hydrolysis of i.e., PWs and WAS can be accelerated. Notably, the hydrolysis of proteins and carbohydrates was closely related with two kinds of hydrolytic enzymes, a-glucosidase and protease [37]. The activities of the two enzymes were all enhanced by co-fermenting with PWs, especially for the SSR test ( Figure 2E). It is common knowledge that "hydrolysis-acidogenesis-methanogenesis" typically occur sequentially. Under proper conditions, the produced VFAs could be easily metabolized to CH 4 [38]. It is possible to get the hydrolysis constant k h by applying the CH 4 yield and fit it to Equation (1). As presented in Figure 2F, all the coefficients of determination (R 2 ) were > 0.85. The k h values obtained during the three co-digestions were higher than in the mono-digestion, and peaked in the SSR test (0.0059 h −1 ). That is, the external addition of PWs could indeed boost the sole WAS digestion, especially in the case of SSR. In other words, the enhancement of the high hydrolysis in SSR test was mainly attributed to the increase of C/N of the feedstock, which led to the release of hydrolytic enzymes by hydrolytic bacteria, which accelerated the conversion of the macromolecular organics (i.e., proteins and polysaccharides) to soluble organics (i.e., amino acid and monosaccharide), and provided the abundant substrates for the subsequent acidification process.
This result was further supported by EEM fluorescence spectra analysis (Figure 3). The peak locations of the four DOM-EEM profiles were similar, but with different values of fluorescence intensities (FI). After the calculation using the PARAFAC algorithm, tryptophan, tyrosine, protein and fulvic-like substrates were identified ( Figure S4). The FI of DOM components were obviously accelerated by co-fermenting with PWs ( Figure S5). Com. 1 and Com. 2 were the dominant components for the SSR tests, with the sum FIs of 53, compared to 23, 21 and 17 for the VR, SR and the control tests. This implied that the dissolving of proteins was accelerated by co-fermentation, which was inconsistent with process measurements ( Figure 2C). Small amounts of Com. 4 were also shown in the DOM, which were introduced by the biological metabolism of dead organics and mostly provided by external PWs addition in the co-fermentation tests.
The three PWs had considerably different effects on the enhancement of WAS acidification, with SSR being the preferable conditioning carbon source, showing a 4.2-fold increase over the sole WAS test. In contrast, the ones achieved in the VR and SR tests had roughly a 2.9-fold increase compared to the control test. Previous studies also observed the diversified conditioning effect on VFAs production, when introducing distinct external carbon sources, with e.g., 12-fold increase over sole sludge fermentation, using perennial ryegrass conditioning, but for instance only 1.4-fold increase with bagasse [39,40]. Guo et al. also observed that straws showed a 3-fold increase in VFAs yield over the sole WAS test, which was 1.2-fold higher than that produced by spent mushroom substrates [41]. That is, the external carbon resource played a vital role in enhancing WAS acidification, with different effects produced by the different lignocellulosic biomass selected. Moreover, SSR seemed to be beneficial to 2-3 carbon atom VFAs accumulation, with 40% and 52% increases over that generated in the SR and VR tests. Previous research also showed the mixture of acetate acid and propionic acid was beneficial for the enrichment of denitrifying bacteria and phosphorus accumulating organisms (PAOs) [42].

315
Based on the fact that anaerobic digestion is strongly influenced by the dominant microbial 316 community, 13 key microbiomes (at genus level) were chosen for CCA analysis (Figure 4).

317
Interestingly, only pH values were positively correlated with the first canonical axis (CCA 1) and 318 only propionic acid and pH correspondingly showed negative interrelations for CCA 2 (Table S3).

319
As expected, pH also showed a negative correlation with the VFAs. Arrows signify the direction

Potential Interrelation Between Process Performance and Key Microflora
Based on the fact that anaerobic digestion is strongly influenced by the dominant microbial community, 13 key microbiomes (at genus level) were chosen for CCA analysis (Figure 4). Interestingly, only pH values were positively correlated with the first canonical axis (CCA 1) and only propionic acid and pH correspondingly showed negative interrelations for CCA 2 (Table S3). As expected, pH also showed a negative correlation with the VFAs. Arrows signify the direction and magnitude of process variables related with functional microbes. The length of arrow implies the linkage strength between the process variables and the functional microbes. Taking into consideration the seven selected performance measurements, acetate acid production was intensively correlated with the functional microbes, followed by Sca, Spr and VFAs (according to the length of the vector). Sca and Spr were associated to VFA production (the intersection angles < 90 • ). As indicated, more fermentative acidogenic bacteria were selectively enriched and thereby led to more VFAs accumulating (especially for the SSR test in this study). This is coherent with the result of the microbial community, showing that Papillibacter (1.8%), Levilinea (7.1%), Aquihabitans (4.2%) and Petrimonas (2.0%) predominated in SSR, and were mainly distributed in the positive direction of axis 1 (Figure 4).  On the basis of the results obtained in this work, the underlying mechanism of VFAs recovery from WAS by conditioning with PWs is proposed ( Figure 5). Four processes, including solubilization, hydrolysis, acidogenesis and acetogenesis, are involved in this system. Certain functional microbiomes, especially the hydrolytic acidogenic bacteria, e.g., Proteiniclasticum, Cloacibacillus and Gemmobacter, were related to the composition of feedstock. In this system, Petrimonas and Levilinea are the main monosaccharide-degrading fermenters, which can utilize various carbohydrates to produce acetate acid and propionic acid. Amino-acid-degrading fermenters, which mainly consists of Proteiniclasticum, Cloacibacillus and Levilinea, can convert amino acid to 2-3 carbon atom VFAs. In this sense, a better understanding and control for the shaping of key functional microbiomes, as a reaction to feedstock conditioning, would give precious insight into the value-added biomolecules recovery from environmental biorefinery of WAS and other biomass.

370
Ex/Em (300/440). DOM was obtained from the co-digesting WAS and three PWs; Figure S5: fluorescence 371 spectral intensity of components in the DOM samples obtained from the co-digesting WAS and three PWs. Table S1: α-diversity parameters, Table S2: genera (relative abundance > 1% in at least one samples), Table S3: 373 the eigenvalues of first two canonical axes and their relationships with each environmental factor.

Conclusions
This study investigated the conditioning effects of three cellulosic PWs on acetogenesis enhancement during WAS co-fermentation. A comprehensive study on the "hydrolysis-acidogenesis-methanogenesis" processes was undertaken to explore the underlying mechanism. SSR was more beneficial to VFAs accumulation from co-digestion with WAS, especially 2-3 carbon atom VFAs accumulation. Hydrolysis of particulate organics was significantly prompted by PWs conditioning, as confirmed by the hydrolysis constant and EEM spectra analyses. Results from the fiber analyzer and infrared spectra also revealed that the lower crystallinity of SSR could enhance the biological accessibility for the subsequent AD process. The adjustment effect of feedstock composition on selective enrichment of functional microbial consortia was also noteworthy.

Supplementary Materials:
The following are available online at http://www.mdpi.com/1996-1073/13/14/3630/s1, Figure S1: rarefaction curves of bacterial communities obtained from the co-digesting WAS and three PWs based on sequencing of 16S rRNA gene, Figure S2: overlap of the four bacterial communities based on OTU (3% distance) (A). The shared OTUs were analyzed at the phylum level (B), principal component analysis (C) and hierarchical cluster analysis (D) of bacterial communities obtained from the co-digesting WAS and three PWs based on sequencing of 16S rRNA gene, Figure S3: taxonomic classification of pyrosequences from bacterial communities of five samples at the phylum (A), class (B) and genus (C) levels, Figure S4: four components of DOM decomposed by the PARAFAC approach: (Com.1) tryptophan-like substances, Ex/Em (270/350); (Com.2) tyrosine-like substances, Ex/Em (270/300); (Com.3) protein-like substances, Ex/Em (225/300); (Com.4) fulvic-like substrates, Ex/Em (300/440). DOM was obtained from the co-digesting WAS and three PWs; Figure S5: fluorescence spectral intensity of components in the DOM samples obtained from the co-digesting WAS and three PWs. Table S1: α-diversity parameters, Table S2: genera (relative abundance > 1% in at least one samples), Table S3: the eigenvalues of first two canonical axes and their relationships with each environmental factor.