莠去津降解菌泛基因組測序及比較基因組分析

王婭麗 朱姍姍 楊峰山 馬玉堃 付海燕 劉春光

（1. 黑龍江大學農業微生物技術教育部工程研究中心，哈爾濱 150500；2. 黑龍江大學生命科學學院 黑龍江省寒地生態修復與資源利用重點實驗室，哈爾濱 150080；3. 黑龍江大學生命科學學院 黑龍江省普通高校微生物重點實驗室，哈爾濱 150080）

莠去津又名阿特拉津，是一種選擇性內吸傳導型苗前苗后除草劑［1］，主要通過植物根部吸收并向上傳導，抑制雜草的光合作用使其枯死［2］，由于其成本較低，使用方便，價格低廉效果好，已成為世界范圍內廣泛使用的三嗪類除草劑，在美國的年使用量高達36 000 t［3］，主要用于玉米、甘蔗、高粱田和草坪雜草的防除。據報道，2004年莠去津被歐盟正式禁止使用，但美國、印度和中國等農業大國仍在使用［4-6］。莠去津施用后大部分進入土壤，其在環境中結構穩定，半衰期為56-364 d，被微生物代謝分解的過程緩慢，在墾區的一些農場發生莠去津后茬敏感作物減產甚至絕產的藥害情況，嚴重影響了農業種植結構的調整和農業生產的安全［7-10］。因此尋找一種有效解決莠去津殘留污染的方法迫在眉睫。生物修復是目前用來消除環境污染物的有效治理方法，已成為當前國內外環境科學研究的熱點。國內外科研人員相繼篩選出多種降解莠去津的土壤微生物［11］，結果發現細菌、真菌、放線菌和藻類中均有能降解莠去津的菌株［12］。

1995年，Mandelbam［13］從長期施用莠去津的土壤中分離到第一株以莠去津為唯一氮源的革蘭氏陰性假單胞菌Pseudomonassp. ADP［14］，降解基因闡明后，成為莠去津降解微生物的模式菌株。此后，研究人員又發現了一株革蘭氏陽性菌Arthrobacter aurescensTC1［15］，并克隆了其降解基因，然而這兩株菌的降解基因并不完全一致。近期有一些莠去津降解菌的相關報道，無論是利用普通PCR擴增還是高通量測序均發現不同種屬細菌的降解相關基因型也并不完全一致，這可能與不同種屬菌株的降解特性有關［16-18］。

本課題組從黑龍江省長期施用莠去津的玉米田土壤中成功分離出以該除草劑為唯一碳源生長的高效降解菌株Microbacteriumsp. HBT4，該菌108 h對于100 mg/L的莠去津的降解率能夠達到99%以上，為黑土農田莠去津污染土壤的修復研究工作提供了良好微生物資源。據報道，Microbacterium屬的菌株多見于轉化蔗糖、合成低聚果糖、產絮凝劑、產纖維素酶、降解多環芳烴、赤霉烯酮、重金屬等研究［19-23］，然而作為莠去津的降解菌株卻少有報道，且未見該菌比較基因組學分析相關報道，可能與文獻報道的基因型不同，也可能是由于適應環境而產生若干突變，存在基因水平轉移現象，有待進一步研究。為更加深入了解該菌株與已報道的兩株模式菌株基因組的異同，本研究對降解菌HBT4進行泛基因組測序，與GenBank中已登錄的兩株模式菌株進行比較基因組分析，旨在找到3株菌株參與莠去津降解的代謝途徑，為明確菌株HBT4降解莠去津的相關分子機理奠定基礎。

1 材料與方法

1.1 材料

供試菌株：菌株HBT4是從黑龍江省長期施用莠去津的玉米田土壤中分離純化出的降解菌。菌株Arthrobacter aurescensTC1和Pseudomonassp. ADP的基因組數據來源于NCBI數據庫（https：//www.ncbi.nlm.nih.gov/）。

1.2 方法

