巨 艷，李子華,王婭娜,趙嘉慶，朱明星，李君良，趙 巍

細粒棘球絳蟲原頭蚴mRNA測序及表達譜分析

巨 艷1,2,3，李子華1,3,王婭娜1,3,趙嘉慶1,3，朱明星1,3，李君良1,3，趙 巍1,3

目的 通過對細粒棘球絳蟲原頭蚴的mRNA的測序及表達譜分析，初步建立起細粒棘球絳蟲原頭蚴的表達譜數(shù)據(jù)庫，了解細粒棘球絳蟲原頭蚴基因表達及蛋白構成情況，為全面了解細粒棘球絳蟲原頭蚴生物學特征及寄生蟲與宿主之間的關系奠定基礎并為新的診斷方法、篩選新的藥物靶點和疫苗候選分子選擇提供理論依據(jù)。方法 用TRIZOL法提取人源細粒棘球絳蟲原頭蚴的總RNA，構建細粒棘球蚴的轉錄組測序文庫，Illumina的solexa測序平臺對RNA進行測序并進行生物信息學分析。結果 測序結果去雜后得到2G數(shù)據(jù)，通過從頭拼接我們得到18 569個contig，這些contig的總長度為71 329 bp，contig平均長度為384 bp，最小的contig長度為201 bp，最大contig長度為4 618 bp，N50(覆蓋50%所有核苷酸的最大序列重疊群的長度)為384 bp。預測得到unigene為9 029條，將這9 029條基因與NCBI的nr數(shù)據(jù)庫做blast比對，最終有7 441條unigene具有同源比對信息。結論 根據(jù)GO分析可以發(fā)現(xiàn)，共有10 550條unigene與數(shù)據(jù)庫中的基因有較高同源性，且較多的unigene可以與多條基因相對應，一共建立了10 550條對應關系。通過與KEGG數(shù)據(jù)庫進行比對分析，細粒棘球絳蟲原頭蚴的轉錄組中有4 731條unigene得到注釋，這4 731個得到注釋的基因位于241條代謝通路中，這些代謝通路分別與代謝過程，基因信息過程，環(huán)境相關過程，細胞過程及與人類疾病相關。

原頭蚴；轉錄組學；生物信息分析；表達譜

細粒棘球蚴病(echinococcosis)也被稱為囊性包蟲病(hydatid disease)是一種人獸共患的寄生蟲病。包蟲病影響范圍廣泛，具有全球性分布趨勢，在各個大陸基本都有包蟲病被發(fā)現(xiàn)，有報道顯示100個國家都存在包蟲病患者[1]。包蟲病嚴重危害人畜健康，該病廣泛流行于世界各地的畜牧業(yè)發(fā)達地區(qū)。包蟲生活史復雜，在細粒棘球絳蟲由六鉤蚴發(fā)育到原頭蚴時會形成包囊，包囊的囊壁由多層的非細胞層和細胞層構成[2]，對原頭蚴形成了良好的保護，這導致宿主的免疫攻擊及藥物都不能對其造成傷害，臨床的觀察也發(fā)現(xiàn)一些細粒棘球絳蟲原頭蚴在進入宿主體內后會被宿主的免疫系統(tǒng)抑制并殺傷，但也有一些棘球蚴能夠逃逸宿主的免疫攻擊而生存下來，我們對這一期間細粒棘球絳蟲的發(fā)育及免疫機制認識有限，這給包蟲病的有效診斷及靶分子的篩選帶來了挑戰(zhàn)。

不論是尋找具有特異性的診斷抗原分子還是篩選出有效的靶蛋白分子，亦或是研發(fā)具有真正療效的新型藥物都需要對細粒棘球絳蟲原頭蚴有更深刻和全面的認識，對細粒棘球絳蟲原頭蚴的發(fā)育機制或免疫機制都需要更深入的研究，而這些研究的前提是對細粒棘球絳蟲原頭蚴的基礎有深入全面的認識，對細粒棘球絳蟲原頭蚴的基因構成，蛋白組成情況及其功能研究是上述一切研究的基礎及理論支持。而這正是包蟲病研究的薄弱環(huán)節(jié)，因此研究者期望從整體出發(fā)全面研究細粒棘球絳蟲原頭蚴的基因構成及蛋白組成情況。

1 材料與方法

1.1 細粒棘球絳蟲原頭蚴 實驗所需原頭蚴來自于包蟲病人的包囊，包囊由外科手術剝離所得，包囊主要由寧夏醫(yī)科大學附屬醫(yī)院提供。無菌條件下抽取包囊中的囊液，分離原頭蚴，再用PBS洗滌3次。

1.2 細粒棘球絳蟲原頭蚴總RNA提取 總RNA提取方法參照Invitrogen公司的 Trizol Reagent 說明書進行。

