基于RNA-Seq的羊種布魯氏菌新轉錄本與非編碼RNA鑒定

郭英飛,王玉飛,龔春麗,楊明娟，袁久云,莊妤冰,柯躍華,杜昕穎,汪舟佳,陳澤良

基于RNA-Seq的羊種布魯氏菌新轉錄本與非編碼RNA鑒定

郭英飛1,3,王玉飛2,龔春麗1,楊明娟1，袁久云1,莊妤冰1,柯躍華1,杜昕穎1,汪舟佳1,陳澤良1

目的 對羊種布魯氏菌的轉錄本進行測序，鑒定基因組中新的轉錄本和非編碼RNA。方法 提取羊布魯氏菌的總RNA，去除rRNA后連接接頭逆轉錄成cDNA，PCR擴增后進行測序，以羊布魯氏菌16M的基因組序列為參照對測序得到的reads進行比對作圖，通過生物信息學方法進行新轉錄本和非編碼RNA的鑒定。結果 測序數據分析顯示，reads在基因組上覆蓋度較好。與現有注釋基因相比，有773個基因的5′或3′端在基因組的原有位置基礎上發生了延伸，并發現了16個新的轉錄本。根據測序結果，共鑒定出241條候選的非編碼RNA(sRNA)，進一步的RT-PCR驗證結果顯示，預測的sRNA在體外條件下有表達。結論 布魯氏菌基因組中存在除了預測以外的新的轉錄本，布魯氏菌中存在非編碼RNA且在不同條件下有差異表達。

布魯氏菌；RNA-Seq；轉錄組；sRNA

布魯氏菌病(Brucellosis, 簡稱布病)是由布魯氏菌引起的一種人獸共患傳染病，危害嚴重，在世界各地都有廣泛流行，該病不僅危害人類健康，還影響畜牧業的發展。目前布魯氏菌可分為9種，即羊種、牛種、豬種、犬種、綿羊附睪種、沙林鼠種、田鼠種、鯨型海洋種和鰭型海洋種[1-2]。在我國流行的主要是羊、牛和豬種布魯氏菌，其中羊種布魯氏菌流行最為廣泛。布魯氏菌是一種胞內寄生菌，胞內生存和復制是布魯氏菌的主要毒力特征。布魯氏菌主要是通過改變吞噬體的成熟過程來避免吞噬體與溶酶體的融合，從而使得細菌可以到達其胞內復制部位內質網，開始大量復制[3]。基因組測序與基因的功能注釋，對于深入理解布魯氏菌的基因功能及其在致病性中的作用具有重要意義。

轉錄組是特定細胞或組織在特定時間或功能狀態下轉錄出來的所有RNA的集合，是連接基因組遺傳信息與蛋白質組生物功能之間的紐帶[4-5]。隨著新一代高通量測序技術的發展，轉錄組測序(RNA sequencing，RNA-Seq)技術已成為轉錄組研究的重要手段。RNA-seq能夠對特定細胞或組織在某一環境或生理條件下的幾乎所有轉錄本進行檢測，它在分析轉錄本的表達水平、研究結構變異的同時，還能夠發現未知轉錄本和非編碼RNA[6-8]。與傳統的芯片雜交技術相比，RNA-Seq無需設計探針，也不需要事先知道參考基因組序列的注釋信息，并且具有分辨率高、背景噪音低等優勢，目前已廣泛應用于原核和真核生物的轉錄組學研究。本研究以羊布魯氏菌為例，應用RNA-Seq技術對其進行高通量轉錄組測序分析，描繪出羊布魯氏菌的轉錄組圖譜，鑒定新的轉錄本以及非編碼RNA，為深入理解布魯氏菌與宿主細胞相互作用的機制奠定基礎，同時也為進一步完善布魯氏菌基因結構信息及挖掘潛在的新基因和非編碼RNA提供了有價值的數據。

1 材料與方法

1.1 菌株、培養基及主要試劑

