艾美游仆蟲大核基因組與轉錄組測序及結構特征

何 燕 尹家奇 生 欣

(遵義醫科大學生物化學與分子生物學教研室, 遵義 563000)

游仆蟲是進化上較為高等的腹毛目纖毛蟲, 為一種常見的水生浮游生物, 以水中微生物為食。具有十分復雜的皮層微管骨架系統和背腹分化的纖毛器。常被作為研究微管裝配特征、纖毛結構與功能、纖毛形態發生、纖毛基因定位和功能研究的重要模式生物。此外, 游仆蟲不僅同時具有負責營養的大核與負責生殖的小核, 而且具有其獨特的大核基因結構特征, 如微染色體、特殊的終止密碼子和高頻率的編程性核糖體移碼(Programmed ribosomal frameshifting, PRF)等。因此, 近年來, 游仆蟲受到研究者的廣泛關注, 目前, 八肋游仆蟲(Euplotes octocarinatus)、厚游仆蟲(Euplotes crassus)和扇形游仆蟲(Euplotes vannus)等的基因組數據庫已經建立, 為進一步揭秘這種水生浮游生物的基因、細胞和環境生物學特征提供了不可或缺的基礎資料。

艾美游仆蟲(Euplotes amieti)為一類較大型淡水種游仆蟲。活體為(88—125) μm×(55—78) μm,呈不對稱的卵圓形, 背腹扁平, 腹面觀右緣較左緣膨出, 體前有一明顯的領口區開闊呈三角形約占體長的活體時口側膜不易觀察到橫棘毛間有列短而不規則的嵴[1]。本課題組前期對艾美游仆蟲微管類細胞骨架、纖毛形態發生過程和γ-微管蛋白在纖毛裝配過程中的作用進行了研究, 但由于缺乏該種游仆蟲基因組和轉錄組數據, 限制了對其相關基因功能的研究, 因此, 本研究通過提取艾美游仆蟲大核基因組DNA和mRNA, 對其大核基因組DNA和轉錄組進行了測序和基因注釋, 旨在篩選出與微管和纖毛相關蛋白基因, 分析其基因結構, 為進一步探索其功能提供基礎資料。

1 材料與方法

1.1 材料

本實驗所用的艾美游仆蟲(Euplotes amieti)釆自上海市青浦區近郊的農田水塘中。本實驗所用DNA提取試劑盒, RNA提取試劑盒購買于TaKaRa公司。基因組測序及注釋由北京諾禾致源科技有限公司完成。轉錄組測序由上海派森諾生物科技股份有限公司完成。

1.2 方法

艾美游仆蟲的培養及收集 本研究所用的腹毛類纖毛蟲艾美游仆蟲使用長梭綠藻(Chlorogoium elongatum)喂食。在收集艾美游仆蟲細胞前,先進行饑餓處理5—7d。然后用孔徑較小的紗布過濾饑餓后的細胞, 除去較大的雜質。隨后用定性濾紙濃縮收集蟲體。同時除去較小雜質。最后用純凈水將蟲體收集到離心管中。4000 r/min離心5 min收集蟲體, 棄掉上清備用。

艾美游仆蟲DNA的提取及基因組測序 參照TaKaRa公司DNA提取試劑盒說明書, 提取的DNA采用瓊脂糖凝膠電泳分析DNA的純度和完整性。Nanodrop檢測DNA的純度(A260/A280比值);Qubit對DNA濃度進行精確定量。檢驗合格的DNA樣品通過Covaris破碎機隨機打斷成長度為350 bp的片段。采用 NEB Next?Ultra DNA Library Prep Kit(NEB, USA)進行建庫, 構建好的文庫通過Illumina NovaSeq PE150進行測序。然后再采用SPAdes的“careful”模式進行組裝。隨后, 利用CAP3將SPAdes 的拼接結果進行融合。考慮其一些reads質量較低, 采用blat將長度小于500 bp的序列與長度大于500 bp的序列比對, 去掉一致性百分比≥90%, 覆蓋度≥80%的序列、細菌、古菌污染的DNA、無端粒序列及線粒體基因組序列, 11條長度小于100 bp的序列等共去除9508條。最終得到Clean Contigs用于全基因組注釋。

