波紋龍蝦全長轉(zhuǎn)錄組測序分析

梁妃爽，梁華芳，孫榕澤，徐思行，溫崇慶，王 偉 (廣東海洋大學(xué)水產(chǎn)學(xué)院,廣東湛江 524088)

轉(zhuǎn)錄組(Transcriptome)是指特定細(xì)胞或組織中的所有轉(zhuǎn)錄產(chǎn)物，包含信使 RNA、核糖體 RNA、轉(zhuǎn)運 RNA以及非編碼 RNA。1977年，第1代DNA測序技術(shù)(Sanger雙脫氧鏈終止法)開始探索基因結(jié)構(gòu)，但 Sanger 測序方法速度慢、通量低、成本高，難以滿足大量測序要求，因此難以應(yīng)用于組學(xué)測序研究。第2代測序技術(shù)主要基于Roche/454、ABI/Solid、Illumina/Solexa測序平臺，有效解決了速度慢、成本高、通量低的問題，但是其讀長短，拼接得到的轉(zhuǎn)錄本結(jié)構(gòu)不完整。PacBio平臺單分子實時測序技術(shù)，也稱為SMRT(Single Molecule Real-Time)測序，因其超長的測序讀長、超高的測序通量、無GC偏好性、無PCR擴增偏向性、直接檢測堿基修飾等諸多優(yōu)勢而廣泛應(yīng)用于水產(chǎn)領(lǐng)域研究，為基因組學(xué)、轉(zhuǎn)錄組學(xué)及DNA甲基化等研究注入了新活力。張金勇等對金烏賊()采用全長轉(zhuǎn)錄組測序，篩選SSR位點并分析出現(xiàn)頻率較高且類型豐富、多態(tài)性潛能較高等特點。Zhang等利用三代測序技術(shù)篩選施氏鱘 () 性腺中早期配子發(fā)生的相關(guān)基因，探討了其生殖調(diào)控機制。Pootakham等利用三代測序技術(shù)對斑節(jié)對蝦生成第1個全長轉(zhuǎn)錄組。因此，利用三代測序技術(shù)對波紋龍蝦全長轉(zhuǎn)錄組測序，分析環(huán)境因子相關(guān)脅迫差異基因，以期為相關(guān)功能基因的研究以及波紋龍蝦養(yǎng)殖過程中水質(zhì)調(diào)控和抗環(huán)境因子變化脅迫等提供理論支持和數(shù)據(jù)支撐。

波紋龍蝦隸屬于十足目(Decapoda)龍蝦科(Palinuroidea)龍蝦屬()，是我國沿海地區(qū)的重要經(jīng)濟種類，營養(yǎng)豐富，味道鮮美。目前波紋龍蝦基因組尚未測序完成，其基因序列信息比較匱乏，相關(guān)的分子遺傳基礎(chǔ)也很薄弱，分子生物學(xué)水平研究也較少，僅見李斌等研究了波紋龍蝦C-型凝集素PhLecA的基因克隆與表達(dá)，Zhuo等研究了波紋龍蝦蛻皮激素受體的分子克隆、特征和表達(dá)分析，羅嘉俊等研究了波紋龍蝦GIH基因克隆、表達(dá)及其對光周期的響應(yīng)等。筆者采用三代測序技術(shù)的PacBio 單分子實時測序平臺對波紋龍蝦()進(jìn)行全長轉(zhuǎn)錄組測序，通過生物信息學(xué)方法進(jìn)行序列拼接、功能注釋、分類和代謝通路分析，獲取豐富的波紋龍蝦序列信息，旨在為進(jìn)一步挖掘波紋龍蝦相關(guān)功能基因、基因組學(xué)及開發(fā)分子標(biāo)記等研究奠定基礎(chǔ)。

1 材料與方法

