宏轉錄組測序揭示褐土脲酶基因的表達豐度和細菌來源

胡利偉，軒貝貝，戴華鑫，郭建華，牟文君，劉 陽，翟 振，閆 鼎，蔡憲杰，奚家勤，薛超群，宋紀真*

1. 中國煙草總公司鄭州煙草研究院，鄭州高新技術產(chǎn)業(yè)開發(fā)區(qū)楓楊街2 號 450001

2. 上海煙草集團有限責任公司采購中心，上海市楊浦區(qū)長陽路717 號 200082

脲酶是氮素循環(huán)過程中的一種關鍵的含鎳寡聚酶，土壤脲酶直接參與土壤中含氮有機化合物的轉化，其活性可表征土壤中的氮素狀況［1-2］。一般認為，土壤脲酶的酶活性與土壤中有機物的含量和土壤微生物的數(shù)量等正相關［3］，多地實驗表明施加有機肥比無機肥的煙田土壤脲酶酶活更高，施用有機肥使土壤脲酶活性穩(wěn)定并維持較高水平［4-6］。土壤脲酶的活性隨著煙草的生育期變化，不同時期變化速率不同，煙田土壤的脲酶活性從團棵期到旺長期緩慢增強，從旺長期到現(xiàn)蕾期增速加快，但從現(xiàn)蕾期到采烤前期快速下降，之后于采烤的后期上升至旺長期水平［7-8］。脲酶活性的變化反映了植煙土壤氮素的轉化、積累、供應效應與土壤保肥性之間的協(xié)調(diào)性［9］。煙田土壤脲酶活性與煙株發(fā)病程度也有一定的相關性，例如隨著黑脛病發(fā)病程度加重，煙株根際脲酶活性降低［10］。

土壤脲酶來自土壤微生物,如土壤細菌、真菌等。細菌在土壤中的種類和數(shù)量較多, 并具有較強的產(chǎn)酶能力［11-12］，脲酶活性細菌是與土壤脲酶的性質和數(shù)量有直接關系的一類微生物［13-14］。近年來已經(jīng)從土壤中分離出許多產(chǎn)脲酶菌。張知曉等［15］分離的土壤脲酶活性細菌分屬于厚壁菌門（Firmicutes）和變形菌門（Proteobacteria）2 大類群，包括厚壁菌門的桿菌屬（Bacillus）、芽孢八疊球菌屬（Sporosarcina）和芽孢桿菌屬（Lysinibacillus）以及變形菌門的產(chǎn)堿桿菌屬（Alcaligenes）和假單胞菌屬（Pseudomonas）等，認為土壤脲酶活性細菌優(yōu)勢類群是厚壁菌門，同時也存在變形菌門菌株。新西蘭學者的研究表明枸櫞酸桿菌等9 種細菌以及瓜枝孢菌等24 種真菌參與了牧場土壤中尿素的代謝［16］。也有學者從苗圃根際土壤分離到芽孢桿菌（Bacillus）和芽孢八疊球菌（Sporosarcina）等產(chǎn)脲酶菌株［17-18］。這些微生物編碼的脲酶基因包括結構基因 ureA、ureB、ureC 和輔助基因 ureD、ureH、ureE、ureF 和 ureG 等，其中結構基因 ureA、ureB、ureC 分別編碼γ、β、α亞基，組成尿素酶原，而輔助基因編碼蛋白協(xié)助將Ni2+傳遞至脲酶活性中心，激活酶原［19］。

