污水處理廠排放廢水對受納河段沉積微生物群落結(jié)構(gòu)的影響
——以太平河為例

田雨露,張柏桓,鄢仕偉,孫昊田,宋進(jìn)喜,郭家驊

(1.陜西省河流濕地生態(tài)與環(huán)境重點實驗室,陜西 渭南 714099;2.西北大學(xué) 城市與環(huán)境學(xué)院/陜西省地表系統(tǒng)與環(huán)境承載力重點實驗室,陜西 西安 710127)

隨著我國城市化進(jìn)程的加快與人口的飛速增長,污水處理廠的日處理量與負(fù)荷也逐漸增大[1]。這導(dǎo)致污水處理廠出水被直接排入水環(huán)境中,從而增加了河流水環(huán)境容量負(fù)荷[2]。盡管污水處理廠會將污水經(jīng)過一定處理后達(dá)標(biāo)排放,但是污水處理廠的廢水排放依然會對受納水生態(tài)系統(tǒng)的結(jié)構(gòu)和功能產(chǎn)生影響。

沉積物是城市河流系統(tǒng)的重要組成成分, 是流域陸源污染物的“匯”, 也是水體各種污染物的“源”和“匯”[3]。 相較于浮游微生物, 沉積微生物受水體流動的影響較小, 其群落的結(jié)構(gòu)也相對穩(wěn)定, 能夠較好地反映出環(huán)境因子對沉積微生物種群的長期影響[4]。 近些年,隨著高通量測序技術(shù)的不斷發(fā)展, 人們可以快捷地獲取環(huán)境樣品中的16S rRNA等基因序列, 這些基因序列信息通過與數(shù)據(jù)庫中的已知信息進(jìn)行比對, 可揭示不同環(huán)境介質(zhì)中微生物群落的特點[5]。 Xie等基于DNA元條形碼技術(shù)研究原生生物群落演變特征, 可為生態(tài)壓力源的監(jiān)測和評估提供有效支撐[6]; Yang等基于生態(tài)基因組技術(shù)對沉積物中有毒壓力源進(jìn)行生態(tài)風(fēng)險評估, 發(fā)現(xiàn)此技術(shù)可高效監(jiān)測環(huán)境微生物群落的組成變化[7]; Guo等基于高通量測序技術(shù)研究藻類生長, 得到了綠藻群落結(jié)構(gòu)對抗生素的響應(yīng)機(jī)制[8]; Xie等使用eDNA技術(shù)研究石油污染下沉積物中微生物群落結(jié)構(gòu)的變化, 發(fā)現(xiàn)此技術(shù)可有效地評估人為因素對環(huán)境的影響[9]; 楊琴等使用16S rRNA 技術(shù)揭示了陜北地區(qū)石油污染和未污染土壤中細(xì)菌群落結(jié)構(gòu)的差異[10]。