波紋龍蝦為購自海南省瓊海市青葛村，體質(zhì)量為21～57 g的幼龍蝦，采用干運法運回廣東海洋大學(xué)東海島生物研究基地，在20 m的水泥池中暫養(yǎng)，用黑色幕布遮光，24 h不間斷充氣。試驗海水來源于自然海區(qū)，鹽度為26‰～30‰，pH 為8.1～8.3，每天換水50%，投喂菲律賓蛤仔()等主要鮮活餌料，養(yǎng)殖21 d后開始試驗。取正常條件養(yǎng)殖下波紋龍蝦的肝胰腺、鰓、肌肉、性腺、眼柄和80 mg/L亞硝酸鹽(NaNO)脅迫7 d的肝胰腺(做3+2項目，取亞硝酸鹽脅迫混合作全長轉(zhuǎn)錄組，主要是為了進(jìn)行二代測序時方便篩選差異基因用于后續(xù)試驗)，迅速放入RNAlater中保存，后用干冰保存寄往北京諾禾致源科技股份有限公司進(jìn)行RNA提取和相關(guān)測序。

利用Trizol法分別對建庫及定量試驗所用的波紋龍蝦組織進(jìn)行RNA提取，使用1%瓊脂糖凝膠電泳分析RNA純度和完整性，Agilent 2100軟件準(zhǔn)確檢測RNA完整性，用Qubit 2.0軟件精確定量RNA濃度，Nanodrop軟件檢測RNA純度，并選擇符合測序標(biāo)準(zhǔn)的RNA等量混合用于文庫建設(shè)?；旌虾蟮腞NA樣品純度和完整性：6例混樣的Qubit濃度為112 μL,Qubit體積為39 μL,Qubit總量為4.368 μg,樣品純度1.949,1.303,Nano濃度為236.417 ng/μL,RIN值為3.8,NC/QC為2.11。結(jié)果表明，該例樣本基線平整，符合三代建庫標(biāo)準(zhǔn)。波紋龍蝦是水產(chǎn)低等動物，RNA只有單峰值，檢測時會出現(xiàn)其他峰值干擾從而導(dǎo)致RIN值較低的情況，但6例混樣的組織RIN值均符合二代測序要求。

文庫的構(gòu)建流程如圖1，構(gòu)建好的文庫進(jìn)行質(zhì)量檢測后，在Illumina高通量測序平臺NovaSeq 6000進(jìn)行測序(諾禾致源科技有限公司，天津)。

圖1 轉(zhuǎn)錄組文庫構(gòu)建Fig.1 Transcriptome library construction

測序完成后對原始下機數(shù)據(jù)進(jìn)行去接頭和低質(zhì)量reads，得到高質(zhì)量數(shù)據(jù)，采用軟件SMRTlink v8.0對高質(zhì)量數(shù)據(jù)進(jìn)行序列組裝，得到波紋龍蝦的轉(zhuǎn)錄組基因數(shù)據(jù)(Unigene)庫。序列分析步驟分別是使用下機數(shù)據(jù)中subreads.bam文件通過CCS算法，對單分子多測序序列進(jìn)行自我更正，獲得CCS (circular consensus sequence)序列；通過檢測CCS是否包含poly-A、5′-primer、3′-primer,對CCS進(jìn)行分類并找出FLNC(full-length non chimera：全長非嵌合序列)序列和nFL(Non-Full-Length:非全長非嵌合序列)序列；將同一轉(zhuǎn)錄本的FLNC序列使用hierarchical n*log(n)算法聚類，得到consensus序列；最后對由此產(chǎn)生的全長序列進(jìn)行polish，獲得Polished consensus序列進(jìn)行后續(xù)分析。具體流程見圖2。

全長轉(zhuǎn)錄組序列的功能注釋及結(jié)構(gòu)分析去冗余后的序列使用 CD-HIT 軟件進(jìn)行基因注釋，使用的數(shù)據(jù)庫包括：NR， KOG/COG，NT，Pfam，KEGG，Swiss-Prot，GO。

2 結(jié)果與分析

轉(zhuǎn)錄組測序數(shù)據(jù)質(zhì)量分析。波紋龍蝦正常條件下的肝胰腺、鰓絲、眼柄、肌肉和性腺，以及亞硝酸鹽脅迫下的肝胰腺的轉(zhuǎn)錄組測序數(shù)據(jù)見表1。堿基質(zhì)量及組成分析顯示，GC 含量區(qū)間為35.34%～48.26%，各組織樣品Q20堿基百分比不小于97.56%，Q30 堿基百分比不小于92.92%。這說明測序產(chǎn)出質(zhì)量符合要求，能用于后續(xù)組裝分析。轉(zhuǎn)錄本校正分析得出平均序列長度4 147 bp,N50為4 671 bp,注釋率為93.74%。

