周錢森，任憲云，史鯤鵬，于振興，劉 萍，李 健

（1.水產(chǎn)科學(xué)國家級實驗教學(xué)示范中心（上海海洋大學(xué)），上海 201306；2.中國水產(chǎn)科學(xué)研究院黃海水產(chǎn)研究所，農(nóng)業(yè)農(nóng)村部海洋漁業(yè)可持續(xù)發(fā)展重點實驗室，青島 266071；3.青島海洋科學(xué)與技術(shù)試點國家實驗室，海洋漁業(yè)科學(xué)與食物產(chǎn)出過程功能實驗室，青島 266071）

隨著高通量測序技術(shù)的發(fā)展，轉(zhuǎn)錄組測序在基因表達(dá)水平分析和差異表達(dá)分析、新基因的挖掘、尋找單核苷酸多態(tài)性及應(yīng)用、基因功能注釋等方面都有所應(yīng)用，在生物研究中發(fā)揮重要作用［1］。目前為止，大多數(shù)轉(zhuǎn)錄組測序研究利用以illumina平臺為代表的二代測序技術(shù)，其測序長度遠(yuǎn)遠(yuǎn)低于真核生物RNA長度［2－4］。以PacBio平臺為代表的第三代測序技術(shù)因其強大的讀長優(yōu)勢，有效地解決了二代測序技術(shù)由于讀長限制引起拼接不完整的問題［5］。PacBio平臺單分子實時測序技術(shù)，又稱SMRT（single molecule realtime）測序，SMRT cell可以在高GC含量區(qū)域完全跨過，確保序列的均勻覆蓋度，并且DNA建庫過程中，不需要進(jìn)行PCR擴增，避免了PCR冗余和覆蓋度不均一的問題，數(shù)據(jù)精準(zhǔn)度可以達(dá)到99.9%［6－7］。PacBio SMRT測序借助其諸多優(yōu)勢，近些年在水產(chǎn)領(lǐng)域研究中得到廣泛應(yīng)用，張金勇等［8］對金烏賊（Sepiaesculenta）采用全長轉(zhuǎn)錄組測序獲得高質(zhì)量的生物學(xué)數(shù)據(jù)信息，篩選SSR位點并分析了其分布及組成特征?？讎[蘭等［9］通過對竹筴魚（Trachurusjaponicus）高通量測序，篩選獲得27個具有多態(tài)性的微衛(wèi)星標(biāo)記，利用一個竹筴魚群體對通過篩選的微衛(wèi)星標(biāo)記的種群遺傳學(xué)特征進(jìn)行評價。ZHANG等［10］通過三代與二代測序技術(shù)相結(jié)合挖掘施氏鱘（Acipenserschrenckii）性腺中早期配子發(fā)生的相關(guān)基因，更好地了解其生殖調(diào)控機制。本研究采用三代測序技術(shù)的PacBio SMRT平臺對錦繡龍蝦（Panulirusornatus）進(jìn)行全長轉(zhuǎn)錄組測序，通過生物信息學(xué)方法對所測的序列進(jìn)行拼接、分類、功能注釋和代謝途徑分析，獲得錦繡龍蝦豐富的序列信息，旨在為進(jìn)一步挖掘錦繡龍蝦相關(guān)功能基因、基因組學(xué)及開發(fā)分子標(biāo)記等研究奠定基礎(chǔ)。

錦繡龍蝦隸屬于節(jié)肢動物門（Arthropoda），甲殼綱（Crustacea），十足目（Decapoda），龍蝦科（Palinuridae），龍蝦屬，在已知的19個種類的龍蝦中錦繡龍蝦是其中較為珍貴的優(yōu)質(zhì)水產(chǎn)品，其體型較大、營養(yǎng)豐富、脂肪含量低、味道鮮美［11－12］，主要分布在熱帶印度洋東部、東南亞、澳大利亞和西太平洋地區(qū)［13］，生長速度快于波紋龍蝦（P.homarus）、中國龍蝦（P.stimpsoni）、雜色龍 蝦（P.versicolor）和 黃 斑 龍 蝦（P.polyphagus）等［14］。作為寶貴的海洋漁業(yè)資源，錦繡龍蝦因遭受高強度捕撈，資源量急劇減少。目前，錦繡龍蝦育苗尚未取得成功，亟需開展人工繁殖相關(guān)研究以保護(hù)和增殖資源［15］。錦繡龍蝦基因組測序亦未完成，其基因序列信息還比較匱乏，相關(guān)的分子遺傳基礎(chǔ)也很薄弱，因此本研究利用高通量測序技術(shù)對錦繡龍蝦全長轉(zhuǎn)錄組測序，以期為后期相關(guān)功能基因的研究以及種質(zhì)繁育等提供理論支持和數(shù)據(jù)支撐。

1 材料與方法

1.1 材料

實驗所用錦繡龍蝦購于海南省瓊海市博鰲鎮(zhèn)永賀生物科技有限責(zé)任公司，隨機挑選3只體表完整無明顯傷痕、活力好、規(guī)格整齊的龍蝦進(jìn)行取樣，分別取鰓、肝胰腺、肌肉、胃、腸、腦和心臟組織，迅速放入液氮中保存，用于后續(xù)的RNA提取。

1.2 RNA提取與質(zhì)控

采用Trizol法［16］分別提取錦繡龍蝦鰓、肝胰腺、肌肉、胃、腸、腦和心臟組織總RNA，利用瓊脂糖凝膠電泳分析RNA是否存在污染；通過Nanodrop測取OD260／280比值檢測RNA濃度；Qubit精確定量RNA濃度；Agilent 2100判斷RNA的完整性。選取符合標(biāo)準(zhǔn)的各組織RNA等量混合構(gòu)建測序文庫。

1.3 測序文庫構(gòu)建

將各組織混合后的RNA用Oligo（dT）的磁珠富集含有polyA的mRNA，然后通過SMARTer PCR cDNA Synthesis Kit反轉(zhuǎn)錄為cDNA，循環(huán)優(yōu)化選取最合適的條件，通過PCR大規(guī)模擴增利用磁珠篩選的片段，獲得一定數(shù)量的cDNA，全長cDNA進(jìn)行末端修復(fù)、損傷修復(fù)、連接SMRT啞鈴型接頭，構(gòu)建全長轉(zhuǎn)錄組文庫。核酸外切酶消化，剔除cDNA左右端未連接的序列，最終通過DNA聚合酶結(jié)合引物，形成完整的文庫。測序由北京諾禾致源科技股份有限公司完成。

1.4 全長轉(zhuǎn)錄組測序及序列分析

采用SMRTlink（PacBio）處理原始下機數(shù)據(jù)獲得子序列，校正子序列獲得環(huán)形一致性序列（circular consensus sequence，CCS），然后根據(jù)環(huán)形一致性序列是否含有3′端引物、5′端引物以及polyA尾將其分為全長序列與非全長序列；再對全長序列進(jìn)行聚類，獲得聚類一致序列［17］；最后對得到的全長序列進(jìn)行改良，獲得高質(zhì)量的一致序列用于后續(xù)分析。

1.5 全長轉(zhuǎn)錄組序列的功能注釋及結(jié)構(gòu)分析

利用CD-HIT軟件對獲得的改良一致序列去冗余，將錦繡龍蝦組裝所得單基因簇與公共數(shù)據(jù)庫進(jìn)行比對，通過序列相似程度得出具有最高相似性的蛋白，從而對該單基因簇進(jìn)行功能注釋，利用NCBI非冗余蛋白序列數(shù)據(jù)庫（nonredundant protein sequences，NR）、NCBI核酸序列數(shù)據(jù)庫（nucleotide sequences，Nt）［18］、蛋白質(zhì)真核同源數(shù)據(jù)庫（eukaryotic orthologous groups，KOG）［19］、蛋白質(zhì)序列數(shù)據(jù)庫（SwissProt protein sequence database，SwissProt）［20］、蛋白質(zhì)家族數(shù)據(jù)庫（protein families database，Pfam）［21］、基因功能描述分類系統(tǒng)（gene ontology，GO）［22］及京都基因與基因組百科全書（Kyoto encyclopedia of genes and genomes，KEGG）［23］。

2 結(jié)果與分析

2.1 全長轉(zhuǎn)錄組測序與組裝分析

2.1.1 測序結(jié)果與數(shù)據(jù)組裝

采用PacBio Sequel測序平臺對錦繡龍蝦鰓、肝胰腺、肌肉、胃、腸腦和心臟等組織混樣進(jìn)行全長轉(zhuǎn)錄組測序，共獲得39 828 440個子序列（71.1 Gb），平均子序列長度為1 786 bp，N50為2 492 bp。通過每個ZMW孔中的子序列得到CCS聚類之后得到的序列數(shù)為1 018 599個，序列平均長度為2 377 bp，N50為2 872 bp。同時含有3′引物和5′引物，以及3′引物前含有polyA尾的全長序列（full-length，F(xiàn)L）554 600個，全長非嵌合序列（full-length non-chimeric read，F(xiàn)LNC）543 552個，序列平均長度為2 055 bp，F(xiàn)LNC／CCS為53.36%。全長轉(zhuǎn)錄組得到改良后一致序列28 034個，序列平均長度為2 250 bp。

2.1.2 CDS預(yù)測

利用ANGEL軟件對獲得的基因片段進(jìn)行蛋白編碼區(qū)預(yù)測分析［24］，結(jié)果顯示一共有12 660個基因片段可視為蛋白編碼區(qū)，其序列長度的范圍為0～5 000 bp，主要集中于250～1 250 bp（圖1）。

圖1 CDS長度分布圖Fig.1 Statistics of sequence length of CDS

2.1.3 SSR預(yù)測

簡單重復(fù)序列標(biāo)記（simple sequence repeats，SSR），又稱為短串聯(lián)重復(fù)序列或微衛(wèi)星標(biāo)記。是一類由幾個核苷酸為重復(fù)單位組成的長達(dá)幾十個核苷酸的重復(fù)序列，長度較短，且廣泛均勻分布于真核生物基因組中。通過檢索Unigene序列，共發(fā)現(xiàn)14 277個SSR位點，分布在8 813個Unigene中。從表1可以發(fā)現(xiàn)，單核苷酸到三核苷酸重復(fù)最多，占總SSR位點數(shù)量的99.13%，其中單核苷酸重復(fù)以A／T重復(fù)基序最多；二核苷酸重復(fù)以AC／GT和AT／AT重復(fù)基序出現(xiàn)頻率最高；三核苷酸重復(fù)則以ATT／AAT和ATC／GAT為主要類型；四核苷酸重復(fù)則以GAAA／TTTC較多。其他五核苷酸和六核苷酸重復(fù)基元類型較多，數(shù)量非常少。

表1 錦繡龍蝦SSR不同重復(fù)基元分布情況Tab.1 Distribution of different repeat motifs in P.ornatus

2.1.4 lncRNA分析

長鏈非編碼RNA（lncRNA）是一類轉(zhuǎn)錄本長度超過200 nt、不編碼蛋白質(zhì)的RNA分子。由于建庫原理的限制，我們只能獲得含有polyA尾的lncRNA。使用CNCI［25］、CPC2［26］、Pfam［27］以及PLEK［28］對CD-HIT去冗余得到的基因進(jìn)行編碼潛能預(yù)測，經(jīng)數(shù)據(jù)庫比對后預(yù)測其編碼潛能，最終分析6 315個LncRNA序列，其中共有數(shù)目為1 342個（圖2）。

圖2 編碼潛能預(yù)測維恩圖Fig.2 Venn diagram of encoding potential prediction

2.2 Unigene的功能注釋

2.2.1 注釋結(jié)果統(tǒng)計

將錦繡龍蝦轉(zhuǎn)錄組所獲得的單基因簇序列在NR數(shù)據(jù)庫中進(jìn)行對比，共有11 086個Unigene獲得注釋，對比到424個物種，其中對比較多的前20個物種包括鉤蝦（Hyalellaazteca，3 907個）、濕木白蟻（Zootermopsisnevadensis，504個）、大型蚤（Daphniamagna，483個）、凡納濱對蝦（Litopenaeusvannamei，296個）、美洲鱟（Limulus polyphemus，266個）、斑節(jié)對蝦（Penaeusmonodon，257個）、克氏原螯蝦（Procambarusclarkii，196個）、中華絨螯蟹（Eriocheirsinensis，182個）、松葉蜂（Neodiprionlecontei，166個）、擬穴青蟹（Scylla paramamosain，160個）、淡水枝角水蚤（Daphnia pulex，139個）、日 本 囊 對 蝦（Marsupenaeus japonicus，133個）、鴨嘴舌形貝（Lingulaanatina，121個）、中國明對蝦（Fenneropenaeuschinensis，116個）、紅螯螯蝦（Cheraxquadricarinatus，115個）、美洲螯龍蝦（Homarusamericanus，114個）、白氏文昌魚（Branchiostomabelcheri，102個）、三疣梭子蟹（Portunustrituberculatus，85個）、羅氏沼蝦（Macrobrachium rosenbergii，83個）、淡 水 螯 蝦（Pacifastacusleniusculus，81個）。從NR數(shù)據(jù)庫的注釋結(jié)果來看，注釋到鉤蝦的基因最多，推測錦繡龍蝦與鉤蝦同源性較高。隨著錦繡龍蝦生物學(xué)研究的不斷深入，GenBank的基因數(shù)據(jù)越來越豐富，后期將獲得更多錦繡龍蝦基因注釋。

2.2.2 GO功能分類

GO是一套國際標(biāo)準(zhǔn)化的基因功能描述的分類系統(tǒng)，分為細(xì)胞組分（cellular component，CC）、分子功能（molecular function，MF）、生物學(xué)過程（biological process，BP）3大類［22］。將獲得的錦繡龍蝦Unigene與GO數(shù)據(jù)庫進(jìn)行匹配，分別有17 624個、10 398個和9 530個被劃分到BP、CC及MF 3大類，涵蓋上述功能類別的51個亞類（圖3），其中BP 24個亞類，涉及細(xì)胞進(jìn)程、代謝進(jìn)程及單生物進(jìn)程的較多，分別有3 768個、3 456個和2 918個；CC 18個亞類，涉及細(xì)胞、細(xì)胞部分及細(xì)胞器的較多，分別有1 894個、1 894個和1 428個；MF 9個亞類，涉及結(jié)合活性和催化活性的較多，分別有4 789個和3 398個。

圖3 Unigenes的GO功能分類圖Fig.3 GO functional categories of Unigenes

2.2.3 KOG分類統(tǒng)計

COG（cluster of orthologous groups of proteins）是根據(jù)細(xì)菌、藻類和真核生物完整基因組的編碼蛋白系統(tǒng)進(jìn)化關(guān)系分類構(gòu)建而成的蛋白數(shù)據(jù)庫，其中真核生物數(shù)據(jù)庫稱之為蛋白質(zhì)真核同源數(shù)據(jù)庫（eukaryotic orthologous groups，KOG）［19］。通過KOG數(shù)據(jù)庫對比到錦繡龍蝦單基因簇，共有9 433個獲得注釋信息，共分為26個功能組分（圖4），其中獲得注釋較多的功能組分一般功能預(yù)測（1 413個，14.98%）、其次是信號轉(zhuǎn)導(dǎo)機制（1 030個，10.92%）及翻譯后修飾、蛋白轉(zhuǎn)換、伴侶蛋白（917個，9.72%），其中未知蛋白僅有4個。

圖4 Unigenes的KOG功能分類Fig.4 KOG function classification of Unigenes

2.2.4 KEGG功能注釋

KEGG（Kyoto encyclopedia of genes and genomes）是系統(tǒng)分析基因產(chǎn)物和化合物在細(xì)胞中的代謝途徑以及基因產(chǎn)物功能的數(shù)據(jù)庫，已建立了一套完整KO注釋的系統(tǒng)，可完成新測序物種的基因組或轉(zhuǎn)錄組的功能注釋［29］。與KEGG數(shù)據(jù)庫對比，最終錦繡龍蝦單基因簇注釋到5大類34小類（圖5），其中基因數(shù)量較多的代謝通路有信號轉(zhuǎn)導(dǎo)機制通路（1 015個），轉(zhuǎn)運和分解代謝通路（726個），跨膜轉(zhuǎn)運的代謝通路基因最少，僅有9個。這些獲得注釋的單基因簇可為后續(xù)的錦繡龍蝦功能基因研究和應(yīng)用提供數(shù)據(jù)基礎(chǔ)。

2.2.5 轉(zhuǎn)錄因子分析

轉(zhuǎn)錄因子（transcription factor，TF）作為一類特殊的DNA結(jié)合蛋白，能與基因5′端上游的特定序列專一性結(jié)合，從而使目的基因能夠在特定的時空進(jìn)行表達(dá)，通過轉(zhuǎn)錄因子間以及與其他相關(guān)蛋白間的相互作用達(dá)到激活或抑制轉(zhuǎn)錄的效果，發(fā)揮著重要的調(diào)控作用。動物轉(zhuǎn)錄因子鑒定使用動物轉(zhuǎn)錄因子數(shù)據(jù)庫［30］animalTFDB 2.0，預(yù)測到轉(zhuǎn)錄因子有415個，隸屬于29個轉(zhuǎn)錄因子家族（圖6），其中比較多的轉(zhuǎn)錄因子家族：zf-C2H2家族有109個、ZBTB家族有75個，E2F家族、NFYB家族、THR-like家族、HSF家族、IRF家族、MBD家族和STAT家族最少，均只有2個，這些轉(zhuǎn)錄因子家族成員的獲取可為后期錦繡龍蝦生長發(fā)育、代謝調(diào)節(jié)、免疫應(yīng)答等相關(guān)研究奠定基礎(chǔ)。

圖6 錦繡龍蝦轉(zhuǎn)錄因子分析Fig.6 Transcription factor analysis of P.ornatus

3 討論

PacBio Sequel測序儀是美國太平洋生物技術(shù)公司（Pacific Biosciences）基于PacBio RSII平臺最新推出的第三代測序平臺，該系統(tǒng)以其測序通量高、單Gb數(shù)據(jù)成本低、周期短等特點廣受青睞。PacBio第三代測序技術(shù)的最大特點是無需進(jìn)行PCR擴增，可直接讀取目標(biāo)序列，是轉(zhuǎn)錄組從頭測序的首選［31］。對于無基因組參考的水生生物來說，運用三代測序技術(shù)進(jìn)行轉(zhuǎn)錄組分析以及功能基因的挖掘是較好的策略［32－34］。本研究對錦繡龍蝦全長轉(zhuǎn)錄組測序及分析，共獲得39 828 440個子序列，平均子序列長度為1 786 bp，N50為2 492 bp，結(jié)果顯示轉(zhuǎn)錄組所獲得的組裝序列完整性較好。利用NR、Nt、KOG、SwissProt等7大公共數(shù)據(jù)庫進(jìn)行功能注釋分類共獲得序列28 034個，有11 086個獲得注釋，對比到424個物種，有10 463個注釋到350個代謝通路，有9 433個獲得KOG注釋，分為26個功能組分。曾地剛等［35］對凡納濱對蝦肝胰腺采用二代轉(zhuǎn)錄組測序，結(jié)果 獲得500 177個EST（expressed sequence tags），序列平均長度363 bp。拼接獲得20 225個Unigene，序列平均長度507 bp。注釋到NR、COG以及KEGG數(shù)據(jù)庫的分別為13 676個、4 645個、4 104個。李喜蓮等［36］以紅螯螯蝦肝臟、精巢以及卵巢為材料進(jìn)行二代轉(zhuǎn)錄組測序，獲得了6 736個單基因簇，注釋到GO的為16 989個，注釋到COG的為4 697個，注釋到KEGG的為9 842個。對比研究結(jié)果顯示，三代測序獲得的序列質(zhì)量、基因數(shù)以及注釋到的基因信息均優(yōu)于二代測序。

基于PacBio測序平臺對構(gòu)建甲殼動物轉(zhuǎn)錄組文庫序列聚類和改良，ZHANG等［37］通過凡納濱對蝦全長轉(zhuǎn)錄組文庫獲得72 648條高質(zhì)量序列，WANG等［38］通過中國明對蝦進(jìn)全長轉(zhuǎn)錄組文庫獲得10 795條高質(zhì)量序列，POOTAKHAM等［39］通過斑節(jié)對蝦全長轉(zhuǎn)錄組文庫得到22 418條高質(zhì)量序列。本研究通過錦繡龍蝦全長轉(zhuǎn)錄組文庫最終獲得了13 297條高質(zhì)量序列，相較于凡納濱對蝦和斑節(jié)對蝦少，除了物種的差異外，也歸因于轉(zhuǎn)錄組聚類過程中的參數(shù)設(shè)置。同時，測序過程中5′末端的降解也可能導(dǎo)致同一轉(zhuǎn)錄組序列進(jìn)行不同的聚類劃分?？紤]到這一點，本研究將錦繡龍蝦全長轉(zhuǎn)錄組序列聚類，獲得高質(zhì)量的一致序列用于后續(xù)分析。

