張 堃，謝喬光，徐小云，韓曉鴿，郭雄飛，杜 娟，程學慶

（1.廣東輕工職業(yè)技術(shù)學院， 廣東 廣州 510300；2.中國電器科學研究院股份有限公司,廣東 廣州 510300）

地球生態(tài)系統(tǒng)中病毒無處不在，幾乎可以感染所有的細胞生命體，通過裂解宿主細胞、用輔助代謝基因重編程宿主代謝或?qū)λ拗骰蜻M行水平基因轉(zhuǎn)移等來影響宿主功能特性［1-3］。如病毒通過感染裂解宿主細胞，使得胞內(nèi)營養(yǎng)物質(zhì)被釋放到環(huán)境中影響土壤碳循環(huán)［4］；通過影響宿主介導(dǎo)的代謝通路溶解環(huán)境中的有機碳和總氮［5］；有機肥的長期使用增加了細菌與病毒間的相互作用，提高了兩者間功能基因的交換頻率［6］。總之，病毒會對宿主微生物的組成和代謝功能造成壓力并最終影響元素的地球化學循環(huán)。

廢水處理設(shè)施（Wastewater treatment plants，WWTP）采用獨特工程微生物生態(tài)系統(tǒng)的活性淤泥進行城市生活污水或工業(yè)廢水處理，是現(xiàn)代生物技術(shù)的重要應(yīng)用之一［7］?；钚杂倌嘀泻w多種功能性微生物群體，以細菌群體為主進行污水中營養(yǎng)物質(zhì)（氮和磷等）或有機物的降解，而細菌群落結(jié)構(gòu)的組成對系統(tǒng)穩(wěn)定性及淤泥沉降等具有重要的影響［8］。以往研究證明WWTP中含有高濃度的生物量以及營養(yǎng)物質(zhì)，是病毒良好的棲息地，其中病毒數(shù)目可能是自然水生生態(tài)系統(tǒng)中含量的 10 到 1 000 倍［9］。Li等人發(fā)現(xiàn) WWTP 的活性淤泥中病毒的存在比較廣泛，涵蓋一些致病性病毒［10］。Wang等人對2處WWTP的長期研究表明地理隔離的污水處理系統(tǒng)中存在一群核心病毒，對污水處理可能存在未知但穩(wěn)定的影響［11］。生存條件如溫度、pH值和營養(yǎng)狀況等的動態(tài)變化會直接影響環(huán)境中病毒的活性，促進病毒群落的變異［12］。目前為止，對于WWTP連續(xù)運行過程中微生物群落組成和功能的研究很多，但對于其中病毒類群的組成和動態(tài)變化的研究較少。

由于病毒變異是高度動態(tài)的過程［13］，短期研究可能無法全面反映在WWTP的復(fù)雜環(huán)境中病毒類群的動態(tài)變化。本次研究以WWTP中EGSB厭氧反應(yīng)器中的污泥為研究對象，同時為更好地研究病毒的空間分布差異，我們每月從EGSB罐體的上層、中層及下層分別取樣，連續(xù)12個月共計獲得36個樣品，結(jié)合宏基因組高通量測序技術(shù)對樣品內(nèi)病毒的組成和豐度進行研究，主要包括：（1）鑒定WWTP的EGSB體系內(nèi)病毒群落組成和多樣性；（2）挖掘不同月份、罐體不同部位間病毒的動態(tài)變化規(guī)律；（3）鑒定系統(tǒng)內(nèi)病毒-宿主關(guān)聯(lián)性及預(yù)測病毒潛在功能組成。本研究加深了對EGSB反應(yīng)器內(nèi)病毒群落結(jié)構(gòu)及功能變異模式的理解。

1 材料和方法

1.1 污水處理系統(tǒng)及樣品采集

本研究以位于廣東省的一處工業(yè)規(guī)模的生產(chǎn)性膨脹顆粒污泥床（Expanded granular sludge beds，EGSB）反應(yīng)器污水為研究對象。實驗取樣從2020年1月持續(xù)到12月，每個月15號分別從反應(yīng)器的上層、中層和底層采樣管處用無菌針筒吸取污水污泥樣本，并干冰運輸立即轉(zhuǎn)移到實驗室中-40℃保存。

1.2 DNA提取及宏基因組高通量測序

采用QIAGEN強力土壤基因組DNA提取試劑盒（天根生化科技（北京）有限公司，型號2745288）對所有收集到的36個污泥樣品進行總DNA提取，具體提取方法參考試劑盒提供的方法說明書。采用瓊脂糖凝膠電泳檢驗提取到的總DNA的質(zhì)量，采用Qubit3.0對樣品DNA進行準確定量。將符合標準的DNA送至廣東美格基因科技有限公司（廣州，中國）進行Illumina NovaSeq 6000雙端測序。平均每個樣品原始測序數(shù)據(jù)量約 12 Gb，總數(shù)據(jù)量 440 Gb。

1.3 病毒序列分析

