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

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

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

1 材料與方法

波紋龍蝦為購自海南省瓊海市青葛村，體質量為21～57 g的幼龍蝦，采用干運法運回廣東海洋大學東海島生物研究基地，在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項目，取亞硝酸鹽脅迫混合作全長轉錄組，主要是為了進行二代測序時方便篩選差異基因用于后續(xù)試驗)，迅速放入RNAlater中保存，后用干冰保存寄往北京諾禾致源科技股份有限公司進行RNA提取和相關測序。

利用Trizol法分別對建庫及定量試驗所用的波紋龍蝦組織進行RNA提取，使用1%瓊脂糖凝膠電泳分析RNA純度和完整性，Agilent 2100軟件準確檢測RNA完整性，用Qubit 2.0軟件精確定量RNA濃度，Nanodrop軟件檢測RNA純度，并選擇符合測序標準的RNA等量混合用于文庫建設?；旌虾蟮腞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。結果表明，該例樣本基線平整，符合三代建庫標準。波紋龍蝦是水產(chǎn)低等動物，RNA只有單峰值，檢測時會出現(xiàn)其他峰值干擾從而導致RIN值較低的情況，但6例混樣的組織RIN值均符合二代測序要求。

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

圖1 轉錄組文庫構建Fig.1 Transcriptome library construction

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

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

2 結果與分析

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

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

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

CDS預測。CDS(coding sequence)是編碼一段蛋白產(chǎn)物的序列。在全長轉錄組的測序結果中，預測蛋白質編碼區(qū)有助于基因的初步分析，同時也是進行后續(xù)蛋白結構分析的基礎。利用ANGEL軟件進行CDS預測分析，結果顯示共有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是一類轉錄本長度超過200 nt，不編碼蛋白質的RNA分子。由于建庫原理的限制，只能獲得含有polyA尾的lncRNA。使用CNCI、PLEK、CPC2軟件以及Pfam數(shù)據(jù)庫對PacBio測序數(shù)據(jù)進行編碼潛能預測，最終分析得到3 105個LncRNA序列，其中共有數(shù)目為272個(圖4)。

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

轉錄本分析。測序結果與數(shù)據(jù)組裝使用PacBio 測序平臺對波紋龍蝦鰓、肝胰腺、肌肉、性腺和眼柄等組織混樣進行全長轉錄組測序，對原始數(shù)據(jù)進行過濾，共獲得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%。全長轉錄組得到改良后一致序列21 524個，10 425個Unigenes(圖5)，序列平均長度為4 008 bp，N50為4 474 bp。

圖5 波紋龍蝦轉錄本的長度分布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)計結果Fig.6 Annotation statistical results of seven databases

NR分析。NR數(shù)據(jù)庫注釋將波紋龍蝦轉錄組所獲得的單基因簇序列在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 功能注釋結果見圖7～9，共有7 585條Unigenes被注釋分類，從細胞組分可細分為16 類，占比最多的是細胞(cell)和細胞組成(cell part)(44.85%)，其次為細胞器(organelle)(17.57%)(圖7)。從分子功能細分為10類，其中捆綁(binding)包含Unigenes最多(64.31%)，催化活性(catalytic activity)次之(41.08%)(圖8)。生物學過程可細分為24 類，細胞過程(cellular process)包含Unigenes最多(42.50%)，代謝過程(Metabolic process)類次之(40.07%)(圖9)。

圖7 細胞組分Fig.7 Cellular component

圖8 分子功能Fig.8 Molecular function

圖9 生物學過程Fig.9 Biological process

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

KEGG功能注釋分析。KEGG功能注釋結果顯示，共有9 254條 Unigenes 被注釋分類，分布于 347 個已知途徑中，其中前 12 個代謝途徑，注釋基因數(shù)占總量的 25.26%。前 5 個途徑分別是心肌細胞的腎上腺素能信號(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

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

3 討論

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

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

圖11 轉錄因子分析Fig.11 Transcription factor analysis

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

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

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