本研究利用Nr、KEGG等數(shù)據(jù)庫對錦繡龍蝦轉(zhuǎn)錄組序列進(jìn)行功能注釋，結(jié)果顯示有28 034個單基因簇獲得注釋信息，部分基因序列仍沒有獲取注釋，可能由于錦繡龍蝦是未測序物種，公共數(shù)據(jù)庫中缺乏相關(guān)的功能注釋信息，也可能是部分錦繡龍蝦單基因簇序列較短造成。將獲得的蛋白同源序列與NR數(shù)據(jù)庫對比發(fā)現(xiàn)，與鉤蝦對比的同源信息最多，占35%，推測可能是由于鉤蝦與錦繡龍蝦的進(jìn)化史和生活史較為相似。將錦繡龍蝦轉(zhuǎn)錄組單基因簇注釋到GO數(shù)據(jù)庫，可劃分為3大類共計51個分支，與KOG數(shù)據(jù)庫對比可分為26類，GO和KOG注釋功能分類的結(jié)果顯示錦繡龍蝦單基因簇的功能涉及了各類生命活動。KEGG注釋到5大類34小類，其中涉及代謝通路和次生物質(zhì)的生物合成基因較多。從轉(zhuǎn)錄組共檢測出14 277個SSR位點，這些SSR在提高物種遺傳多樣性潛能方面發(fā)揮著重要的作用，為下一階段開發(fā)錦繡龍蝦多態(tài)性SSR分子標(biāo)記提供了基礎(chǔ)數(shù)據(jù)。通過對錦繡龍蝦轉(zhuǎn)錄組進(jìn)行注釋分析，初步闡明了錦繡龍蝦基因的功能、參與的生物過程、所在的代謝途徑或信號通路等，對于后期深入研究關(guān)鍵基因的功能提供了信息，為發(fā)掘錦繡龍蝦功能基因、研究相關(guān)生理功能奠定了基礎(chǔ)。