褐土在河南省境內(nèi)主要分布于沙河與伏牛山一線之北，京廣線以西的偏西北區(qū)域［20］，自然狀態(tài)下褐土的有機質含量多為4%～8%，耕種土壤在2%以下。之前的宏基因組測序結果顯示河南地區(qū)褐土中產(chǎn)脲酶微生物主要為鏈霉菌、節(jié)桿菌和不動桿菌等細菌，足馬杜拉分枝菌等真菌以及亞硝化球菌等古細菌［21］。這些微生物編碼的脲酶基因種類繁多，但是其中哪些脲酶基因表達能最終發(fā)揮作用尚不清楚。 宏轉錄組測序技術（Metatranscriptomic sequencing）通過對特定時期、特定環(huán)境樣品中所有微生物群落的轉錄本進行大規(guī)模高通量測序，可以直接獲得環(huán)境中可培養(yǎng)和不可培養(yǎng)的微生物基因表達信息，目前已經(jīng)越來越多地在環(huán)境微生物學，醫(yī)學微生物等領域應用［22-24］。為此，對來自河南省平頂山、洛陽、許昌、鄭州和安陽等市的5 個高脲酶活性土壤樣本富集培養(yǎng)，進行土壤微生物RNA 提取、建庫及轉錄組測序，分析脲酶基因在轉錄水平的表達豐度，并對微生物來源進行解析，旨在明確褐土尿素代謝中微生物的作用。

1 材料與方法

1.1 土壤樣品的采集

采樣時間為2018 年6～10 月，采樣地點為鄭州市上街區(qū)和滎陽市、洛陽市洛寧縣、平頂山市郟縣、許昌市襄城縣及安陽市林州市的農(nóng)田。取樣深度為0～10 cm 土壤，采用5 點混合法進行土壤取樣，剔除雜物后裝入塑料自封袋，置于10～15 ℃通風條件良好處保存。每個地級市取樣10份，共取土樣50 份，對其脲酶活性進行檢測［21］，選取各市脲酶活性最高的土壤樣品進行下一步實驗。

1.2 產(chǎn)脲酶菌的富集培養(yǎng)

洛陽、許昌、鄭州、安陽和平頂山脲酶活性最高的土壤樣品分別為 XC-10、ZZ-07、AY-06、LY-07和PDS-01，每個樣品取樣3 次（10 g/份），共15 份，分別加入含有50 mL液體培養(yǎng)基（蛋白胨1.0 g/L、氯化鈉1.2 g/L、磷酸二氫鉀 0.8 g/L、葡萄糖 1.5 g/L、硫酸鎳 0.1 g/L、酚紅0.016 g/L、尿素20 g/L，pH 7.0）［25］的250 mL 三角瓶中，在30 ℃搖床中180 r/min振蕩培養(yǎng)24 h，接著將15 瓶富集培養(yǎng)產(chǎn)物進行混合，用液氮冰凍后放入-70 ℃冰箱保存。

1.3 土壤RNA 的提取

使用Powersoil 土壤總RNA 提取試劑盒（美國Mobio 公司）進行土壤RNA 的提取，操作按照試劑盒說明書進行。采用瓊脂糖凝膠電泳檢測RNA的完整度并使用紫外分光光度計（美國Nanodrop公司）檢測RNA 的純度和濃度，將提取的RNA 置于-70 ℃冰箱中保存。

1.4 RNA 的逆轉錄與建庫

利用DNA 探針與RNaseH 酶促反應特點，采用Ribo-off rRNA 去除試劑盒（南京諾唯贊生物科技有限公司）對土壤微生物總RNA 進行處理，以獲得能夠滿足后續(xù)文庫構建的mRNA 及其他非編碼RNA 產(chǎn)物，接著采用VAHTS?Stranded mRNA-seq Library Prep Kit for Illumina?試劑盒（南京諾唯贊生物科技有限公司）對產(chǎn)物進行逆轉錄、雙鏈cDNA 合成、cDNA 片段化修飾、磁珠純化及片段化分選、文庫擴增等處理。文庫用8%的PAGE 膠電泳檢測，利用Qubit 2.0 DNA 檢測試劑盒對回收的DNA 進行精準定量。

1.5 測序及生物信息學分析

