波紋龍蝦全長轉錄組測序分析

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

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

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

1 材料與方法

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

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

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

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

測序完成后對原始下機數據進行去接頭和低質量reads，得到高質量數據，采用軟件SMRTlink v8.0對高質量數據進行序列組裝，得到波紋龍蝦的轉錄組基因數據(Unigene)庫。序列分析步驟分別是使用下機數據中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序列；最后對由此產生的全長序列進行polish，獲得Polished consensus序列進行后續分析。具體流程見圖2。

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

2 結果與分析

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

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

表1 數據產出質量情況

CDS預測。CDS(coding sequence)是編碼一段蛋白產物的序列。在全長轉錄組的測序結果中，預測蛋白質編碼區有助于基因的初步分析，同時也是進行后續蛋白結構分析的基礎。利用ANGEL軟件進行CDS預測分析，結果顯示共有1 043個基因片段可視為蛋白編碼區，其序列長度為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數據庫對PacBio測序數據進行編碼潛能預測，最終分析得到3 105個LncRNA序列，其中共有數目為272個(圖4)。

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

轉錄本分析。測序結果與數據組裝使用PacBio 測序平臺對波紋龍蝦鰓、肝胰腺、肌肉、性腺和眼柄等組織混樣進行全長轉錄組測序，對原始數據進行過濾，共獲得17 044 319個子序列(大小59.47 Gb)，平均子序列長度為3 490 bp，N50為4 037 bp。通過每個ZMW孔中子序列的CCS聚類之后得到的序列數為517 682個，序列平均長度為4 181 bp，N50為4 965 bp。同時含有3′引物和5′引物，以及3′引物前含有polyA尾的全長序列(Full-Length，FL)459 737個，全長非嵌合序列(Full-Length non-chimericRead，FLNC)458 653個，序列平均長度為4 094 bp，N50為4 919 bp ，FLNC/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數據庫注釋到Unigene的數量最多，為9 580個，NT數據庫注釋到的最少，僅3 498個(圖6)。

圖6 七大數據庫注釋統計結果Fig.6 Annotation statistical results of seven databases

NR分析。NR數據庫注釋將波紋龍蝦轉錄組所獲得的單基因簇序列在NR數據庫中比對，共比對到358個物種，其中鉤蝦()的同源序列最多，為3 865 個，占注釋序列總數的 40.34%，推測波紋龍蝦與鉤蝦同源性較高；其次為濕木白蟻()563個，美洲鱟()291個，大型蚤()282個，凡納濱對蝦 ()171個，斑節對蝦()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 個代謝途徑，注釋基因數占總量的 25.26%。前 5 個途徑分別是心肌細胞的腎上腺素能信號(ko04261)289 條、病毒性心肌炎(ko05416)254 條、心臟肌肉收縮(ko04260)236 條、癌癥中的蛋白多糖(ko05205)200 條和局部黏連(ko04510)191 條(表2)。

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

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

3 討論

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

圖10 KOG數據庫注釋統計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等七大公共數據庫進行功能注釋分類，有9 580 個獲得NR數據庫注釋，對比到358個物種，與鉤蝦對比的同源信息最多，占40.34%，推測可能是由于鉤蝦與波紋龍蝦的進化史和繁殖習性較為相似。將獲得的波紋龍蝦單基因簇與GO數據庫進行匹配，有7 585個得到GO注釋，被劃分到BP、CC及MF3大類，涵蓋上述功能類別的50個亞類。通過KOG數據庫對比波紋龍蝦單基因簇，共有7 436個獲得注釋信息，共分為26個功能組分。與KEGG數據庫對比，最終波紋龍蝦單基因簇注釋到6大類43小類，其中基因數量較多的代謝通路有信號轉導機制通路887個。李喜蓮等進行紅螯螯蝦肝臟、卵巢和精巢二代測序獲得了6 736條Unigene，注釋到GO為16 989個，注釋到COG為4 697個，注釋到KEGG為9 842個。陳雪峰等對羅氏沼蝦卵巢4個不同發育期進行2代測序產生了95 379個Unigenes，注釋到GO為6 422個，注釋到KEGG為8 423個。沈曄等進行脊尾白蝦對低鹽脅迫響應的轉錄組學分析，結果獲得了72 734 條Unigenes，注釋到NR為21 931個。由此看出，單分子實時測序技術獲得的序列質量、基因數量和注釋基因信息優于第二代測序。

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

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