1.1.1 菌株 羊布魯氏菌Brucellamelitensis16M為本室保存。

1.1.2 培養基及培養條件 大豆胰蛋白胨瓊脂(TSA)和大豆胰蛋白胨肉湯(TSB)培養基購自生物梅里埃公司。挑取TSA平板上的羊布魯氏菌單菌落在5 mL TSB過夜培養至后，37 ℃培養至穩定期，然后以1∶50轉接到100 mL的TSB培養基中，培養至對數生長中期。

1.1.3 分子生物學試劑 MasterPureTMRNA Purification Kit購自Epicenter公司，DNase購于美國 Qiagen公司；SYBR Green I定量PCR試劑盒購自TaKaRa公司；Trizol 購自Invitrogen公司；MICROBExpress bacterial mRNA enrichment kit購自美國Ambion 公司；RNA連接酶、CIP 酶、PNK酶購于NEB公司；接頭和引物在美國IDT公司合成；切膠回收試劑盒購于美國Qiagen公司；測序試劑購于美國Illumina公司。

1.2 RNA抽提和rRNA去除 MasterPureTMRNA Purification Kit提取羊布魯氏菌的總RNA，并用DNase消化，方法參照試劑盒說明進行。之后，用MICROBExpress bacterial mRNA enrichment kit從總RNA中去除rRNA，純化的mRNA通過Agilent 2100生物分析儀進行質量檢測。

1.3 cDNA測序文庫的構建及RNA-seq高通量測序 mRNA打斷成100-500 nt短片段后逆轉錄成cDNA，3′末端加polyA并連接測序接頭，連接產物經PCR擴增得到測序文庫，用Illumina HiSeq 2000測序儀對文庫進行序列測定。

1.4 轉錄組數據分析 測序所得的原始數據稱為原始序列(raw reads)。原始序列去除接頭序列、空讀序列(N的比例大于10%的讀數)以及低質量序列(質量值Q≤20的堿基數占整個讀數的40%以上)后得到測序序列(clean reads)用于后續分析。參考基因組16M(AE008917.1和AE008918.1)及它們的注釋基因信息均從NCBI數據庫下載。使用短reads比對軟件SOAPaligner/SOAP2軟件將過濾后的clean reads比對到布魯氏菌基因組和相關基因，通過比對結果統計reads在參考基因組及基因序列上的分布情況及覆蓋度，同時利用RPKM(Reads Per Kb per Million reads)法[10]計算基因的表達量。

基因結構優化及新轉錄本的預測如圖1所示，首先利用cufflinks軟件將比對上羊布魯氏菌基因組的測序序列進行組裝拼接。通過比較gene model與現有注釋基因的差別，對基因的5′和3′端進行延伸，由此優化基因結構。如果組裝的轉錄本序列未能與現有基因比對上，而是位于現有基因之間的基因組上，同時滿足下列條件：距離現有的注釋基因200 bp以上；長度不短于150 bp；平均覆蓋度大于2, 則這些序列有可能為潛在的新轉錄本及新基因。

圖1 基因結構優化及新轉錄本鑒定流程

非編碼RNA(small non-coding RNA，sRNA)的預測主要是參考文獻9，從潛在的gene model中挑選出長度大于等于100 bp且平均覆蓋深度不小于的2的gene model，再從中找出位于基因間區的潛在gene model作為新轉錄本。如果其表達量與兩翼編碼區域相差15%以上，而且向上游或下游延伸200 nt有啟動子的存在，則可作為候選的sRNA。

1.5 sRNA的RT-PCR驗證 挑取候選的sRNA分子進行RT-PCR驗證。將培養至對數生長中期的16M菌液離心，PBS漂洗1遍。將菌體重懸到等量的不同培養基中進行刺激處理后，用Trizol試劑提取布魯菌的RNA。酸刺激是將菌體重懸到pH4.0的TSB培養基中，37 ℃作用20 min；高氧刺激是將菌體重懸到含有440 mmol/L H2O2的TSB培養基(pH7.0)中，37 ℃作用20 min；營養缺乏刺激是將菌體重懸到GEM培養基[10]中(pH7.0)，37 ℃作用20 min；同時將菌體重懸于pH7.0的TSB培養基中，37 ℃作用20 min作為未處理對照。RNA反轉錄成cDNA后用針對候選sRNA分子和16sRNA的引物進行RT-PCR驗證。