使用Hiseq 測序儀（美國Illumina 公司）進行雙端測序，雙向各測150 bp。對原始數(shù)據(jù)質量值等信息進行統(tǒng)計，并使用FastQC 軟件對樣本的測序數(shù)據(jù)質量進行可視化評估。使用Trimmomatic 軟件［26］去除帶 N 堿基的序列、reads 中的接頭序列和低質堿基，隨機從Clean 數(shù)據(jù)中抽取10 000 條序列與 NCBI NT 數(shù)據(jù)庫進行 BLASTN 比對，取e-value ≤1e -10 并且相似度＞90%、覆蓋度＞80%的比對結果，計算其物種分布，進行污染檢測。使用Trinity軟件［27］將clean數(shù)據(jù)de novo組裝成轉錄本，參數(shù)min_kmer_cov 2，其余默認。使用TransDecoder軟件預測Unigene 的CDS 并翻譯成蛋白序列。將獲得的基因和蛋白序列與包括NR、COG、eggNOG等在內(nèi)的相關數(shù)據(jù)庫進行比對，獲得基因的功能及物種注釋信息［28-30］。將拼接得到的轉錄本作為參考序列，使用Bowtie2 將質控后的測序序列與參考序列進行比對，并通過RSeQC 軟件統(tǒng)計比對結果［31］。計算每個特異基因的TPM（每百萬序列轉錄本數(shù)，Transcripts per million）值以分析其表達豐度。

2 實驗結果

2.1 褐土產(chǎn)脲酶微生物的宏轉錄組測序

對各地區(qū)脲酶活性最高的土壤樣品XC10、ZZ07、AY06、LY07 和PDS01 進行富集培養(yǎng)，提取土壤RNA 后進行宏轉錄組測序和分析，初步獲得總數(shù)據(jù)量14 187 012 298 bp，拼接后獲得數(shù)據(jù)量271 Mb 及1 106 689 個轉錄本，預測特異基因940 402個，N50 為 262 bp，最大基因長度 18 629 bp，平均長度288.18 bp（表1）。

表1 轉錄組測序的相關結果Tab.1 Results of metatranscriptomic sequencing

2.2 轉錄本的微生物來源

富集后轉錄本與NCBI 數(shù)據(jù)庫進行比對，結果表明序列匹配數(shù)超過10 條的有4 567 種菌株，序列匹配數(shù)在10～99 條之間的菌株種數(shù)最多，為3 636種，序列匹配數(shù)在200～299 條之間的有154 種，序列匹配數(shù)在300～399 條之間的有92 種，序列匹配數(shù)量隨著菌株種數(shù)的增多基本呈現(xiàn)減少的趨勢，最多的匹配數(shù)為11 408 條序列（圖1），其中豐度最高的10 個菌株為Acidobacteria bacterium DSM 100886、 Brevibacillus choshinensis、 Agromyces sp.NDB4Y10、 Candidatus Rokubacteria bacterium CSP1-6、Bacillus fordii、Gemmatirosa kalamazoonesis、Comamonas testosteroni、Arthrobacter sp. Soil736、Cupriavidus necator 和Sporosarcina koreensis 等（圖 2）。

圖1 序列匹配的菌株數(shù)量分布Fig.1 Quantity distribution of bacteria with matched sequences

圖2 土壤RNA 轉錄本的微生物來源Fig.2 Microbe sources of soil RNA transcripts

2.3 脲酶基因的篩選與注釋

通過注釋篩選3 個及以上拷貝的轉錄本，共獲得122 個脲酶結構基因與115 個脲酶輔助基因。其中脲酶結構基因包含93 個ureC（α亞基）、16 個ureB（β亞基）及13 個ureA（γ亞基）；脲酶輔助基因中有 29 個 ureD、12 個 ureE、36 個 ureF、29 個 ureG、5個ureH 和4 個ureJ。利用脲酶結構基因與輔助基因編碼的蛋白序列對NCBI 的非冗余蛋白庫（Non-redundant database，NR）進行 BLASTP 分析，結果（圖3）表明，122 個脲酶結構蛋白和NCBI 的NR 庫收錄序列的同源性為78%～100%，其中脲酶結構蛋白α、β和γ亞基與NR 庫收錄序列的同源性在90%～100%之間的分別有75、5 和12 個。115個脲酶輔助蛋白和NR 庫收錄序列的同源性為64.25%～99.95%，UreD、UreE、UreF、UreG、UreH及UreJ 和NR 庫收錄序列的同源性在90%～100%之間的分別有7、5、9、20、1 和3 個。而UreD、UreF和NR 庫收錄序列相似度低于80%的蛋白較多，分別有 11 個和 13 個。

2.4 脲酶基因的表達情況