目前,關(guān)于污水處理廠廢水排放對河流沉積物中微生物群落影響的相關(guān)研究較少,楊思航等探究了不同氣候類型對污水處理廠排放廢水的影響,發(fā)現(xiàn)污水處理廠排放廢水中微生物群落對溫度、pH和電導(dǎo)率的響應(yīng)最為顯著[11];趙琛研究了受納河流微生物群落對污水處理廠排放廢水的響應(yīng),發(fā)現(xiàn)水體中的微生物群落存在空間性差異[12];Wang等對比了不同污染源對太湖沉積物中微生物群落的影響,發(fā)現(xiàn)污水處理廠排放污染對其影響最為顯著[13]。由此可見,當(dāng)前研究多集中于揭示微生物群落多樣性,而缺少微生物之間相互關(guān)系的網(wǎng)絡(luò)推斷分析和菌群的功能潛能預(yù)測分析。因此,在特定的時空變化和環(huán)境過程的驅(qū)動下,微生物之間的相關(guān)性尚不清楚,無法準(zhǔn)確獲取微生物群落在河流生態(tài)系統(tǒng)中的生態(tài)功能。本研究利用16Sr RNA測序技術(shù),基于應(yīng)用生物信息學(xué)分析平臺(QIIME 2),并結(jié)合多種統(tǒng)計學(xué)分析,研究太平河在西安市第六污水處理廠和第七污水處理廠(以下簡稱“六污”和“七污”)排水的影響下,河流沉積微生物群落的組成多樣性、微生物群落的關(guān)聯(lián)網(wǎng)絡(luò)和功能潛能預(yù)測。在此基礎(chǔ)上,揭示顯著影響河流沉積微生物群落的環(huán)境因子。本研究結(jié)合微生物群落、河流水質(zhì)的理化性質(zhì)與環(huán)境因子,可更為全面地探究污水處理廠排放廢水對河流造成的復(fù)合效應(yīng),同時可為評價西安市其他污水處理廠對受納河流水生態(tài)系統(tǒng)的影響提供科學(xué)依據(jù)和數(shù)據(jù)支撐。

1 材料與方法

1.1 研究區(qū)概況和采樣點布設(shè)

西安市位于關(guān)中平原中部,北瀕渭河,南依秦嶺,屬于暖溫帶半濕潤大陸性季風(fēng)氣候,春季風(fēng)多干燥,夏季炎熱多雨,秋季天高氣爽,冬季寒冷干燥,四季分明但春秋較短,降水分布不均。西安市內(nèi)地形平坦,河流眾多,自古以來素有“八水繞長安”之稱。其東西走向的河流有涇河和渭河,南北走向的河流有澇河、灃河、皂河、浐河、灞河等。作為西安市內(nèi)最大的過境河流,渭河發(fā)源于甘肅省渭源縣鳥鼠山,于陜西潼關(guān)匯入黃河,是黃河最大支流[11,14]。其支流皂河發(fā)源于長安區(qū)杜曲街道,河流長度為35 km,依次流過長安區(qū)、雁塔區(qū)、高新區(qū)、未央?yún)^(qū),于渭河城市運動公園處匯入渭河,是一條主要接納污水的排污渠。

太平河是皂河的主要支流,全長約26 km,始于西安市第七污水處理廠,流經(jīng)長安區(qū)、高新區(qū)、雁塔區(qū)、蓮湖區(qū)、未央?yún)^(qū),于未央?yún)^(qū)八興灘村附近匯入皂河[15]。太平河是一條人工河渠,主要接納河流沿岸污水處理廠排水、工業(yè)排水、雨水,具有接納污水,防洪蓄水等功能。本研究以太平河在七污和六污排水口附近的河段為研究對象,根據(jù)《水質(zhì)采樣方案設(shè)計技術(shù)規(guī)定(HJ 495—2009)》進(jìn)行采樣點的布設(shè)。在七污排水口處(S1)及中下游(S2、S3為中游,S4、S5為下游)共布設(shè)5個采樣點,在六污排水口上中下游(S6、S7為上游,S8、S9為中游,S10、S11為下游)共布設(shè)6個采樣點,S1到S11的位置分布如圖1所示。

圖1 太平河采樣點布設(shè)圖

1.2 河流沉積物和水樣的采集與保存

沉積物樣品:本次研究于2021年1月13日在全部11個采樣點采集沉積物樣品,并使用GPS記錄了每一個采樣點的經(jīng)緯度坐標(biāo)。在河道中使用采泥斗和不銹鋼鏟來采集表層5 cm的沉積物樣品,在1 m2的區(qū)域隨機(jī)抽取3個子樣本進(jìn)行混合,作為一個樣品,將混合之后的沉積物裝入已清洗干凈并貼有編號標(biāo)簽的25 mL離心管中。運回實驗室之后,于-80 ℃低溫保存。24 h內(nèi)將樣品低溫(干冰)運送至派森諾生物公司進(jìn)行高通量測序和微生物群落的分析。