圖2 Iso-Seq分析流程Fig.2 Iso-Seq analysis flow chart

表1 數(shù)據(jù)產(chǎn)出質(zhì)量情況

CDS預(yù)測。CDS(coding sequence)是編碼一段蛋白產(chǎn)物的序列。在全長轉(zhuǎn)錄組的測序結(jié)果中，預(yù)測蛋白質(zhì)編碼區(qū)有助于基因的初步分析，同時也是進(jìn)行后續(xù)蛋白結(jié)構(gòu)分析的基礎(chǔ)。利用ANGEL軟件進(jìn)行CDS預(yù)測分析，結(jié)果顯示共有1 043個基因片段可視為蛋白編碼區(qū)，其序列長度為0～7 500 bp，主要集中于300～2 500 bp(圖3)。

圖3 CDS長度分布Fig.3 The statistics of sequence length of CDS

lncRNA分析。LncRNA是一類轉(zhuǎn)錄本長度超過200 nt，不編碼蛋白質(zhì)的RNA分子。由于建庫原理的限制，只能獲得含有polyA尾的lncRNA。使用CNCI、PLEK、CPC2軟件以及Pfam數(shù)據(jù)庫對PacBio測序數(shù)據(jù)進(jìn)行編碼潛能預(yù)測，最終分析得到3 105個LncRNA序列，其中共有數(shù)目為272個(圖4)。

圖4 編碼潛能預(yù)測維恩圖Fig.4 Encoding potential prediction Venn diagram

轉(zhuǎn)錄本分析。測序結(jié)果與數(shù)據(jù)組裝使用PacBio 測序平臺對波紋龍蝦鰓、肝胰腺、肌肉、性腺和眼柄等組織混樣進(jìn)行全長轉(zhuǎn)錄組測序，對原始數(shù)據(jù)進(jìn)行過濾，共獲得17 044 319個子序列(大小59.47 Gb)，平均子序列長度為3 490 bp，N50為4 037 bp。通過每個ZMW孔中子序列的CCS聚類之后得到的序列數(shù)為517 682個，序列平均長度為4 181 bp，N50為4 965 bp。同時含有3′引物和5′引物，以及3′引物前含有polyA尾的全長序列(Full-Length，F(xiàn)L)459 737個，全長非嵌合序列(Full-Length non-chimericRead，F(xiàn)LNC)458 653個，序列平均長度為4 094 bp，N50為4 919 bp ，F(xiàn)LNC/CCS為88.60%。全長轉(zhuǎn)錄組得到改良后一致序列21 524個，10 425個Unigenes(圖5)，序列平均長度為4 008 bp，N50為4 474 bp。

圖5 波紋龍蝦轉(zhuǎn)錄本的長度分布Fig.5 Length distribution of Panulirus homarus Unigene

Unigene的功能注釋NR數(shù)據(jù)庫注釋到Unigene的數(shù)量最多，為9 580個，NT數(shù)據(jù)庫注釋到的最少，僅3 498個(圖6)。

圖6 七大數(shù)據(jù)庫注釋統(tǒng)計結(jié)果Fig.6 Annotation statistical results of seven databases

NR分析。NR數(shù)據(jù)庫注釋將波紋龍蝦轉(zhuǎn)錄組所獲得的單基因簇序列在NR數(shù)據(jù)庫中比對，共比對到358個物種，其中鉤蝦()的同源序列最多，為3 865 個，占注釋序列總數(shù)的 40.34%，推測波紋龍蝦與鉤蝦同源性較高；其次為濕木白蟻()563個，美洲鱟()291個，大型蚤()282個，凡納濱對蝦 ()171個，斑節(jié)對蝦()162個，鴨嘴舌形貝()142個，淡水枝角水蚤()134個，中華絨螯蟹()128個，葉蟬()103個，日本囊對蝦()96個，白氏文昌魚()92個，囊舌蟲()91個，赤擬谷盜()83個，鞘翅鳥()82個，溫室希蛛()74個，裸長角蟲兆()72個，淡水螯蝦()70個，紅螯螯蝦()68個，克氏原螯蝦()66個，其他物種2 945個。