試驗流程按照Illumina公司提供的標準規程執行，在Illumina Hiseq 4000測序平臺對菌株HBT4進行測序，包括DNA小文庫制備和測序試驗。首先提取樣品DNA，其次使用超聲波隨機打斷，獲得需要的插入片段，然后在片段的粘性末端添加引物構建文庫，將構建好的測序文庫在測序芯片上進行橋式PCR 后進行雙端測序，最后將產生的數據解讀成核酸序列用于下一步信息分析。

對測序產生的原始reads進行評估，再通過過濾低質量、糾錯等方法得到高質量的reads，然后使用Velvet軟件［24］對菌株HBT4的測序數據進行組裝，通過軟件Prodigal［25］對組裝后的基因組進行編碼基因預測，使用Repeat Masker［26］軟件對細菌基因組進行重復序列預測，使用軟件tRNAscan-SE［27］預測基因組中的tRNA，使用軟件Infernal 1.1［28］基于Rfam［29］數據庫預測基因組中的rRNA，通過編碼基因、重復序列和非編碼RNA預測對基因組組分進行分析，使用通用數據庫和專用數據庫對基因功能進行注釋，通過SNP、Small InDel SV檢測水平轉移基因獲得基因間變異，通過核心基因組、物種特有基因和共線性分析與系統進化樹構建進行比較基因組學分析。

通用數據庫注釋方法如下：利用預測得到的基 因 序 列 與 COG［30］、KEGG［31］、Swiss-Prot［32］、TrEMBL［32］、Nr［33］等功能數據庫做 BLAST［34］比對，得到基因功能注釋結果。基于nr數據庫比對結果，應用軟件 Blast2GO［35］進行 GO［36］數據庫的功能注釋。利用軟件HMMER基于Pfam［37］數據庫進行Pfam功能注釋。另外根據注釋結果進行進一步的COG類別分析、KEGG代謝通路分析和GO功能分類分析。專用數據庫注釋方法如下：利用預測得到的基因的蛋白序列與轉運蛋白分類數據庫（TCDB）［38］、病原體-宿主互作因子數據庫（PHI）［39］、抗生素抗性基因數據庫（ARDB）［40］、毒力因子數據庫（VFDB）［41］等功能數據庫做 BLAST［34］比對，得到相應的注釋結果。另外利用軟件HMMER基于碳水化合物相關酶數據庫（CAZyme）［42］進行碳水化合物酶類基因的功能注釋。

應用MUMmer［43］軟件包將參考菌株與測序菌株基因組進行基因組比對，找到SNP和Small InDel；利用MUMmer軟件包將參考菌株與測序菌株基因組進行基因組比對，根據得到的同源共線性區域的信息識別出每一個菌株相比于參考基因組發生結構變異的區域；利用兩種軟件Alien Hunter［44］和HGT-Finder［45］分別對菌株Arthrobacter aurescensTC1和菌株Pseudomonassp. ADP進行水平轉移基因預測；使用軟件Mugsy［46］對將所有測序菌株的基因組序列和參考基因組的基因組序列進行比對，從比對結果中找出各菌株共有的序列為核心基因組序列，剩余序列為非必需基因組序列；使用OrthoMCL［47］對各個菌株預測得到的蛋白序列以及參考基因組的蛋白序列進行基因家族分析，尋找每個菌株共有的和特有的基因家族；使用測序物種及參考物種的單拷貝蛋白序列，利用phyML［48］軟件構建進化樹，研究物種間的進化關系。

2 結果

2.1 測序結果統計

對菌株構建270 bp文庫，結果得到測序數據量為0.92 Gb，測序深度為260X，測序質量值在20以上的堿基比例為96.98%，測序質量值在30以上的堿基比例為92.24%。

2.2 基因組組裝結果

組裝得到的基因組為3.53 Mb，Scaffold片段數量為10，Scaffold N50的長度為0.67 Mb，Contig 片段數量為17，Contig N50的長度為0.62 Mb，GC含量為69.98%。

2.3 基因組組分分析