水樣:在沉積物采樣點處,同時采集距沉積物表面5～10 cm水樣,將其裝入已清洗干凈并貼有編號標(biāo)簽的500 mL采樣瓶中。運回實驗室之后,于-20 ℃低溫保存。1周內(nèi)完成水樣化學(xué)分析。

1.3 水樣的化學(xué)分析

在采樣現(xiàn)場,使用川步tp101工業(yè)溫度計測量水樣溫度(t,單位:℃);使用YD-1H筆式鹽度計測量水樣的鹽度(salinity,單位:10-6);使用JBP-607A溶解氧分析儀來測量水樣的溶解氧含量(DO,單位:mg/L);使用pH-100型pH計來測量水樣的pH值;使用ORPscan10筆式ORP計測量水樣的氧化還原電位(ORP,單位:mV);使用力辰科技筆式電導(dǎo)率儀/TDS筆測量水樣的電導(dǎo)率(conductivity,單位:S/m)和溶解性總固體(TDS,單位:10-6)。

1.4 沉積物樣品中微生物測序

本文采用Illuminm測序[16]方法,對沉積物樣品進(jìn)行基因組DNA抽提,抽提完成后,利用0.8%瓊脂糖凝膠電泳檢測抽提的基因組DNA,利用紫外分光光度計對DNA進(jìn)行定量分析。PCR采用NEB Q5 DNA高保真聚合酶,選用細(xì)菌16S rDNA V3-V4區(qū)特異性引物,正向引物為F:ACTCCTACGGGCAGCA,反向引物為R:TCGGACTACHVGGGTWTCTAAT。PCR反應(yīng)過程見表1。利用Illumina公司的TruSeq Nano DNA LT Library Prep Kit進(jìn)行文庫構(gòu)建。對合格的文庫,在4 MiSeq機(jī)器上利用MiSeq Reagent Kit V3(600cycles)進(jìn)行2×250 bp的雙端測序。

表1 PCR 擴(kuò)增過程

1.5 生物信息學(xué)分析

1.5.1 物種分類學(xué)注釋

對于原始測序數(shù)據(jù),首先進(jìn)行初步篩查,對問題樣本進(jìn)行重測、補測,隨后將初步篩查完成后的樣本進(jìn)行文庫和樣本劃分,并去除barcode序列,按照QIIME2 dada2分析流程對序列進(jìn)行去噪處理,獲得ASVs。采用QIIME2的classify-sklearn算法對每個ASV特征序列進(jìn)行物種分類學(xué)水平注釋;采用QIIME2對ASVs特征序列進(jìn)行抽平處理,便于后續(xù)分析。

1.5.2 分類學(xué)組成分析

采用QIIME2和自編pre腳本統(tǒng)計各樣本的分類單元數(shù),并進(jìn)行分類學(xué)組成分析,采用R語言(ggraph、ggplot、ggtree和phyloseq包)對抽平后的ASVs數(shù)據(jù)進(jìn)行統(tǒng)計,可以實現(xiàn)各樣本在門、綱、目、科、屬、種6個分類水平上組成分布的可視化。本研究基于門水平進(jìn)行分類學(xué)組成分析,并以柱狀圖的形式展現(xiàn)結(jié)果。

1.5.3 Alpha多樣性分析

采用QIIME2和R語言(ggplot包)計算Alpha多樣性指數(shù),使用R語言腳本繪制箱線圖,直觀地呈現(xiàn)出不同樣本間Alpha多樣性指數(shù)的差異。本研究以7種不同的Alpha指數(shù)來表征微生物群落不同的指標(biāo)。多樣性指數(shù)標(biāo)簽下的數(shù)字為Kruskal-Wallis 檢驗的p值,當(dāng)p≤0.05 時,代表其具有較高的統(tǒng)計學(xué)差異可信度[17]。