GO 功能注釋。GO 功能注釋結(jié)果見圖7～9，共有7 585條Unigenes被注釋分類，從細(xì)胞組分可細(xì)分為16 類，占比最多的是細(xì)胞(cell)和細(xì)胞組成(cell part)(44.85%)，其次為細(xì)胞器(organelle)(17.57%)(圖7)。從分子功能細(xì)分為10類，其中捆綁(binding)包含Unigenes最多(64.31%)，催化活性(catalytic activity)次之(41.08%)(圖8)。生物學(xué)過程可細(xì)分為24 類，細(xì)胞過程(cellular process)包含Unigenes最多(42.50%)，代謝過程(Metabolic process)類次之(40.07%)(圖9)。

圖7 細(xì)胞組分Fig.7 Cellular component

圖8 分子功能Fig.8 Molecular function

圖9 生物學(xué)過程Fig.9 Biological process

KOG 功能分類。KOG功能注釋結(jié)果共有7 436條 Unigenes 被注釋分類，分布于26 類(圖10)。其中只是一般功能預(yù)測類共有1 421條注釋信息，占比最大 (19.11%)，其次為信號轉(zhuǎn)導(dǎo)機制，有1 238條注釋信息(16.65%)，未知蛋白僅有7條(7.13%)。

KEGG功能注釋分析。KEGG功能注釋結(jié)果顯示，共有9 254條 Unigenes 被注釋分類，分布于 347 個已知途徑中，其中前 12 個代謝途徑，注釋基因數(shù)占總量的 25.26%。前 5 個途徑分別是心肌細(xì)胞的腎上腺素能信號(ko04261)289 條、病毒性心肌炎(ko05416)254 條、心臟肌肉收縮(ko04260)236 條、癌癥中的蛋白多糖(ko05205)200 條和局部黏連(ko04510)191 條(表2)。

表2 前12個代謝途徑基因數(shù)量Table 2 Number of genes in the first 12 metabolic pathways

轉(zhuǎn)錄因子分析。轉(zhuǎn)錄因子(transcription factor,TF)作為一類特殊的DNA結(jié)合蛋白，可與基因5′末端上游的特定序列結(jié)合，使目的基因可以特定時空表達(dá)，通過轉(zhuǎn)錄因子與其他相關(guān)蛋白質(zhì)的相互作用來激活或抑制轉(zhuǎn)錄效果，發(fā)揮著重要的調(diào)控作用。動物轉(zhuǎn)錄因子鑒定使用動物轉(zhuǎn)錄因子數(shù)據(jù)庫—animal TFDB 2.0預(yù)測到轉(zhuǎn)錄因子家族共有543個，屬于29個家庭(圖11)，其中轉(zhuǎn)錄因子家族較多的有：zf-C2H2家族有190個、ZBTB家族有110個，TEA家族最少，只有1個，這些轉(zhuǎn)錄因子家族成員的獲取可為后期波紋龍蝦生長發(fā)育、代謝調(diào)節(jié)、免疫應(yīng)答等相關(guān)研究奠定基礎(chǔ)。

3 討論

轉(zhuǎn)錄組測序技術(shù)是一種成本低，能快速獲取大量轉(zhuǎn)錄數(shù)據(jù)并對研究生物體生物學(xué)特性、基因功能、相關(guān)代謝途徑和信號通路等具有重要作用的測序技術(shù)。轉(zhuǎn)錄組測序分析技術(shù)廣泛應(yīng)用于水產(chǎn)養(yǎng)殖相關(guān)研究中，并已成為研究環(huán)境脅迫對甲殼動物免疫、生長、繁殖以及蛻殼等過程的影響的重要手段之一。

圖10 KOG數(shù)據(jù)庫注釋統(tǒng)計Fig.10 KOG database annotation statistics

圖11 轉(zhuǎn)錄因子分析Fig.11 Transcription factor analysis