菌株HBT4預測到的基因個數為3 397，所有基因的長度和為3.26 Mb，平均每個基因的長度為959 bp。菌株HBT4預測出的重復序列總長度為46 967 bp，基因組中重復序列含量為1.33%。非編碼RNA即不編碼蛋白質的RNA，針對非編碼RNA的結構特點，采用不同的策略預測不同的非編碼RNA。預測到的tRNA數量為48，預測到的種類數量為43；預測到的rRNA數量為6，預測到的種類數量為3；預測除tRNA和rRNA之外的其它ncRNA，預測到的rRNA數量為9，預測到的種類數量為8。

2.4 基因功能注釋

2.4.1 通用數據庫注釋 通用數據庫功能注釋統計信息如表1所示。

2.4.2 專用數據庫注釋 注釋結果的統計如表2所示。

2.5 菌株間變異差異分析

結果顯示，菌株Arthrobacter aurescensTC1的SNP堿基個數為3 114，Small InDel堿基個數為132，菌株Pseudomonassp. ADP的SNP堿基個數為624，Small InDel堿基個數為48；菌株Arthrobacter aurescensTC1和菌株Pseudomonassp. ADP的Deletion及Insertion的結構變異數量均為0；菌株Arthrobacter aurescensTC1兩種軟件均預測是水平轉移基因的基因數量為77個，兩種軟件預測到的水平轉移基因的并集為1 509，菌株Pseudomonassp. ADP兩種軟件均預測是水平轉移基因的基因數量為24，兩種軟件預測到的水平轉移基因的并集為1 301。

表1 基因功能注釋統計信息

表2 基因功能注釋統計信息

2.6 莠去津降解菌間比較基因組分析結果

2.6.1 核心基因組與非必需基因組 比對結果統計顯示核心基因組序列占總基因組大小的比例為0.02%，基因組序列長度為2 769 bp；非必需基因組序列占總基因組大小的比例為99.98%，基因組序列長度為12.45 Mb。任意兩個菌株共有而其他基因組沒有的基因組序列占總基因組大小的比例為0.75%，基因組序列長度為93 264 bp；單個菌株特有的基因組序列占總基因組大小的比例為99.25%，基因組序列長度為12.35 Mb。菌株HBT4占總基因組大小的比例為0；Arthrobacter aurescensTC1占總基因組大小的比例為58.44%，基因組序列長度為7.22Mb；Pseudomonassp. ADP占總基因組大小的比例為41.56%，基因組序列長度為5.13 Mb。

核心基因組與非必需基因組比例分布圖如圖1所示（圖中Germs代表菌株HBT4）。多個菌株的泛基因組中各組分的比例分布圖，左側餅圖展現的是泛基因組中核心基因組占泛基因組的比例以及非必需基因組占泛基因組的比例。中間餅圖展現的是非必需基因組中，不同個數菌株共有（但不是所有菌株共有）的非必需基因組以及單個菌株特有的非必需基因組各占非必需基因組的比例。右側矩形展現的是各個菌株分別特有的基因組占物種特有基因組的比例。

2.6.2 核心基因家族和物種特有基因 各菌株共有的基因家族即為核心基因家族，核心基因家族和非必需基因家族的統計結果如圖2所示。Microbacteriumsp. HBT4與Arthrobacter aurescensTC1有1891個相同基因家族，Pseudomonassp. ADP與Microbacteriumsp. HBT4有1115個相同基因家族，Arthrobacter aurescensTC1與Pseudomonassp. ADP有1341個相同基因家族，它們之間共同的基因家族有986個。

圖1 基因組與非必需基因組比例分布圖

圖2 基因家族韋恩圖