1.3 細粒棘球絳蟲原頭蚴的轉錄組測序 細粒棘球絳蟲原頭蚴的轉錄組測序委托北京基因組完成

1.4 序列的拼接及注釋 序列分析利用 Illumina 平臺將測序所得的圖像數(shù)據(jù)轉化為相應的序列數(shù)據(jù)，對所產生的原始序列文件進行質量評估和可信度分析，并去除測序過程中低質量的序列和不確定的序列(Q<20)。之后利用Trinity分析軟件對所獲得的高質量序列進行轉錄組 de novo 拼接， 以得到contig序列。先用repeatmasker去參考序列中的簡單重復序列，augustus軟件，Pasa軟件，cufflinks進行基因預測并將上述結果整合。生成的uni轉錄物與NCBI的nr數(shù)據(jù)庫進行比對注釋，然后將生成的uni轉錄物與GO數(shù)據(jù)庫進行比對注釋，并獲得相應的GO功能分類。最后與KEGG數(shù)據(jù)庫進行比對獲得代謝通路數(shù)據(jù)。

2 結 果

2.1 細粒棘球絳蟲原頭蚴mRNA的測序及轉錄組數(shù)據(jù)的組裝結果分析 mRNA測序完成后，我們對結果進行去雜處理共得到2G的clean reads。將clean reads通過從頭拼接我們得到18 569個contig，這些contig的總長度為71 329 bp，contig平均長度為384 bp，最小的contig長度為201 bp，最大contig長度為4 618 bp，N50(覆蓋50%所有核苷酸的最大序列重疊群的長度)為384 bp。對轉錄組contig的長度特征進行分析，其中長度在200～300 bp范圍內的contig數(shù)目有10 265條，占總數(shù)的55%，300～400 bp內的contig有3 488條占總數(shù)19%，≥400 bp的contig共有14 816條，比例為36%(見表1)。

表1 細粒棘球絳蟲原頭蚴轉錄組contig數(shù)據(jù)組裝質量統(tǒng)計

Tab.1 Data assembly for contig digital transcriptome of protoscolex

2.2 細粒棘球絳蟲原頭蚴轉錄物功能注釋及分類 經過augustus、pasa、cufflinks三種軟件對contig進行基因預測，最終預測得到unigene為9 029條，我們將這9 029條基因與NCBI的nr數(shù)據(jù)庫做blast比對，最終有7 441條unigene具有同源比對信息，在這些比對信息中，其中2 984條unigene與中華肝吸蟲(Clonorchissinensis)同源占總數(shù)的40.71%，其次與曼氏血吸蟲和日本血吸蟲的同源性也較高，分別占總數(shù)的27.01%和8.8%。不知道功能的假設蛋白和未命名蛋白分別為1 228條和68條。在這些匹配的同源基因中與細粒棘球絳蟲相吻合的98條基因，還有69條基因與多房棘球絳蟲同源。

根據(jù)GO分析可以發(fā)現(xiàn)，共有10 550條unigene與數(shù)據(jù)庫中的基因有較高同源性，且較多的unigene可以與多條基因相對應，一共建立了10 550條對應關系，從而能夠得到盡可能全面的注釋和分類。細粒棘球絳蟲原頭蚴的轉錄組中的Unigene根據(jù)GO功能大致可分為生物過程、細胞組分和分子功能3大類44分支。由表2可知,在所有轉錄本中，有4 912個轉錄本(50.75%)的GO注釋歸類為生物學過程，2 820個(29.14%)歸為細胞組分，1 947個(20.12%)被歸為分子功能。在細胞組分類型中，細胞和細胞組分所含比例較高，分別占34.82%和30.60%。與突觸相關轉錄本最少僅有3條。在分子功能中細胞進程相關基因最多有967條，起連接作用和有催化活性的轉錄本都較多占總數(shù)的39.29%和48.79%，與核酸連接調控轉錄的和與蛋白相連調控轉錄的unigene最少，分別只有4條和5條。

2.3 細粒棘球絳蟲原頭蚴轉錄組的pathway注釋分析 通過與KEGG數(shù)據(jù)庫進行比對分析，細粒棘球絳蟲原頭蚴的轉錄組中有4 731條unigene得到注釋，這4 731個得到注釋的基因位于242條代謝通路中，這些代謝通路分別與代謝過程，基因信息過程，環(huán)境相關過程，細胞過程及與人類疾病相關。其中代謝map00230即嘌呤代謝涉及的基因個數(shù)最多有112個(圖1)。

3 討 論