該研究對波紋龍蝦全長轉(zhuǎn)錄組進(jìn)行測序及分析，共獲得21 524個轉(zhuǎn)錄本，10 425個Unigenes，轉(zhuǎn)錄本校正分析得出平均序列長度為4 147 bp，N50為4 671 bp，注釋率為93.74%，測序結(jié)果可以看出，組裝得到序列完整性較好。利用NR、NT、KOG、KEGG等七大公共數(shù)據(jù)庫進(jìn)行功能注釋分類，有9 580 個獲得NR數(shù)據(jù)庫注釋，對比到358個物種，與鉤蝦對比的同源信息最多，占40.34%，推測可能是由于鉤蝦與波紋龍蝦的進(jìn)化史和繁殖習(xí)性較為相似。將獲得的波紋龍蝦單基因簇與GO數(shù)據(jù)庫進(jìn)行匹配，有7 585個得到GO注釋，被劃分到BP、CC及MF3大類，涵蓋上述功能類別的50個亞類。通過KOG數(shù)據(jù)庫對比波紋龍蝦單基因簇，共有7 436個獲得注釋信息，共分為26個功能組分。與KEGG數(shù)據(jù)庫對比，最終波紋龍蝦單基因簇注釋到6大類43小類，其中基因數(shù)量較多的代謝通路有信號轉(zhuǎn)導(dǎo)機制通路887個。李喜蓮等進(jìn)行紅螯螯蝦肝臟、卵巢和精巢二代測序獲得了6 736條Unigene，注釋到GO為16 989個，注釋到COG為4 697個，注釋到KEGG為9 842個。陳雪峰等對羅氏沼蝦卵巢4個不同發(fā)育期進(jìn)行2代測序產(chǎn)生了95 379個Unigenes，注釋到GO為6 422個，注釋到KEGG為8 423個。沈曄等進(jìn)行脊尾白蝦對低鹽脅迫響應(yīng)的轉(zhuǎn)錄組學(xué)分析，結(jié)果獲得了72 734 條Unigenes，注釋到NR為21 931個。由此看出，單分子實時測序技術(shù)獲得的序列質(zhì)量、基因數(shù)量和注釋基因信息優(yōu)于第二代測序。

李斌等用波紋龍蝦肝胰腺和卵巢組織的mRNA表達(dá)譜進(jìn)行了2代轉(zhuǎn)錄組測序，測序結(jié)果與該研究3代測序結(jié)果相差較大，2代測序總Unigene 74 124個，而該研究為10 425個，造成3代測序的Unigene比2代少的原因是3代測序進(jìn)行了無參轉(zhuǎn)錄組對照測序，只能比對到其他數(shù)據(jù)庫，這樣測序出來的結(jié)果就遠(yuǎn)比2代的少。2代測序基因功能注釋率為33.80%，該研究為93.74%；2代測序僅14.00%注釋到GO數(shù)據(jù)庫，22.70%注釋到KEGG代謝途徑，12.00%注釋到COG蛋白數(shù)據(jù)庫。造成這種結(jié)果的原因可能是對波紋龍蝦開展研究少，國內(nèi)外對波紋龍蝦的研究報道也相對較少，在NCBI數(shù)據(jù)庫中找到龍蝦屬的核酸序列不足1 000條；也有可能是個體間的差異較大和測序上樣量的不同，2代測序和3代測序條件和方法有所差異。生物分子數(shù)據(jù)庫的完善對波紋龍蝦的研究、養(yǎng)殖和保護(hù)起著重要作用，因此加大波紋龍蝦的分子生物學(xué)研究力度至關(guān)重要。

Zhang等通過凡納濱對蝦全長轉(zhuǎn)錄組文庫獲得72 648條高質(zhì)量序列，Wang等通過對中國對蝦進(jìn)行全長轉(zhuǎn)錄組文庫測序獲得10 795條高質(zhì)量序列。該研究通過波紋龍蝦全長轉(zhuǎn)錄組文庫最終獲得10 425條高質(zhì)量序列，較斑節(jié)對蝦、凡納濱對蝦和中國對蝦少，這可能是物種之間的差異。這些數(shù)據(jù)為進(jìn)一步了解波紋龍蝦的生物學(xué)特性、基因功能、相關(guān)代謝途徑和信號通路等提供理論基礎(chǔ)，為后續(xù)研究提供一定參考。