使用R語言的vegan包,繪制稀疏曲線與物種累積曲線。通過繪制稀疏曲線的方式來探究Alpha多樣性指數(shù)與抽平深度的關(guān)系,繪制物種累積曲線以探究樣本數(shù)量是否能夠有效反映群落豐富度。

1.5.4 Beta多樣性分析

采用QIIME2和R語言(ape、vegan、ggtree包)計算距離矩陣,進(jìn)行主座標(biāo)(PcoA)分析,并將結(jié)果繪制成散點圖,以表征Beta多樣性。

1.5.5 物種差異與標(biāo)志物種分析

采用Python LEfSe包進(jìn)行LEfSe (LDA effect size)分析,通過分類學(xué)分支圖的形式,對所有分類水平同時進(jìn)行差異分析,以快速找出引起組內(nèi)差異的物種。使用基于多組學(xué)數(shù)據(jù)的機(jī)器學(xué)習(xí)方法,進(jìn)行隨機(jī)森林分析,通過熱圖與柱狀圖組合的形式呈現(xiàn)引起組建差異的標(biāo)志物種。

以上生物信息學(xué)分析均在派森諾基因云平臺上進(jìn)行(www.genescloud.cn)。

1.6 統(tǒng)計學(xué)分析

1.6.1 關(guān)聯(lián)網(wǎng)絡(luò)分析

將相關(guān)閾值(r)和卡方檢驗值(p)構(gòu)建相關(guān)性矩陣,使用隨機(jī)矩陣?yán)碚摯_定閾值,再采用igraph[18]構(gòu)建關(guān)聯(lián)網(wǎng)絡(luò)數(shù)據(jù),并繪制豐度大小前15的物種模塊關(guān)聯(lián)網(wǎng)絡(luò)圖。使用指示關(guān)鍵物種的拓?fù)渲笖?shù)Zi值(within-module connectivity,節(jié)點與其所屬模塊內(nèi)其他節(jié)點的連接情況)和Pi值(among-module connectivity,節(jié)點與其他模塊內(nèi)節(jié)點的連接情況)進(jìn)行分析,將2個值可視化,便于尋找關(guān)鍵物種。

1.6.2 功能潛能預(yù)測

使用PICRUSt2軟件[19]對微生物群落的代謝通路(pathways)進(jìn)行比較分析。通過柱狀圖的形式展現(xiàn)太平河上下游的微生物群落功能差異。再使用分層的樣本代謝通路豐度表進(jìn)行通路的物種組成分析,探究承擔(dān)差異功能的菌種。

1.6.3 RDA冗余分析

RDA 能夠有效地對多個指標(biāo)進(jìn)行統(tǒng)計學(xué)檢驗, 并確定對水質(zhì)環(huán)境變化具有最大解釋能力的最小變量[20]。本研究使用R語言軟件包vegan對 ASVs 的豐度數(shù)據(jù)(門水平)和水樣的水質(zhì)指標(biāo)數(shù)據(jù)進(jìn)行線性回歸整合(見表2),以三序圖的形式呈現(xiàn)微生物群落與環(huán)境因子之間的關(guān)系。

表2 樣本測序量統(tǒng)計表

2 結(jié)果與討論

2.1 微生物的測序結(jié)果

本次測序一共得到2 403 313對序列,去除低質(zhì)量序列之后得到2 161 929對序列,去噪處理后得到有效序列2 072 120對,經(jīng)過拼接和去除嵌合體后得到1 443 218對高質(zhì)量序列,再將序列中樣本豐度為1的單元去除,得到1 399 712對序列(見表2)。本次測序共得到23個門、51個綱、113個目、181個科、276個屬和86個種,經(jīng)過處理得到的ASVs中,物種分類學(xué)注釋能夠注釋到的域、門、綱、目、科、屬、種的占比分別為1%、0%、4%、3%、27%、59%和6%(見圖2)。