物種特有基因——泛基因基因家族隨物種的添加變化情況為菌株HBT4隨著物種由上至下的添加核心基因組中基因家族的數量為1 891，隨著物種由上至下的添加泛基因組中基因家族的數量為5 372；Arthrobacter aurescensTC1隨著物種由上至下的添加核心基因組中基因家族的數量為4 072，隨著物種由上至下的添加泛基因組中基因家族的數量為4 072；Pseudomonassp. ADP隨著物種由上至下的添加核心基因組中基因家族的數量為986；隨著物種由上至下的添加泛基因組中基因家族的數量為8 599。將數據做成直方圖方便比較二者之間的關系如圖3所示。該圖是由上面的統計表格繪制而成的，沿著縱軸逐漸添加菌株，核心基因家族的數量逐漸減少（磚紅色矩形），泛基因家族的數量逐漸增加（藍色矩形）。

圖3 泛基因家族隨物種的添加量變化

對于每個菌株的特有基因（物種特有的基因家族的基因）進行GO數據庫注釋分析，結果如圖4所示。細胞膜上的特有基因在細胞組分GO分類中所占比例最大，具有催化功能的特有基因在分子功能GO分類中所占比例最大，代謝途徑中的特有基因在生物進程GO分類中所占比例最大。

經過泛基因組測序和比較基因組分析，發現共有基因和特有基因主要集中在三羧酸循環KEGG代謝通路中，如圖5所示。三羧酸循環KEGG代謝通路中的共有基因和特有基因功能注釋如表3所示。

2.6.3 系統發育樹構建 系統進化關系結果如圖6所示。能夠看出Microbacteriumsp. HBT4與Arthrobacter aurescensTC1進化程度相近，而Pseudomonassp. ADP進化程度較高。

2.6.4 基因共線性分析Microbacteriumsp .HBT4，Arthrobacter aurescensTC1與Pseudomonassp. ADP共線性結果如圖7所示。圖中左邊是相應的菌株名，右邊的長矩形框組成的每一行表示一個菌的基因組，每個矩形框表示一個scaffold，不同基因組之間彩色的連線表示上下兩個基因組之間同源基因在各個不同細菌基因序列上的不同位置信息得到核酸水平上的共線性關系。從圖中可以看出Microbacteriumsp.HBT4與Arthrobacter aurescensTC1基因共線性情況明顯高于Pseudomonassp. ADP與Microbacteriumsp.HBT4，所以Microbacteriumsp. HBT4與Arthrobacter aurescensTC1親緣關系更近一些。

圖4 特有性基因的GO二級節點注釋分類

圖5 共有和特有基因的KEGG代謝通路富集圖

表3 三羧酸循環KEGG代謝通路中的共有基因和特有基因功能注釋

圖6 各菌株與參考菌株的系統進化樹

3 討論

據報道細菌對多環芳烴的降解途徑最終通過進入三羧酸循環得以完全降解，推測菌株的降解特性與三羧酸循環具有相關性。史延華等［49］初步推測菌株Pseudomonas stutzeriYC-YH1對萘的降解是通過水楊酸途徑，萘首先被其代謝為1，2-二羥基萘，而后被轉化為水楊酸和鄰苯二酚，最后進入三羧酸循環被徹底降解。張丹等［50］綜述了惡臭假單胞菌、分枝桿菌屬、芽胞桿菌屬、微球菌屬及糞產堿菌等細菌對菲的降解是通過水楊酸途徑或者鄰苯二甲酸途徑降解，經過一系列反應進一步氧化開環后，最終進入三羧酸循環徹底降解為二氧化碳和水。另外，高亞平等［51］研究表明，莠去津的長期作用改變了鰻草組織氮含量和碳氮比，降低了糖和三羧酸循環中間產物的水平，這說明長期的莠去津暴露會降低能量供應并改變碳氮代謝。