同一有機體在不同的時期所表達的基因是不完全相同的，具有特定的空間性和時間性[3]。與基因組所具有的靜態(tài)實體的特點不同，轉錄組是受多種因子調控的，因此轉錄組是有機體基因組和外部物理特征的動態(tài)聯(lián)系[4]。本研究利用二代高通量測序對細粒棘球絳蟲原頭蚴轉錄組進行深入測序，豐富了細粒棘球絳蟲的基礎資料。

最后測序得到的經過過濾(將錯誤或存在雜質的數(shù)據(jù)去除)的數(shù)據(jù)達到2G。我們對數(shù)據(jù)進行從頭拼接,通過拼接獲得18 569個contig，這些contig的總長度為71 329 bp，contig平均長度為384 bp，最小的contig長度為201 bp，最大contig長度為4 618 bp，N50(覆蓋50%所有核苷酸的最大序列重疊群的長度)為384 bp。對轉錄組contig的長度特征進行分析，其中長度在200～300 bp范圍內的contig數(shù)目有10 265條，占總數(shù)的55%，從上述結果可以看出拼接的片段長度與預測相符，N50的值也大于300 bp說明測序結果較為完整，這次的測序結果能夠較準確的說明細粒棘球絳蟲在原頭蚴階段的表達情況，另一方面本次表達譜測序也能為將來的基因組注釋提供基因鑒定的基礎。

我們用從頭拼接的方法得到9 027條unigene，這9 027條unigene中有7 441條unigene獲得同源比對信息，在這些比對信息中，其中2 984條unigene與中華肝吸蟲(Clonorchissinensis)同源占總數(shù)的40.71%，其次與曼氏血吸蟲和日本血吸蟲的同源性也較高，分別占總數(shù)的27.01%和8.8%。值得注意的是在這些基因注釋信息中有98條與NCBI數(shù)據(jù)庫中的細粒棘球絳蟲原頭蚴信息吻合，這其中有62條(63.27%)匹配率達90%以上，這說明基因拼接過程基本準確。比對的信息中有21條未知蛋白，這些蛋白中17條和血吸蟲相匹配，且同源性高，這些蛋白可能是寄生蟲的特有蛋白，值得更深一步研究。還有1 588條轉錄本未匹配到基因信息可能是由于轉錄組序列偏單過短，在魏利斌等的芝麻發(fā)育轉錄組研究中未匹配到信息的轉錄組中長度范圍在100～500 bp的占98.07%，在水稻[5]、茶樹[6]和粉虱[7]的轉錄組研究中有相同的現(xiàn)象出現(xiàn)。由此可以看出轉錄組注釋信息的多少和轉錄組拼接后unigene的長短密切相關，unigene越長可注釋的信息越多，注釋結果也更據(jù)可靠性細粒棘球絳蟲原頭蚴基因信息的匱乏，目前細粒棘球絳蟲原頭蚴基因組的研究才剛剛起步，數(shù)據(jù)庫中與細粒棘球絳蟲相關的基因缺乏，也是造成注釋困難的主要原因之一，盡快開展細粒棘球絳蟲的基因組學研究有助于我們對細粒棘球絳蟲生活史的了解。細粒棘球絳蟲原頭蚴轉錄組的特異性，與果蠅、擬南芥、線蟲等模式生物相比，細粒棘球絳蟲的基因組研究落后，可供參考的基因組信息較少，且細粒棘球絳蟲在生物學上與日本血吸蟲有一定距離，所以這其中也有一些基因可能是細粒棘球絳蟲原頭蚴的特屬基因。

我們在GO分析中發(fā)現(xiàn)共有10 550條unigene與數(shù)據(jù)庫中的基因有較高同源性，這其中存在一個unigene對應多個數(shù)據(jù)庫基因的現(xiàn)象。在分子功能中大部分的分子都是起連接作用和催化活性，這可能是由于在原頭蚴階段，寄生蟲已經穩(wěn)定寄生于宿主體內，原頭蚴階段需要形成與肝臟相連接的包囊，且寄生蟲與外界穩(wěn)定連接并獲取營養(yǎng)，寄生蟲在這一階段最主要的生活是寄生，需要從宿主獲取營養(yǎng)并進行一系列的代謝活動所以七催化活性的蛋白也占大多數(shù)。而起發(fā)育作用的基因只占0.3%。另一需要注意的是其轉運功能的基因有122條，在本研究組前期對細粒棘球絳蟲原頭蚴的囊液蛋白質組研究就發(fā)現(xiàn)囊液中含有大量的轉鐵蛋白和白蛋白，前者轉運鐵離子，后者可以與難溶的小分子結合并作為營養(yǎng)物質。寄生蟲靠汲取宿主的養(yǎng)分生活，若能截斷這種轉運則很有可能截斷寄生蟲的營養(yǎng)源，所以從這個角度看，這些具有轉運功能的蛋白具有很大研究價值。