利用 Trimmomatic（Version 0.39）軟件對下機數(shù)據(jù)進行質(zhì)控；采用Megahit（Version 1.2.5-beta）軟件進行數(shù)據(jù)組裝；CD-HIT（Version 4.6.7）進行序列去冗余（相似性＞90%，覆蓋度＞80%）。VIBRANT（Version 1.2.1）對篩選保留到的長度5 Kb以上的序列進行病毒鑒定；采用CheckV（Version 0.8.1）對預(yù)測到的病毒信息進行質(zhì)量評估，去除質(zhì)量評級中未明確的序列。采用Prodigal（Version 2.6.3）對鑒定獲得的病毒序列進行ORF預(yù)測，并結(jié)合DRAMv（Version 1.2.4）進行功能注釋。采用VPF-tools對病毒序列進行物種歸類及宿主預(yù)測，其中對病毒物種注釋域水平上隸屬度（Membership ratio）與置信度（Confidence score）閾值設(shè)定0.2以上；科水平二者閾值設(shè)定為0.3；屬水平二者閾值設(shè)定為0.25；而對病毒-不同分類層級宿主（域、門和科）預(yù)測中隸屬度（Membership ratio）與置信度（Confidence score）閾值均設(shè)定 0.2 以上［14］。采用Bowtie 2計算病毒在各樣品中的相對豐度。采用 QIIME2（Version 2020.8）進行病毒群落α和β多樣性分析。

1.4 統(tǒng)計分析

樣品間共存和特有病毒類型花瓣圖采用R語言包“Venn diagram”進行分析繪制；采用R（Version 4.2.1）對不同處理組間病毒群落顯著性差異統(tǒng)計分析，P＜0.05判定為組間存在顯著性差異；從NCBI refseq下載已知的噬菌體terL氨基酸序列，與本次注釋到的terL基于最大似然法（Maximum likelihood）采用 MEGA（Version X）構(gòu)建蛋白樹，并采用iTOL在線網(wǎng)頁進行樹圖美化；采用本地Python編制的腳本（Local script）進行病毒與宿主關(guān)聯(lián)?；鶊D繪制等。

2 結(jié)果與討論

2.1 病毒序列鑒定及豐度動態(tài)變化

36個樣品中共獲得3 360條非冗余病毒序列，僅7.5%的序列為溶原性病毒（Prophage），一般認為溶原性病毒對宿主的代謝以及環(huán)境元素循環(huán)的影響遠低于烈性噬菌體［15］。樣品中Chao1指數(shù)在 1 503.60 到 3 046.80 間；Shannon 指數(shù)在8.20到9.58間，說明持續(xù)運行過程中污泥內(nèi)病毒群落存在較高的動態(tài)變化。其次，對組間α多樣性的顯著性差異分析發(fā)現(xiàn)Chao1指數(shù)在連續(xù)月份間呈現(xiàn)規(guī)律性變動趨勢：多樣性先從一二月份顯著升高至三四月份后急劇降低，并在十月份開始再次發(fā)生顯著升高至十二月；而Shannon指數(shù)則無明顯規(guī)律（見表1）。由于病毒豐度不均勻，導(dǎo)致病毒豐富度最高的三四月份整體多樣性并不最高，可能存在一些優(yōu)勢病毒類型影響樣品中病毒整體多樣性的體現(xiàn)。反應(yīng)器不同部位的病毒群落均不存在顯著性差異（P＜0.05），可能差異明顯的采樣月份是導(dǎo)致反應(yīng)器不同部位不存在顯著差異的主要原因。

表1 組間Chao1和Shannon指數(shù)差異顯著性分析

采用Bray-Curtis PCoA法進行不同分組間病毒類群整體離散性分析（見圖1）。樣品按照采樣月份能各自聚類，表明相應(yīng)月份內(nèi)的病毒整體組成上較為相似；而在采樣部位主導(dǎo)下，分組內(nèi)離散性較高，不存在明顯聚類模式。有研究表明AS體系中微生物群落結(jié)構(gòu)出現(xiàn)明顯季節(jié)性演替，春季和冬季采集的樣品細菌群落更相似［16］。有趣的是，本體系中病毒群落在春季和冬季也出現(xiàn)同樣的規(guī)律，可能源于病毒主要依靠侵染的宿主存活，其更傾向于直接受到病毒群落結(jié)構(gòu)的影響［17］。

圖1 基于Bray-Curtis PCoA的不同分組下樣品散點聚類分析

2.2 病毒物種組成及遺傳多樣性

所有樣品中長尾噬菌體科（Siphoviridae）（35.91%±7.84%）、短尾病毒科（Podoviridae）（13.35%±4.54%） 和肌 病毒 科（Myoviridae）（10.44%±2.54%）相對豐度最高，均屬于噬菌 體 目（Caudovirales）（ 見 表2）。 與 以 往 對WWTP中病毒類群的研究結(jié)果一致［18］。樣品間病 毒 屬 以 T4virus（4.74%±2.41%）、P22virus（2.94% ±1.45%）、Lambdavirus（2.63%±0.68%）及P12024virus（2.21%±0.40%）等為主，其中豐T4virus屬的成員被認為具有裂解性但不會產(chǎn)生毒性蛋白［19］。