2 結 果

2.1 轉錄組的測序數據評估 利用瓊脂糖凝膠電泳和Agilent 2100對提取的羊布魯氏菌總RNA質量進行評估，檢測結果顯示5S、16S和23S RNA條帶完整，23S/16S >1.5，OD260/280≥1.8，OD260/230≥1.8。結果表明RNA純度較高，質量滿足下一步建庫要求。此外，我們通過檢測reads在基因組上的分布來評價測序文庫的隨機性。統計參考基因上不同位置比對上的Reads數量，因為不同的參考基因有不同的長度，我們把 Reads在基因上的位置標準化到相對位置(Reads在基因上的位置與基因長度的比值)。結果表明reads在基因組上分布均勻(圖2A)，說明片段的隨機性好。基因覆蓋度統計表明，覆蓋率達90%～100%的基因數有2 779條，約占88%；覆蓋率達80%～90%的基因數有253條，約占8%(圖2B)。此結果表明，我們的測序結果隨機性好，并且對基因的覆蓋度良好，數據可以用于后續的轉錄組分析。

2.2 基因結構的優化 基因編碼區的起始與終止位置，是通過基因組注釋得到的，但真正的起始與終止位置有待驗證。通過轉錄組測序，能夠對基因的結構進行進一步的確認與優化。基因結構優化分析顯示，羊布魯氏菌基因組上共有773個基因的5′和或3′端在原有基礎上發生了延伸，不同延伸類型的統計結果見表2。從表中可以看出，I號染色體和II號染色體上的延伸分布有差異。部分基因的5′端和3′端的延伸情況見表3。

表1 樣品和參考基因組及參考基因比對的統計結果

表2 基因起始與終止位置延伸的基因分布

2.3 新轉錄本的鑒定 通過RNA-seq，本研究共鑒定出16個新轉錄本，長度分布為166～748 bp，其中13個位于I號染色體上，3個位于II號染色體上，具體信息見表4。

2.4 sRNA的預測及鑒定 通過轉錄組測序，在羊布魯氏菌16M的兩條染色體上共篩選出241條候選的sRNA分子，大小在100～834 nt之間(表5)，其中168個位于I號染色體上，73個位于II號染色體上。我們從中隨機抽取了10個sRNA候選分子用RT-PCR進行了驗證，結果表明其中7個sRNA分子在羊布魯氏菌中存在轉錄本，而且大多數sRNA在不同的環境壓力條件下表達量不同(圖3)，提示這些sRNA有可能在布魯氏菌適應環境壓力中發揮一定的作用。

圖2 Reads隨機性分析(A)及基因覆蓋度統計(B)

表3 部分基因5′端和3′端的延伸

表4 鑒定的新轉錄本信息

表5 部分鑒定的新sRNA信息

T4: TSB 4.0, H: H2O2，G：GEM 7.0，T7: TSB 7.0.

3 討 論

布魯氏菌是一種胞內寄生菌，適應胞內環境是布魯氏菌致病與生存的關鍵，因此了解布魯氏菌在胞內極端苛刻環境下的生存機制對理解其致病機制至關重要。目前，已完成布魯氏菌7個種18株菌的全基因組測序工作，為進一步在分子水平研究布魯氏菌與宿主相互作用奠定了基礎。基于Solexa、454和SOLID平臺的RNA-Seq作為近年來新發展起來的高通量轉錄組測序技術，為轉錄組學研究提供了一種方便、有效的手段。目前該技術已廣泛用于基因表達水平比較分析，同時，該技術在進一步完善基因組結構信息、鑒定新轉錄本及非編碼RNA方面的應用也受到大家的關注[11-14]。本研究應用RNA-Seq技術對羊布魯氏菌進行高通量轉錄組測序分析，為更全面地理解其發病機制提供必需的基礎資料。