通過與KEGG數(shù)據(jù)庫進行比對分析，細粒棘球絳蟲原頭蚴的轉錄組中有4 731條unigene得到注釋，這4 731個得到注釋的基因位于241條代謝通路中，這些代謝通路分別與代謝過程，基因信息過程，環(huán)境相關過程，細胞過程及與人類疾病相關，其中代謝過程的代謝最多。

[1]Thompson RCA. Biology and systematics ofEchinococcus[J]. In Thompson RCA, Limbery AJ (eds):Echinococcus and Hydatid Disease. Wallingford: CAB International, 1995:1-50.

[2]Zhao WX. Human parasitology[M]. Beijing:People’s Medical Publishing House, 1983: 514. (in Chinese) 趙慰先.人體寄生蟲學[M].北京:人民衛(wèi)生出版社,1983:514.

[3]Lander ES, Linton LM, Birren B, et al. Initial sequencing and analysis of the human genome[J]. Nature, 2001, 409(6822):860-921. DOI:10.1038/35057062

[4]Wu Q, Sun C, Chen SL. The application of transcriptome in medicinal plant research[J]. World Sci Technol/Modernizat Traditional Chin Med Materia Medica, 2010, 12(3):457-462. (in Chinese) 吳瓊,孫超,陳士林，等. 轉錄組學在藥用植物研究中的應用[J].世界科學技術中醫(yī)藥現(xiàn)代化,2010,12(3):457-462.

[5]Lu TT, Lu GJ, Fan DL, et al. Function annotation of the rice transcriptome at single-nucleotide resolution by RNA-seq[J]. Genome Res, 2010, 20: 1238-1249. DOI:10.1101/gr.106120

[6]Shi CY, Yang H, Wei CL, et al. Deep sequencing of theCamelliasinensistranscriptome revealed candidate genes for major metabolic pathways of tea-specific compounds[J]. BMC Genomics, 2011, 12:131. DOI:10.1186/1471-2164-12-131

[7]Wang XW, Luan JB, Li JM, et al. De novo characterization of a whitefly transcriptome and analysis of its gene expression during development[J]. BMC Genomics, 2010, 11:400. DOI:10.1186/1471-2164-11-400

mRNA sequencing and transcriptome characteristic ofEchinococcusgranulosus

JU Yan1,2,3，LI Zi-hua1,3，WANG Ya-na1,3，ZHAO Jia-qing1,3，ZHU Ming-xing1,3，LI Jun-liang1,3，Zhao Wei1,3

(1.CenterofMedicalScience,NingxiaMedicalUniversity,Yinchuan750004,China;2.NingxiaCentersforDiseasesPrevention&Control,Yinchuan750004,China;3.DepartmentofGeneticsandCellBiology,NingxiaMedicalUniversity,Yinchuan750004,China)

The objective of this study was to establish preliminarily expression profile database of protoscolex and comprehend gene expression and protein composition of protoscolex by mRNA sequencing of protoscolex and analysis of expression profile, and laid the foundation for the research of relationship between the parasites and hosts. It would be the theoretical basis for finding the new methods of diagnosis, screening new drug targets and selecting the vaccine candidate. TRIZOL method was used to extract total RNA from human protoscolex. Illumina’s solexa sequencing platforms for RNA sequencing gene annotation and the analysis of metabolic pathways were performed by comparing unigene with non redundant database of NCBI, the gene ontology database and KEGG database gene. As a result, we obtained 18 569 contigs, the total length of the contig was 71 329 bp and the average length was 384 bp, the smallest contig was 201 bp and the largest one was 4 618 bp. The N50 was 384 bp and the predicted unigene was up to 9 029. Comparing the unigene with non redundant database of NCBI, the gene ontology database and KEGG database gene, we got 7 441 annotation comparison information in NR database of NCBI, 10 550 information from gene ontology and 4 731 annotation in KEGG. The 4 731 annotated genes from KEGG was mapped in 241 metabolic pathways, which was associated with metabolism process, the genetic information process, environment relative processes, cellular processes, and human disease respectively.

protoscolex; transcriptome; biological information analysis; expression profile

Zhao Wei, Email:zw-6915@163.com

趙巍，Email： zw-6915@163.com

1.寧夏醫(yī)科大學醫(yī)學科學技術研究中心，銀川 750004； 2.寧夏疾病預防控制中心，銀川 750004； 3.寧夏醫(yī)科大學醫(yī)學遺傳學與細胞生物學系，銀川 750004

10.3969/cjz.j.issn.1002-2694.2015.01.005

R383.3

A

1002-2694(2015)01-0021-05

2014-03-11；

2014-09-22

寧夏教育廳基金(NGY2010041)資助

Supported by the Key Fund Project of Department of Education of Ningxia Hui Autonomous Region (No. NGY2010041)