圖2 物種分類學(xué)注釋

2.2 微生物群落多樣性分析

2.2.1 微生物群落的物種組成

圖3為基于門水平的微生物群落的物種組成,由圖可知,太平河中優(yōu)勢菌種為變形菌門(proteobacteria)。變形菌種類繁多,在污水處理廠的活性污泥中有較高豐度[21],可以作為河流中的優(yōu)勢菌種。在S1點發(fā)現(xiàn)藍(lán)細(xì)菌門(cyanobacteria),且豐度很高,這是由于S1點距離七污排污口很近,廢水排出導(dǎo)致此處溶解氧和營養(yǎng)物質(zhì)含量較高(溶解氧8.8 mg/L,總磷1.92 mg/L,均為所有樣本中最高含量,見表3),適合藍(lán)細(xì)菌的生長[22]。

表3 水樣的水質(zhì)指標(biāo)數(shù)據(jù)

圖3 基于門水平的分類學(xué)組成分析

2.2.2 微生物群落的Alpha多樣性

圖4 Alpha多樣性指數(shù)檢驗

圖5為物種稀疏曲線,由圖可知,每一個樣本的曲線都在抽平深度30 000～40 000 bp趨于平緩。圖6為物種累積曲線,由圖可知,隨著樣本數(shù)量的增加,曲線逐漸趨于平緩。由此可見,本次測序深度和樣本數(shù)量都足夠反映微生物群落的豐富度變化情況。

圖5 稀疏曲線

圖6 物種累積曲線

2.2.3 微生物群落的Beta多樣性

圖7為PCoA分析結(jié)果,由圖可知,在置信度較高的PCo1維度上(置信度35.3%,大于PCo2維度的14.8%),受七污影響的河段與受六污影響的河段微生物群落之間有較大的差異,以七污排污口和中游與六污中下游的差異最為明顯。

圖7 PCoA分析結(jié)果

2.3 物種差異與標(biāo)志物種分析

圖8為LEfSe分析結(jié)果,可以看出,受七污影響引起差異的物種為藍(lán)細(xì)菌門、α-變形菌門(alphaproteobacteria)和變形菌門;受六污影響的河段引起差異的物種為酸桿菌門(acidobacteria)、厚壁菌門(firmicutes)和放線菌門(actinobacteria)。

圖8 LEfSe分析結(jié)果

圖9為隨機(jī)森林分析結(jié)果,可以看出,引起樣本之間差異最顯著的物種為藍(lán)細(xì)菌門和厚壁菌門。藍(lán)細(xì)菌門在七污口處尤為顯著,是因為七污口處溶解氧和營養(yǎng)元素含量較高;而厚壁菌門在六污影響河段中最為顯著,六污以改良型AAO工藝為主,陳重軍等研究發(fā)現(xiàn)厚壁菌常存在于厭氧氨氧化反應(yīng)器的污泥中[25],厚壁菌隨著六污工藝的廢水排放進(jìn)入太平河中。

圖9 隨機(jī)森林分析結(jié)果

2.4 關(guān)聯(lián)網(wǎng)絡(luò)分析

圖10為關(guān)聯(lián)網(wǎng)絡(luò)繪制圖,由圖可知,藍(lán)細(xì)菌門與大部分的變形菌門呈正相關(guān)關(guān)系,與小部分呈負(fù)相關(guān)關(guān)系;擬桿菌門與藍(lán)細(xì)菌門呈較強的負(fù)相關(guān)關(guān)系;衣原體門與絕大部分變形菌門呈負(fù)相關(guān)關(guān)系;小部分變形菌門模塊與其他變形菌門模塊呈負(fù)相關(guān)關(guān)系。圖11將Zi值和Pi值可視化,得到太平河微生物群落中的關(guān)鍵物種為變形菌門中的α-變形菌綱。

圖10 關(guān)聯(lián)網(wǎng)絡(luò)繪制圖