對于RNA-Seq來說，沒有基因組DNA污染的高質量RNA是保證其后續測序結果質量的關鍵環節。由于mRNA僅占到總RNA的3%～5%，因此當總RNA提取完成后需要進行mRNA富集。將富集后的mRNA打斷成小片段，并對末端進行接頭修飾，為防止RNA自連，我們對RNA的5′ 末端進行去磷酸化，然后再連接接頭逆轉錄成cDNA，PCR擴增后進行測序。在構建cDNA文庫時，采用特定的處理使其在后續的信息分析中能夠區分轉錄方向，并且確定轉錄本的邊界，從而可以更為準確的計算基因的表達量，同時挖掘sRNA及其它非翻譯區的信息，獲得基因表達譜芯片技術無法獲得的信息。經Illumina高通量深度測序后, 我們得到了一個包含6 666 668條原始測序讀數的測序文庫。比對分析顯示，能比對上參考基因組的reads所占比例為90.28%，其中比對到唯一位置的reads所占比例為51.45%。單位置比對(unique mapped)的reads主要為mRNA，而多位置比對(multi-position mapped)的reads基本上是rRNA和tRNA，因此我們在后續分析中使用的都是unique mapped reads，以保證分析結果的可信度。羊布魯氏菌16M兩條染色體上共有3 199條基因，3 170條基因被覆蓋比對上，其中高峰度表達的基因(比對到基因中的reads數多于10的基因)有3 157條。此結果表明，正常羊布魯氏菌至少表達轉錄了3 170條基因，而且絕大多數是高峰表達的。

RNA-Seq在進一步完善基因組結構信息和鑒定新轉錄本方面也發揮著重要作用。將RNA-Seq測序結果與羊布魯氏菌16M現有基因組上基因進行比對后，發現共有773個基因的5′或3′端在原有基礎上發生了延伸。該結果表明，原基因的5′或3′UTR區的預測存在一定的偏差，而這些延伸進一步優化了相關基因的結構。我們還對沒有與現有基因比對上，而是位于現有基因之間的基因組上的轉錄本序列進行了分析，共發現了16個新的轉錄本，長度為166～748 bp，后續的深入研究需要對這些新轉錄本進行進一步的實驗驗證。

轉錄組測序的另一個重要作用就是挖掘和發現新的非編碼小RNA。細菌基因組中除了編碼基因外，還有一些不編碼基因也可被轉錄出來。由于這類RNA長度較短，具有調控功能，也被稱為非編碼小RNA(small non-coding RNA，sRNA)。細菌sRNA的長度在50～500堿基之間，不編碼蛋白質，常位于基因間區，也有少數位于反義鏈上[15]。目前已在多種細菌中發現了sRNA的存在，其中研究最廣泛的是大腸桿菌和沙門氏菌。sRNA在細菌的生物學功能中發揮了重要的調節作用。目前認為大部分的細菌sRNA在應對環境變化的基因表達調控中發揮重要作用，是細菌適應環境壓力的重要調控子。它們的主要功能是感應環境的變化，調控細胞的代謝途徑、生長方式使之與環境相適應，以及控制細菌毒力基因的表達[16]。我們通過RNA-Seq在布魯氏菌16M的兩條染色體上共發現了241條候選的sRNA分子，大小在100～834 nt之間。這說明sRNA在布魯氏菌中是普遍存在的。布魯氏菌是一種胞內寄生菌，它們主要寄生于機體的單核細胞(主要是巨噬細胞)內。巨噬細胞內微環境包含有多種殺(抑)菌物質，比如酸、過氧化物、溶菌酶等。布魯氏菌必須適應一系列不同的環境狀態才能在宿主體內長期生存，而這種適應主要是通過調控基因的協調表達來實現的[17]。我們隨機挑取了10個sRNA候選分子用RT-PCR進行了驗證，結果表明其中7個sRNA分子在羊布魯氏菌中均存在轉錄本，表明通過RNA-Seq的方法來發現sRNA是非常有效的。此外，我們還分析了這7個sRNA在不同的環境壓力條件下的表達，RT-PCR結果表明除了ncRNACandidate_33外，其它sRNA在不同的刺激條件下表達量不同，說明這些sRNA有可能在布魯氏菌適應胞內環境中發揮一定的作用。