圖3 脲酶亞基和NR 數(shù)據(jù)庫比對的同源性分布Fig.3 Distribution of homology of urease subunits and those in NR database

圖4 脲酶基因的表達情況Fig.4 Expression levels of urease genes

如圖4 所示，ureC、ureB、ureA、ureD、ureE、ureF、ureG、ureH 和ureJ 等脲酶基因中，編碼α亞基的ureC基因數(shù)量最多，表達值TPM 在0～1 之間的有66個，1～2 之間的有 21 個，大于 2 的有 6 個，最高的DN897996_c0_g1 為 3.16。編碼β亞基的 ureB 基因有 16 個，表達值 TPM 在 0～1 之間的有 12 個，1～2之間的有4 個。編碼γ亞基的 ureA 基因有13 個，表達值 TPM 在 0～1 之間的有 9 個，1～2 之間的有 4個。另外ureD、ureE 和ureF 的表達值最高的分別為 5.31、2.00 和 2.71，其中 ureD 表達值 TPM 超過 2的有 3 個轉錄本，分別為 DN911511_c0_g1（2.03）、DN893873_c0_g1（3.31）和 DN907166_c3_g2（5.31）。ureH 和ureJ 數(shù)量較少，表達值也較低。

2.5 脲酶基因的微生物來源

脲酶基因的屬來源如表2 所示。芽孢桿菌屬Bacillus、芽孢八疊球菌屬Sporosarcina、節(jié)桿菌屬Arthrobacter、鏈霉屬Streptomyces 和類芽孢桿菌屬Paenibacillus 等屬細菌具有較多的脲酶結構基因和輔助基因在轉錄水平表達，芽孢桿菌屬Bacillus 中有39 個脲酶基因在轉錄水平表達，其中結構基因ureC、ureB 和 ureA 各有 8、3 和 1 個表達，ureD、ureE、ureF、ureG 和ureH 等脲酶輔助基因分別有10、3、6、6 和 2 個表達在mRNA 水平被檢測到。芽孢八疊球菌屬Sporosarcina 中有24 個脲酶結構基因和輔助基因在 RNA 水平表達，包括 ureC、ureB、ureD、ureE、ureF 和 ureG 等，其中脲酶輔助基因 ureC 和ureF 數(shù)量最多，分別為9 個和7 個。節(jié)桿菌屬Arthrobacter 中 ureC、ureB、ureD、ureF 和 ureG 等脲酶基因，分別有 5、2、4、6 和 3 個基因在轉錄水平表達。 鏈霉屬Streptomyces細菌在RNA水平表達的有ureC、ureB、ureA、ureD、ureF 和ureG 等基因，其中ureC較多，有10 個基因。此外，還發(fā)現(xiàn)了Paenibacillus、Candidatus、 Bradyrhizobium、 Pseudarthrobacter、Ralstonia 及Cupriavidus 等屬細菌的脲酶結構基因和輔助基因轉錄的mRNA。

表2 細菌脲酶基因的屬來源Tab.2 Genus sources of urease genes in bacteria （個）

3 討論