RNA 的提取及轉錄組測序 RNA 的提取參照TaKaRa公司RNA提取試劑盒說明書, 將提取的RNA采用離子打斷的方式打斷成300 bp左右的片段。構建好大小為450 bp文庫。采用第二代測序技術, 基于 Illumina HiSeq 測序平臺, 對這些文庫進行雙末端 (Paired-end, PE)測序。Cutadapt 去除 3′端帶接頭的序列和平均質量分數低于Q20 的Reads, 得到的CleanData 用Trinity 軟件與基因組比對, 然后根據比對結果, 使用Trinity的Genomeguid模式進行從頭拼接、聚類。挑選最長的轉錄本作為Unigene, 并進行后續基因功能注釋和基因結構預測等。

基因注釋 重復序列注釋方法使用Repeatmasker和repeatproteinmas軟件對重復序列數據庫RepBase庫進行同源序列比對, 識別與已知重復序列相似的序列。從頭預測使用LTR_FINDER和RepeatScout, RepeatModele等軟件首先建立de novo重復序列庫, 將de novo預測出來的重復序列庫與同源重復序列數據庫Repbase進行整合, 再用Repeat-Masker軟件對艾美游仆蟲基因組進行repeat注釋。

基因的結構預測, 主要通過同源預測,de novo預測和其他證據支持的預測。同源預測的方法是將已知的同源物種的編碼蛋白序列與新物種的基因組序列進行比對(同源物種個數通常8—10個), 通過blast, genewise等比對軟件預測基因組中的基因結構。de novo預測使用依賴于基因組序列數據統計學特征(如密碼子頻率和外顯子-內含子分布)的軟件來預測基因結構, 常用的軟件有Augustus、GlimmerHMM和SNAP等。再結合轉錄組比對數據, 使用EVidenceModeler(EVM)整合軟件,將各種方法預測得到的基因集整合成一個非冗余的, 更加完整的基因集。最后, 使用PASA(http://pasa.sourceforge.net/), 結合轉錄組組裝結果, 對EVM的注釋結果進行校正, 加入UTR及可變剪切等信息, 得到最終的基因集。

基因功能注釋, 是將基因結構注釋得到的基因集, 利用比對軟件與已知蛋白數據庫(SwissProt、Nr、Pfam、KEGG和InterPro)等比對, 得到基因的功能信息。注釋物種包括三偽尖毛蟲(Oxytricha trifallax,O. trifallax)、線蟲(Caenorhabditis elegans,C. elegans)、浮萍棘尾蟲(Stylonychia lemnae,S.lemnae)、扇形游仆蟲(Euplotes vannus,E. vannus)、草履蟲(Paramecium duboscqui,P. duboscqui)、八肋游仆蟲(Euplotes octocarinatus,E. octocarinatus)、四膜蟲(Tetrahymena utriculariae,T. utriculariae)、杜氏利什曼蟲(Leishmania donovani,L. donovani)、有孔蟲(Reticulomyxa filosa,R. filosa)和艾美游仆蟲(Euplotes amieti,E. amieti)等[2—15]。

非編碼RNA的注釋包括tRNA、rRNA、miRNA和snRNA。根據tRNA的結構特征, 利用tRNAscan-SE(http://lowelab.ucsc.edu/tRNAscan-SE/)軟件來尋找基因組中的tRNA序列; 由于rRNA具有高度的保守性, 因此可以選擇近緣物種的rRNA序列作為參考序列, 通過blast比對來尋找基因組中的rRNA; 利用Rfam家族的協方差模型, 采用Rfam自帶的INFERNAL(http://infernal.janelia.org/)軟件可預測基因組上的miRNA和snRNA序列信息。最終得到艾美游仆蟲基因組的ncRNA信息。

2 結果

2.1 艾美游仆蟲基因組和轉錄組測序數據分析

鑒于游仆蟲基因組特點, 利用二代測序技術Illumina測序平臺對艾美游仆蟲基因組的200 bp小片段數據質控和建庫進行測序, 共得到10.92 Gb 原始數據。用 SPAdes拼接軟件進行組裝, 初始組裝產生了64836條Contigs。由于在提取游仆蟲基因組DNA的過程中, 無法將其線粒體基因組分離去除,所以拼接結果中也可能混有線粒體DNA的污染和其他雜菌污染, 初始數據經去除細菌、線粒體基因組DNA污染, 去除＜100 bp短序列后, 最終得到50287條艾美游仆蟲基因組序列。這些基因序列的GC含量較低, 為 31%, 與其他纖毛蟲類似。此外,兩端同時具有端粒的微染色體數量為27542條, 占54.76%, 只含有一端端粒的基因數量為6118條(表 1)。與其他游仆蟲相比, 艾美游仆蟲大核基因組均大于八肋游仆蟲(88.9 Mb)和扇形游仆蟲(85.1 Mb)。但微染色體數量小于八肋游仆蟲(29413)與扇形游仆蟲(37501)。

表1 艾美游仆蟲基因組測序數據統計Tab. 1 Genomic sequencing data statistics of Euplotes amieti

轉錄組測序結果去除低質量的Raw Reads, 最終得到68900504條Clean Reads, 經Trinity拼接后得到60691條Transcripts, 其平均長度和N50值分別為1326.6和 1759 bp。所有的transcripts進一步聚類后獲得38588條Unigenes, Unigenes的平均長度和N50值分別為1189.9和1643 bp(表 2), 將38588條轉錄組序列通過NCBI中進行BLASTX進行比對, 并進一步通過查找開放閱讀框, 以E值2×10-5作為標準, 發現其中2%—3%基因發生了編程性移碼。

表2 艾美游仆蟲轉錄組拼接序列Tab. 2 Transcriptome splicing sequence of Euplotes amieti

用alveolata_odb10 進行BUSCO評估測序結果的完整度, 其完整性評估基于其他原生動物的數據庫建立, 且通過與目前已知的八肋游仆蟲和扇形游仆蟲基因組比較, 在基因組大小、微染色體數量、基因的GC含量、N50和基因長度分布等方便均與已報道游仆蟲類似, 表明基因組測序質量較好。

基因組通過結構預測得到了27650個基因, 其中96.5%的基因能夠被預測出功能。分別將基因組27650個基因與轉錄組38588條Unigenes對比到已知的NR、GO、KEGG、SwissProt、Pfam、InterPro和eggNOG數據庫, 其中有26673條大核基因獲得注釋, 在各數據庫中注釋成功的基因數與轉錄本數分別見圖 1A和1B。功能分析顯示, 在eggNOG數據庫中注釋最多的功能類為信號轉導機制(1790條)和轉錄后修飾、蛋白折疊和分子伴侶(1157條)。而GO注釋主要富集在生殖詞條中(1599條); KEGG 分析富集主要在外界環境信息的信號轉導(962條)、細胞內物質運輸與代謝分解、分泌、細胞增殖與死亡和環境適應等。這些功能與游仆蟲捕食浮游生物并進行細胞內消化與營養物質的運輸、適應外環境變化、分裂與結合生殖和個體間信息傳遞等生命活動密切相關。

圖1 基因組及轉錄組測序結果Fig. 1 Genome and transcriptome sequencing results

2.2 基因結構分析

將艾美游仆蟲基因組測序得到的50287條大核基因組進行重復序列的注釋得到含有19.31%的repeat序列。結合轉錄組測序結果進行基因結構注釋;注釋物種包括O. trifallax、C. elegans、S. lemnae、E. vannus、P. duboscqui、E. octocarinatus、T. utriculariae、L. donovani、R. filosa和E. amieti等。結果顯示艾美游仆蟲的平均transcript長度為680.21 bp,其中, 平均CDS為560.11 bp, 是幾種纖毛蟲中最小的(圖 2A); 而每個基因的平均外顯子數為1.8, 是除八肋游仆蟲外最大的(圖 2B)。此外, 其平均外顯子與內含子大小分別為311.69和150.70 bp。與以上幾種纖毛蟲相比, 盡管艾美游仆蟲的平均trancript、CDS和外顯子長度均較小, 而平均內含子長度與平均外顯子數量卻較大(圖 2C)。

圖2 基因結構特征比較Fig. 2 Comparison of gene structures

通過對轉錄組38588條Unigenes分析發現艾美游仆蟲中同樣存在較高頻率的編程性核糖體移碼現象, 發生率為2%—3%。其中絕大多數為+1PRF基因, 發生移碼的位置多見于滑動序列的終止密碼子TAR(R為A或G); 也有部分為+2PRF基因, 多發生在滑動序列的終止密碼子TAR(R為A或G), 同時跳過T和A, 這些可以證明艾美游仆蟲也存在編程性移碼突變的現象。此外, 艾美游仆蟲的終止密碼子既可以做終止密碼子也可以編碼蛋白質。將艾美游仆蟲與其他幾種纖毛蟲的終止密碼子做比較, 艾美游仆蟲與八肋游仆蟲一樣, UAA 和 UAG常作為終止密碼子, 而UGA編碼半胱氨酸和硒代半胱氨酸。而在四膜蟲草履蟲和尖毛蟲中, 只有 UGA 作為終止密碼子, 而UAA和 UAG 則編碼谷氨酰胺(圖 3)。

圖3 艾美游仆蟲與多種真核生物基因組比較(進化樹基于18S rRNA繪制)Fig. 3 CComparison of representative eukaryotic genome. The tree was constructed based on the sequences of 18S rRNA genes

2.3 非編碼RNA及轉錄因子分析

利用tRNAscan-SE和Rfam等軟件對基因組中的非編碼RNA序列信息注釋; 結果顯示miRNA的基因有23個, 平均長度為125.78 bp; tRNA有105個,平均長度為74.32 bp; rRNA有56個, 其中, 18S 26個和28S 2個。平均長度分別為104.68和104.81 bp;snRNA 21個其中包括CD-box 1個, splicing 9個, 平均長度分別為263和145.22 bp(表 3)。轉錄因子(Transcription factor, TF)是一類能與基因5′端上游特定序列專一性結合, 并與RNA聚合酶Ⅱ形成轉錄起始復合體, 共同參與轉錄起始的過程的蛋白質分子。將植物和動物與PlantTFDB和AnimalTFDB數據庫比較, 從而預測得到轉錄因子及轉錄因子所屬的家族信息。在所有預測的轉錄因子家族中, zf-C2H2的Count數最多, 達139條, 其次分別為ZBTB 116條、MYB 89條、bHLH 55條等(圖 4)。

圖4 轉錄因子家族統計圖Fig. 4 Transcription factor family statistic

表3 非編碼RNA序列注釋Tab. 3 Noncoding RNA sequence annotation

2.4 基因組及轉錄組基因驗證

將基因組及轉錄組測序的結果序列中隨機取50個基因進行PCR驗證。基因組驗證結果如圖 5A和5B所示, 所用Marker大小為2000 bp, 圖中只呈現出1—48個基因的驗證條帶, 基因組測序結果圖 5A和5B可以看出只有1個基因RABEP2未P出, 其余各個基因條帶大小與測序的基因大小一致。轉錄組驗證結果如圖 5C和5D所示, 圖中除了POC1A外, 其余各基因條帶均與測序結果一致。

圖5 基因組及轉錄組結果驗證Fig. 5 Genome and transcriptome validation

3 討論

3.1 艾美游仆蟲具有游仆蟲屬共有的基因組特征

原生動物纖毛蟲種類繁多, 不同類型的纖毛蟲在基因表達調控、大分子相互作用和應激反應等方面呈現極大的多樣性特征。如游仆蟲中存在的編程性移碼現象[16], 尖毛蟲中出現的由長非編碼RNA指導的大規模基因亂序和重排等[17]。因此, 揭示這些存在于不同類群纖毛蟲中的分子與細胞生物學特征對于研究纖毛蟲的系統發育和進化關系具有重要的意義。然而, 基因組信息的缺失也極大地限制了纖毛蟲分子和細胞生物學水平的多樣性研究。目前為止, 已經獲得基因組測序資料的纖毛蟲有20多種, 而游仆蟲只有3種, 分別為厚游仆蟲、八肋游仆蟲和扇形游仆蟲, 這些研究揭示了游仆蟲的共有特征: 含有微染色體, 即一個基因一個染色體; 終止密碼密碼子的重新分配; 程序性核糖體框架轉移和對環境應激源的強烈抵抗力等[18]。

本研究對艾美游仆蟲基因組和轉錄組測序組裝后, 最終成功注釋了50287條大核基因組序列, 其N50為2774 bp, GC含量為31%, 54.8%的基因含有兩端端粒, 與以上3種游仆蟲大核基因組的組裝結果相似[19,20]。這表明艾美游仆蟲測序與拼裝結果較好, 且符合游仆蟲較低GC含量和微染色體的基因組基本結構特點[19]。同時也存在一定的差異, 其GC含量和微染色體百分比均介于厚游仆蟲與八肋游仆蟲之間, 且艾美游仆蟲中長度小于500 bp的基因與含有一個端粒或無端粒的基因數量均較八肋游仆蟲和扇形游仆蟲多。排除細菌與小核基因組污染, 推測可能是艾美游仆蟲中還含有大量“聯合微染色體”, 這些染色體中含有2個以上基因[20]。盡管目前有研究“聯合微染色體”上不同基因的表達方向的可能相同或相反, 但微染色體所占的比例是否與生物進化和物種親緣性相關還未見報道, 但可以確定的是, 不同游仆蟲之間存在多樣化的基因組特征將為確定種屬之間的進化地位提供新的依據。另一方面, 基因結構分析顯示基因組序列中含

有19.31%的重復序列, 這種基因的高度串聯重復性最早在1983年在鬃棘尾蟲小(Stylonychia pustulata)核中被報道[21—24], 此后這種現象在纖毛蟲大核基因組中被廣泛報道, 可見高度重復性是纖毛蟲的基因組特征之一。另外, 與幾種已知的腹毛類纖毛蟲基因結構相比, 艾美游仆蟲也具有內含子較短的特征,平均長度僅為150 bp, 與八肋游仆蟲的189 bp類似,而艾美游仆蟲平均CDS區和外顯子大小均小于其他物種, 表明艾美游仆蟲基因同樣具有高度片段化特征。

3.2 艾美游仆蟲基因組同樣具有編程性移碼

編程性核糖體移碼是一種重編碼事件, 它是指翻譯中的核糖體能夠在mRNA 上的特定位置, 從起始的0讀框轉換到+1或者-1讀框, 然后繼續進行翻譯的現象。這種現象的發生是可調控的, 發生頻率高達80%。目前已經報道的核糖體編程性移碼信號常見于3個主要元件, 即七核苷酸滑動序列(5′—AAA-UAR-V-3′ R為A 或G; V≠U), 及其上游都有SD相似序列CAAGAA, 和5—12個核苷酸組成的間隔序列以及假結(Pseudonotes)或莖環; 其可產生牽引效應將XXX(AAA、GGG或UUU)和ZZZ(AAA或UUU)引入P和A位點然后引起核糖體結構重排[16]。本研究中2%—3%的編程性移碼基因, 與目前報道的八肋游仆蟲和扇形游仆蟲中3.5%和2.8%的編程性移碼基因發生率相當[10,25], 且其中絕大多數為+1PRF。此外, 與其他纖毛蟲相比, 艾美游仆蟲也發生了終止密碼子的重新分配, 與八肋游仆蟲和扇形游仆蟲一樣, UAA 和 UAG作為終止密碼子, 而UGA編碼半胱氨酸和硒代半胱氨酸。相比之下, 在四膜蟲、草履蟲和尖毛蟲中, 只有 UGA 作為終止密碼子, 而UAA和 UAG 則編碼谷氨酰胺。以上結果進一步證實游仆蟲中存在著高頻率的編程性移碼現象, 這種現象是在翻譯水平上進行基因表達調控的獨特方式, 對于游仆蟲有限的大核基因信息來說, 通過編程性移框與終止密碼子的重新分配, 能夠產生多樣化的蛋白表型, 是其適應外界環境的變化的分子基礎。目前, 部分研究顯示游仆蟲中的編程性移碼在不同類型之間還存在轉變現象, 這種現象存在簡約和精致的調控機制[16], 這種調控機制的揭示將有利于闡述游仆蟲高頻率編程性移碼與基因表達調控和環境適應性之間的關系[22—24]。

3.3 艾美游仆蟲的基因表達與適應環境密切相關

另外, 轉錄組測序獲得了60691個轉錄本, 分別是扇形游仆蟲與八肋游仆蟲的1.5倍和的3倍。其中, 38588個基因被成功注釋, mRNA序列平均長度為1189.99 bp, 96.5%的基因被成功預測功能。通過基因功能分析顯示, 基因組絕大多數基因富集到信號轉導、轉錄后修飾和蛋白質折疊, 而GO與KEGG功能分析顯示其轉錄本顯著富集于生殖、單細胞過程、外界環境信息的信號轉導、細胞內物質運輸與代謝分解、分泌、細胞增殖與死亡和環境適應等。與八肋游仆蟲和扇形游仆蟲的功能富集上具有較大差別。這提示艾美游仆蟲在營養生長期,除了編碼游仆蟲屬特有的功能蛋白外, 還編碼大量蛋白質以適應外界環境變化。盡管本研究并未對不同溫度、離子濃度及pH等應激狀態下的游仆蟲進行轉錄組測序, 從扇形游仆蟲在營養條件、極端溫度和鹽濃度等條件下的基因表達情況可見, 游仆蟲中存在具有應對不同理化刺激的基因簇, 這些基因的表達對于其適應外界環境具有重要意義。而不同游仆蟲中應對不同刺激基因簇是否一致還有待進一步探索。此外, 非編碼RNA與轉錄因子分析顯示, 在艾美游仆蟲中存在一定量的非編碼RNA,其中miRNA的數量為23個, 未見長非編碼RNA, 這在游仆蟲中還未見報道。有研究顯示長非編碼RNA在尖毛蟲中具有指導大規模基因亂序與重排等作用[25], 因此, 這些miRNA是否也參與了游仆蟲基因高度片段化和移碼突變等過程還有待進一步研究。而在預測的轉錄因子家族中, 鋅指蛋白基因ZF-C2H2與ZBTB的數量最多, 盡管有研究顯示轉錄因子E2Fl1在嗜熱四膜蟲中參與了減數分裂[26], 但這2種轉錄因子在纖毛蟲中的功能研究還未見報道,在植物與細菌中的研究顯示, ZF-C2H2與ZBTB參與了真核生物生長發育及逆境脅迫的耐受等[27], 表明艾美游仆蟲可能通過這兩類轉錄因子對基因表達調控, 從而適應外界環境的變化。