綜上所述，我們將RNA-Seq技術成功地用于布魯氏菌的轉錄組研究中，該結果為描繪布魯氏菌的轉錄組圖譜奠定了基礎。通過本研究，還優化了布魯氏菌的基因結構、發現了新的轉錄本和新的sRNA，這些結果不僅增加了對布魯氏菌現有基因組信息的注釋，也為更全面地理解其致病機制奠定了基礎。

[1]Chain PS, Comerci DJ, Tolmasky ME, et al. Whole-genome analyses of speciation events in pathogenic Brucellae[J]. Infect Immun, 2005, 73(12): 8353-8361. DOI: 10.1128/ IAI.73.12.8353-8361.2005

[2]Corbel MJ. Brucellosis: an overview[J]. Emerg Infect Dis, 1997, 3(2): 213-221. DOI: 10.3201/eid0302.970219

[3]Gorvel JP, Moreno E.Brucellaintracellular life: from invasion to intracellular replication[J]. Vet Microbiol, 2002, 90(1-4): 281-297. DOI: 10.1016/S0378-1135(02) 00214-6

[4]Hou ZW, Wang Z, Gao H, et al. The principle of dRNA-seq and its applications in prokaryotic transcriptome analyses[J]. Hereditas, 2013, 35(8): 983-991. (in Chinese) 侯志偉, 王贊, 高宏,等. dRNA-seq原理及其在原核生物轉錄組學研究中的應用[J]. 遺傳, 2013, 35(8): 983-991.

[5]Gustincich S, Sandelin A, Plessy C, et al. The complexity of the mammalian transcriptome[J]. J Physiol, 2006, 575(Pt 2): 321-332. DOI: 10.1113/ jphysiol.2006.115568

[6]Wang Z, Gerstein M, Snyder M. RNA-Seq: a revolutionary tool for transcriptomics[J]. Nat Rev Genet, 2009, 10(1): 57-63. DOI: 10.1038/nrg2484

[7]Tariq MA, Kim HJ, Jejelowo O, et al. Whole-transcriptome RNAseq analysis from minute amount of total RNA[J]. Nucleic Acids Res, 2011, 39(18): e120. DOI: 10.1093/nar/gkr547

[8]Grabherr MG, Haas BJ, Yassour M, et al. Full-length transcriptome assembly from RNA-Seq data without a reference genome[J]. Nat Biotechnol, 2011, 29(7): 644-652. DOI: 10.1038/nbt.1883

[9]Yoder-Himes DR, Chain PS, Zhu Y, et al. Mapping theBurkholderiacenocepacianiche response via high-throughput sequencing[J]. Proc Natl Acad Sci U S A, 2009, 106(10): 3976-3981. DOI: 10.1073/pnas.0813403106

[10]Wang YF, Qiao F, Zhong ZJ, et al. Transcriptional analysis ofBrucellavirulence regulation genes under stress conditions and during cell infection[J]. Chin J Microbiol Immunol, 2008, 28(10): 919-924. (in Chinese) 王玉飛， 喬鳳， 鐘志軍. 布魯菌重要毒力調控基因在環境刺激和細胞侵染過程中的轉錄研究[J]. 中華微生物學和免疫學雜志, 2008, 28(10): 919-924.

[11]Leyn SA, Kazanov MD, Sernova NV, et al. Genomic reconstruction of the transcriptional regulatory network inBacillussubtilis[J]. J Bacteriol, 2013, 195(11): 2463-2473. DOI: 10.1128/JB.00140-13