張知曉等［15］研究發(fā)現(xiàn)從云南昆明市農(nóng)田土壤分離出的9 種產(chǎn)脲酶菌分別來自厚壁菌門的桿菌屬、芽孢八疊菌屬、芽孢桿菌屬以及變形菌門的產(chǎn)堿桿菌屬、假單胞菌屬等。沈瓊雯［32］從貴州喀斯特地區(qū)的160 份土壤中成功鑒定了343 株產(chǎn)脲酶細菌，有25 屬79 種，大部分為革蘭氏陽性菌，其中優(yōu)勢菌屬是芽孢桿菌屬、腸產(chǎn)氣桿菌屬、假單孢菌屬和類芽孢桿菌屬，分別占比15.19%、12.66%、12.66%和5.06%。之前的研究采用宏基因組測序方法對河南省代表性褐土中的脲酶微生物來源進行了分析，發(fā)現(xiàn)褐土脲酶中主要細菌來源有芽孢桿菌屬、節(jié)桿菌屬、鏈霉菌屬、芽孢八疊球菌屬、不動桿菌屬、類芽孢桿菌屬和假單胞菌屬等，產(chǎn)脲酶真菌與足馬杜拉分枝菌、藍莓枝枯病菌和嗜熱毛殼菌等親緣關系較近，古細菌中亞硝化球菌、氨氧化古菌和奇古菌等可能也參與了褐土中尿素的降解［21］。本研究通過 BLAST 比對分析，在 mRNA 水平表達的脲酶結構基因和輔助基因與芽孢桿菌屬Bacillus、芽孢八疊球菌屬Sporosarcina、節(jié)桿菌屬Arthrobacter、鏈霉屬 Streptomyces、類芽孢桿菌屬Paenibacillus、念珠菌屬Candidatus、類芽孢桿菌屬Bacillales、慢生根瘤菌Bradyrhizobium、假節(jié)桿菌屬Pseudarthrobacter、勞爾氏菌 Ralstonia、貪銅菌屬Cupriavidus、叢毛單胞菌屬Comamonas、賴氨酸芽孢桿菌屬Lysinibacillus、細桿菌屬Microbacterium 和中華根瘤菌屬Sinorhizobium 等細菌的親緣關系最近。未檢測到編碼脲酶的真菌和古細菌的存在，可能是因為細菌在褐土氮素供應狀態(tài)轉化中作用更 為 關 鍵 。 Bacillus、Sporosarcina、Arthrobacter、Streptomyces 和 Paenibacillus 等 5 個屬 mRNA 水平表達的脲酶結構基因和輔助基因總共108 個，占所有表達脲酶基因的45.6%，這些屬的細菌在褐土脲酶活性的發(fā)揮中起著重要的作用，下一步將對它們的特異抑制劑進行重點研究。

細菌脲酶是由2～3 個不同亞基組成的雜聚肽，其中α亞基在其中起到至關重要的作用。褐土的宏轉錄組測序分析顯示編碼脲酶結構蛋白α亞基的ureC 基因數(shù)量最多，為93 個，這些基因在RNA 水平也維持較高的表達豐度，說明ureC 在褐土中尿素的代謝中有著較為重要的作用，其高表達也與褐土的高脲酶活性一致。UreF 和UreG 形成一個復合物，使酶原對鎳離子處于感受態(tài)，促進鎳離子有效地結合到活性位點。有研究表明ureF基因的敲除能夠降低脲酶活性［33］，這里褐土中ureF 基因表達豐度較高，其豐度最高的轉錄本TPM 值達到了5.37，ureF 的高表達可能是褐土脲酶活性較高的原因之一。UreD 在脲酶金屬中心組裝過程中是鎳插入位點的促進劑也是招募其他輔助蛋白的支點，有研究表明產(chǎn)氣克雷伯氏菌的UreD 是可溶的，并結合2 個當量的鎳離子，UreD不存在時無法激活脲酶［2］。褐土微生物中ureD 基因的數(shù)量和表達值都相對較高，可能與其在脲酶活性發(fā)揮中的重要角色有關。而其他基因ureB、ureA、ureE、ureH 及 ureJ 數(shù)量相對較少，說明其在脲酶復合體的組成中占比可能較低，在脲酶酶活性中的作用相對較小。

4 結論

通過褐土微生物RNA 的宏轉錄組測序，共篩選到轉錄水平表達的122 個脲酶結構基因（93 個ureC、16 個 ureB 和 13 個 ureA）以及 115 個脲酶輔助基因（29 個 ureD、12 個 ureE、36 個ureF、29 個 ureG、5 個 ureH 和 4 個 ureJ）。ureB、ureA、ureG、ureH 和ureJ 等基因在RNA 水平的表達相對較低，TPM 值都在 2 以下，ureC、ureD、ureE 和 ureF 等基因的 TPM值絕大多數(shù)都在2 以下，個別高于2。細菌來源的芽孢桿菌屬Bacillus，芽孢八疊球菌屬Sporosarcina、節(jié)桿菌屬Arthrobacter、鏈霉屬Streptomyces 和類芽孢桿菌屬Paenibacillus 等5 個屬在褐土脲酶酶活的發(fā)揮中貢獻度較大。