細(xì)粒棘球絳蟲原頭蚴mRNA測序及表達(dá)譜分析

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

細(xì)粒棘球絳蟲原頭蚴mRNA測序及表達(dá)譜分析

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

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

原頭蚴；轉(zhuǎn)錄組學(xué)；生物信息分析；表達(dá)譜

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

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

1 材料與方法

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

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

1.3 細(xì)粒棘球絳蟲原頭蚴的轉(zhuǎn)錄組測序 細(xì)粒棘球絳蟲原頭蚴的轉(zhuǎn)錄組測序委托北京基因組完成

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

2 結(jié) 果

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

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

Tab.1 Data assembly for contig digital transcriptome of protoscolex

2.2 細(xì)粒棘球絳蟲原頭蚴轉(zhuǎn)錄物功能注釋及分類 經(jīng)過augustus、pasa、cufflinks三種軟件對contig進(jìn)行基因預(yù)測，最終預(yù)測得到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%。不知道功能的假設(shè)蛋白和未命名蛋白分別為1 228條和68條。在這些匹配的同源基因中與細(xì)粒棘球絳蟲相吻合的98條基因，還有69條基因與多房棘球絳蟲同源。

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

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

3 討 論

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

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

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

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

通過與KEGG數(shù)據(jù)庫進(jìn)行比對分析，細(xì)粒棘球絳蟲原頭蚴的轉(zhuǎn)錄組中有4 731條unigene得到注釋，這4 731個(gè)得到注釋的基因位于241條代謝通路中，這些代謝通路分別與代謝過程，基因信息過程，環(huán)境相關(guān)過程，細(xì)胞過程及與人類疾病相關(guā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) 趙慰先.人體寄生蟲學(xué)[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) 吳瓊,孫超,陳士林，等. 轉(zhuǎn)錄組學(xué)在藥用植物研究中的應(yīng)用[J].世界科學(xué)技術(shù)中醫(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ī)科大學(xué)醫(yī)學(xué)科學(xué)技術(shù)研究中心，銀川 750004； 2.寧夏疾病預(yù)防控制中心，銀川 750004； 3.寧夏醫(yī)科大學(xué)醫(yī)學(xué)遺傳學(xué)與細(xì)胞生物學(xué)系，銀川 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)