[12]Kroger C, Dillon SC, Cameron AD, et al. The transcriptional landscape and small RNAs ofSalmonellaentericaserovarTyphimurium[J]. Proc Natl Acad Sci U S A, 2012, 109(20): E1277-1286. DOI: 10.1073/pnas.1201061109

[13]Mutz KO, Heilkenbrinker A, Lonne M, et al. Transcriptome analysis using next-generation sequencing[J]. Curr Opin Biotechnol, 2013, 24(1): 22-30. DOI: 10.1016/j.copbio.2012.09.004

[14]Reddy JS, Kumar R, Watt JM, et al. Transcriptome profile of a bovine respiratory disease pathogen: Mannheimia haemolytica PHL213[J]. BMC Bioinformatics, 2012, Suppl 15: S4. DOI: 10.1186/1471-2105-13-S15-S4

[15]Gottesman S. Micros for microbes: non-coding regulatory RNAs in bacteria[J]. Trends Genetics, 2005, 21: 399-404. DOI: 10.1016/j.tig.2005.05.008

[16]Toledo-Arana A, Repoila F, Cossart P. Small noncoding RNAs controlling pathogenesis[J]. Curr Opinion Microbiol, 2007, 10(2): 182-188. DOI: 10.1016/ j.mib.2007.03.004

[17]Celli J. Surviving inside a macrophage: the many ways ofBrucella[J]. Res Microbiol, 2006, 157(2): 93-98. DOI: 10.1016/j.resmic.2005.10.002

Identification of novel transcripts and sRNA ofBrucellamelitensisby RNA-Seq

GUO Ying-fei1,3,WANG Yu-fei2,GONG Chun-li1,YANG Ming-juan1,YUAN Jiu-yun1,ZHUANG Yu-bing1, KE Yue-hua1,DU Xin-ying1,WANG Zhou-jia1,CHEN Ze-liang1

(InstituteofDiseaseControlandPrevention,AcademyofMilitaryMedicalSciences,Beijing100071,China)

To identify novel transcripts and sRNA in genome ofB.melitensisby transcriptome sequencing, total RNA were extracted fromB.melitensisculture and rRNA were removed. After the addition of adaptor, RNA was reversely transcribed into cDNA, which were then subjected to PCR amplification and sequencing. The generated reads were mapped to genome sequence ofB.melitensisstrain 16M. With the mapping results, novel transcripts and sRNA were identified by bioinformatics methods. Sequencing results analysis showed that genome sequence was covered with the reads with good quality. A total of 773 genes were extended in their 5′ and/or 3′ ends of their original locations. Sixteen novel transcripts and 241 sRNAs candidates were identified. RT-PCR showed that some of the sRNAs were differentially expressed under stress conditions. InB.melitensisgenome, there is novel transcript which is not predicted. The sRNA does exist inB.melitensisand were expressed under different conditions.

Brucellamelitensis; RNA-Seq; transcriptome; sRNA

s: Wang Yu-fei, Email: yufeiwang21@yahoo.com; Chen Ze-liang, Email: zeliangchen@yahoo.com

10.3969/cjz.j.issn.1002-2694.2015.03.006

國家自然科學基金(81171530)，北京市科技新星資助項目(Z131102000413062)，北京市自然基金項目(6122030)

王玉飛，Email：yufeiwang21@yahoo.com 陳澤良，Email：zeliangchen@yahoo.com

1.軍事醫學科學院疾病預防控制所，北京 100071; 2.武警總醫院檢驗科，北京 100069; 3.空軍司令部門診部，北京 100843

R378.5

A

1002-2694(2015)03-0216-06

Supported by the National Natural Science Foundation of China (No. 81171530), the Beijing Novo Program (No. Z131102000413062), and the Natural Science Foundation of Beijing (No. 6122030)

2014-08-11；

2014-12-05