表2 病毒科層級物種組成信息統(tǒng)計表

噬菌體目（Caudovirales）在不同科中末端酶大亞基（terminase large subunit，terL）具有一定保守性，利用terL基因的系統(tǒng)發(fā)育分析，可以了解噬菌體的進化地位和遺傳距離［20］。本次提取到227個terL氨基酸序列中存在部分與參考蛋白差異過大，并在系統(tǒng)遺傳樹中分支為一些獨立的簇（見圖2），表明本研究體系中未知噬菌體科存在較高多樣性和新穎性，與近期Shi等人的研究一致［21］。

圖2 基于terL氨基酸序列的噬菌體科水平系統(tǒng)發(fā)育樹。紅色的樹枝以及紅色的圓點代表未能準確定位的噬菌體蛋白。

2.3 病毒-宿主關(guān)聯(lián)預(yù)測

體系中鑒定到 2 785 條噬菌體（Phage或Prophage）和3條噬真菌體（Mycophage）。由于病毒的復(fù)雜侵襲能力，還存在532條病毒序列可以同時侵染細菌、真菌或古菌。Gulino等人對紐約14處污水處理廠的研究發(fā)現(xiàn)活性淤泥中噬菌體的比例高達90%［22］， 與本研究結(jié)果一致。

環(huán)境條件影響噬菌體對宿主的選擇范圍［23］，本體系中噬菌體主要侵襲的宿主細菌包括變形桿菌門（Proteobacteria）、放線桿菌門（Actinobacteria）及厚壁桿菌門（Firmicutes）等；主要宿主屬包括分支桿菌屬（Mycobacterium）、假單胞桿菌屬（Pseudomonas）及芽孢桿菌屬（Bacillus） 等（ 見 圖3）。 其 中， 短 尾 病毒科（Podoviridae）預(yù)測到的宿主大部分屬于變形桿菌門（Proteobacteria）；長尾噬菌體科（Siphoviridae）的宿主則以放線桿菌門（Actinobacteria）的分支桿菌屬（Mycobacterium）為主。噬菌體-宿主感染的模式會對潛在的菌群進化和元素循環(huán)產(chǎn)生影響，而噬菌體-宿主的感染可以跨越多種門和屬的范圍［24］。

圖3 病毒科與宿主屬或門間的潛在關(guān)聯(lián)桑基圖

2.4 病毒促氮代謝功能潛力

廢水可能是環(huán)境中活性氮的最大來源之一，會導(dǎo)致水體和生態(tài)系統(tǒng)富營養(yǎng)化；而WWTP中的工程微生物可促進氮代謝，經(jīng)過包括反硝化、硝化反應(yīng)及氨同化作用等途徑促進活性氮的轉(zhuǎn)變［25］。本次共檢測到7條噬菌體含有氮相關(guān)功能基因（見圖4），其中ctg1051569為沙粒病毒屬（Omegavirus），被檢測含有亞硝酸鹽還原酶基因（nitrite reductase，nirK），可以在反硝化通路中促進亞硝酸鹽還原為氧化一氮［26］。而ctg467926（Cd119virus）含有羥胺還原酶基因（hydroxylamine reductase，hcp），可以將羥胺還原為銨根離子進入?yún)捬醢毖趸緩交虬蓖饔猛緩剑?7］。其余5條病毒序列被預(yù)測到含有氨同化作用過程中的天冬酰胺合成酶（asparagine synthase，ASN），研究發(fā)現(xiàn)該酶可能會提高宿主耐銨性［28］。

圖4 7條含有氮代謝相關(guān)基因的病毒結(jié)構(gòu)基因組序列圖

此外，ctg1051569（Omegavirus）在10月份和12月份時豐度較高，提示可能該月份下反硝化作用細菌豐度可能較高；ctg467926（Cd119virus）在8月份樣品中最高，提示該月份下厭氧氨氧化途徑或氨同化作用途徑參與的宿主可能豐度較高；擁有天冬酰胺合成酶（asnB）的病毒在月份間分布較為分散，提示各個月份下促進氨同化作用的病毒豐度差異較?。徽w上十月和十一月?lián)碛械x相關(guān)AMG功能的病毒豐度最高，說明這兩個月份下可能菌群代謝更為活躍（見圖5）。

圖5 7條關(guān)鍵病毒在12個月份間的豐度聚類熱圖

3 結(jié)論

（1）連續(xù)運行12個月的WWTP活性淤泥內(nèi)病毒群落在月份間存在較高的動態(tài)變化，Chao1指數(shù)出現(xiàn)有規(guī)律的動態(tài)變化，即呈先升高后降低并再次升高的趨勢。

（2）各自月份下的病毒群落相似性較高，但采樣罐部位的影響不顯著；春天和冬天體系下的樣品中病毒群落組成更為接近。

（3）淤泥體系中烈性噬菌體占比最高，且存在一些高度新穎且多樣化的病毒科。

（4）體系中含有氮代謝相關(guān)AMG的噬菌體，且在十及十一月份相對豐度更高。