本研究中3個菌株的基因組經比較獲得的共有和特有基因KEGG代謝通路，在由丙酮酸生成乙酰輔酶A時，菌株HBT4、ADP和TC1首先在丙酮酸脫氫酶作用下生成2 -羥基乙基ThPP，然后在丙酮酸脫氫酶作用下生成乙酰二氫丙醇酰胺，最后在二氫硫辛酰胺S-乙酰轉移酶作用下生成乙酰輔酶A，而不是由丙酮酸在丙酮酸脫氫酶（1.2.7.1）作用下直接脫氫生成乙酰輔酶A，這與經典三羧酸循環稍有差異。在第一次脫氫由異檸檬酸生成α-酮戊二酸時，3株菌都是先在異檸檬酸脫氫酶作用下轉化成草酰琥珀酸鹽，再經過脫氫生成α-酮戊二酸，而不是直接由異檸檬酸在異檸檬酸脫氫酶（1.1.1.41/1.1.1.286）作用下生成α-酮戊二酸。第2次脫氫時，3株菌都是先在氧戊二酸脫氫酶作用下生成羧基羥丙基ThPP，然后脫氫生成琥珀酰二氫硫辛酰胺，最后在二氫硫基賴氨酸殘基琥珀酰轉移酶作用下生成琥珀酰輔酶A，而不是α-酮戊二酸在α-酮戊二酸脫氫酶（1.2.7.3）作用下直接生成琥珀酰輔酶A，其中特有基因（2.3.1.61）就位于α-酮戊二酸與琥珀酰輔酶A的代謝通路之間。

圖7 各菌株與參考菌株的共線性圖

另外，由KEGG通路圖可知，α-酮戊二酸與抗壞血酸鹽和醛酸鹽代謝、丙氨酸、天冬氨酸和谷氨酸鹽代謝、D-谷氨酰胺和D-谷氨酸代謝之間存在雙向的間接影響；纈氨酸、亮氨酸和異亮氨酸的降解對琥珀酰輔酶A代謝產生間接影響；酪氨酸代謝、精氨酸和脯氨酸代謝對延胡索酸產生間接影響；丙氨酸、天冬氨酸和谷氨酸代謝、乙醛酸鹽和二羧酸鹽代謝與草酰乙酸之間存在雙向的間接影響。纈氨酸、亮氨酸和異亮氨酸的降解還有脂肪酸代謝對乙酰輔酶A有間接影響；乙酰輔酶A對脂肪酸合成和線粒體中脂肪酸延長有間接影響；草酰乙酸在磷酸烯醇丙酮酸羧激酶作用下生成的磷酸烯醇丙酮酸鹽對糖酵解/糖異生有間接影響，糖酵解/糖異生對丙酮酸代謝也存在間接影響。

三株菌共有基因所在的代謝途徑有可能是其參與莠去津降解的原因，而三羧酸循環中特有基因編碼的二氫硫基賴氨酸殘基琥珀酰轉移酶可能是Microbacteriumsp. HBT4與Pseudomonassp. ADP和Arthrobacter aurescensTC1的不同之處。三羧酸循環的中間產物，從理論上講可以循環不消耗，但是由于循環中的某些組成成分還可參與合成其他物質，而其他物質也可不斷通過多種途徑而生成中間產物，所以說三羧酸循環組成成分處于不斷更新之中。

4 結論

本研究通過泛基因組測序和比較基因組分析，得到該菌株基因組大小約為3.53 Mb，預測到菌株HBT4編碼基因3 397個、重復序列含量為1.33%、非編碼RNA 63個，通用數據庫基因功能注釋共3 324個，專用數據庫基因功能注釋共1 149個，通過菌株間差異變異分析發現SNP、Small InDel和水平轉移基因，未發現結構變異基因，獲得該菌株特有基因中GO注釋到的基因在細胞組分、分子功能和生物學進程中的數量和比例，從KEGG 代謝通路富集圖中發現特有基因編碼的二氫硫基賴氨酸殘基琥珀酰轉移酶位于三羧酸循環中α-酮戊二酸和琥珀酰輔酶A的代謝通路之間。獲得3個菌株核心基因組與非必需基因組比例分布、系統進化樹和共線性關系，發現三者之間共有基因家族986個、菌株HBT4特有基因家族1 171個。該菌株與兩株模式菌株相比，其基因家族之間既有相同之處，又有較大差異。