圖11 綱水平Zi-Pi圖

2.5 功能潛能預(yù)測

圖12為太平河上下游代謝通路的差異,通過對比 KEGG 數(shù)據(jù)庫得到太平河上游顯著代謝通路為L-賴氨酸發(fā)酵醋酸鹽和丁酸鹽(P163-PWY);下游顯著代謝通路為dTDP-N-乙酰托莫胺生物合成(PWY-7315)。圖13為代謝通路物種組成分析的結(jié)果,由圖可知,承擔(dān)太平河上游微生物群落功能的主要物種為生氧光細(xì)菌(oxyphotobacteria);承擔(dān)下游微生物群落功能的主要物種為放線菌(actinobacteria)。

圖12 代謝通路差異分析

(a) P163-PWY 的物種組成 (b) PWY-7315 的物種組

2.6 環(huán)境因子與微生物群落的相關(guān)性

圖14(a)為采樣點微生物群落與環(huán)境變量之間的關(guān)系,受七污排水影響的河段與受六污排水影響的河段在軸1上具有顯著差異,其中:受七污排水影響的河段對總磷、溶解氧和總固體等因子的響應(yīng)最顯著;受六污排水影響的河段對氨氮、氧化還原電位、溫度等因子的響應(yīng)較為顯著。從物種角度來看,藍(lán)細(xì)菌門與總磷和溶解氧呈正相關(guān),而與鹽度和溫度呈負(fù)相關(guān);厚壁菌門與溫度、氧化還原電位、氨氮等呈正相關(guān),與硝酸鹽氮、溶解氧等呈負(fù)相關(guān)。表4為不同環(huán)境因子的解釋率,其中溫度(30.2%)、溶解氧(25.3%)和鹽度(9.8%)為解釋率最高的環(huán)境因子,本次采樣時間在冬季,排污口排出的廢水溫度會隨著河流流動而迅速降低,所以溫度因子的影響最為顯著。圖14(b)為它們與微生物群落之間的關(guān)系,由圖可知,總固體(2.5%)、總氮(1.4%)、CODMn(0.8%)等環(huán)境因子對于整個太平河微生物群落的影響不顯著。

(a) 采樣點微生物群落與環(huán)境因子的關(guān)系 (b) 微生物群落與起顯著影響環(huán)境因子的關(guān)系

3 結(jié)論

本次測序共得到1 443 218對高質(zhì)量序列,有多半的ASVs可以注釋到屬以上,證明了本次測序的準(zhǔn)確性。在物種組成上,太平河沉積微生物群落優(yōu)勢種為變形菌,隨著采樣點與排污口距離的增加,微生物群落的豐富度逐漸增加,說明了污水處理廠排放的廢水對河流沉積微生物群落存在消極的影響,引起差異的標(biāo)志物種為藍(lán)細(xì)菌門和厚壁菌門。承擔(dān)群落功能的關(guān)鍵物種為變形菌門中的α-變形菌綱。太平河上游與下游的差異功能是L-賴氨酸發(fā)酵醋酸鹽和丁酸鹽和dTDP-N-乙酰托莫胺生物合成,它們分別是由生氧光細(xì)菌和放線菌來承擔(dān)的。溫度、鹽度和溶解氧是對菌群影響程度最大的環(huán)境因子,其中,藍(lán)細(xì)菌門分別和總磷與鹽度呈正相關(guān)與負(fù)相關(guān),厚壁菌門則分別和溫度與溶解氧呈正相關(guān)與負(fù)相關(guān)。

綜上所述,污水處理廠排放廢水會改變河流沉積微生物群落結(jié)構(gòu),而河流沉積微生物群落結(jié)構(gòu)的變化可以用來反映污水處理廠廢水對河流水質(zhì)的影響,進(jìn)而可以表征受納河流生態(tài)系統(tǒng)的健康情況。