基于PacBio 平臺的黃顙魚全長轉錄組測序及分析

王家琪，熊陽，韓慶慶，皇培培，梅潔

（華中農業大學水產學院，武漢 430070）

近年來，轉錄組學技術廣泛應用于水產動物繁育、營養、發育和免疫等各研究［1］。目前，轉錄組測序應用最廣的是二代測序技術（RNA-Seq），二代轉錄組測序具有測序通量高、成本低的優勢。盡管二代測序技術讀取準確率高，但讀長相對較短，給后續序列組裝、拼接以及注釋等帶來困難［2］。而基于PacBio 平臺的單分子實時測序技術（Single molecule real time，SMRT）的第三代全長轉錄組測序技術，在測序過程不需要打斷RNA 片段，可超長讀取包含單個完整轉錄本序列信息、后續無需序列拼接與組裝，極大地提高基于功能注釋的準確性［3］。此外，三代全長轉錄組測序可直接讀取目標序列，無須PCR 擴增等步驟，大大降低了假陽性率，同時也避免偏置及堿基替換等問題，精準度高達99.9%［3，4］。

黃顙魚（Pelteobagrus fulvidraco）是中國一種重要的小型淡水經濟魚類，肉質鮮美且營養價值高。經全國水產原種和良種審定委員會審定的黃顙魚新品種有全雄黃顙魚“全雄1 號”（GS-04-001-2010）和雜交黃顙魚“黃優1 號”（GS-02-001-2018）［5，6］。遺傳育種技術的創新和新品種的培育推動了黃顙魚產業的發展，2019 年黃顙魚年產量高達53.69 萬t［7］。隨著研究技術的不斷革新，關于黃顙魚分子遺傳育種的研究逐漸增多。從傳統的基因克隆、單核苷酸多態性（SNP）［8］、簡單重復序列（SSR）［9］和擴增片段長度多態性（AFLP）［10］等分析，到轉錄組［11］、基因組［12］和蛋白質組［13］等組學大數據分析，為黃顙魚種質資源挖掘和遺傳育種奠定了基礎。關于黃顙魚不同組織的二代轉錄組分析在近些年也有不少報道，馮美惠等［14］對飼料中添加維生素D3的黃顙魚腸和腎臟組織進行轉錄組分析；Wu等［15］對XY 黃顙魚和YY 超雄黃顙魚的精巢組織進行轉錄比較分析；Chen等［11］對黃顙魚的卵巢、精巢、肝臟、腎臟、肌肉、腦、脾和心臟8 種組織進行454 焦磷酸測序法混樣測序分析。目前在黃顙魚中還沒有第三代全長轉錄組測序結果的報道。

本研究采用PacBio 平臺的第三代測序技術對黃顙魚10 種組織的RNA 混樣進行全長轉錄組測序，測序結果與黃顙魚基因組進行對比分析，挖掘出新基因和已知基因的同源異構體，進行序列分析、功能注釋和基因結構分析，為黃顙魚分子遺傳育種提供科學理論。

1 材料與方法

1.1 組織樣本采集

三齡性成熟的黃顙魚購買于武漢百瑞生物有限公司。使用MS-222 將黃顙魚麻醉后，分離出肝臟、腎、背部肌肉、腦、脾、心臟、皮膚、血液、鰓、性腺等組織，用液氮速凍后送至武漢菲沙基因信息有限公司進行RNA 提取、質量和濃度測定及測序分析。

1.2 文庫構建

利用Trizol 法提取黃顙魚組織的總RNA。通過Nanodrop 檢測RNA 的純度（OD260nm/280nm）和濃度，Agilent 2100 對RNA 的完整性進行精確檢測；使用瓊脂糖凝膠電泳檢測有無基因組DNA 污染。以上各組織RNA 檢測合格后進行等量混勻，使用SMARTer PCR cDNA Synthesis Kit 合成全長cDNA，全長cDNA片段通過BluePippin 篩選共獲得3 個文庫（1～2 kb，2～6 kb（a）和2～6 kb（b））；通過PCR技術對全長cDNA進行擴增；對全長cDNA 進行末端修復，加上SMRT啞鈴型接頭和使用核酸外切酶消化；通過BluePippin 進行二次篩選，獲得測序文庫。使用Qubit 2.0 和Agilent 2100 對構建的文庫進行質量檢測，檢測結果達到要求后進行上機測序。

1.3 全長轉錄組測序數據分析

使用PacBio 平臺對檢測合格的文庫進行測序。對測序下機原始輸出數據使用SMRT Link v5.0 進行處理，獲得Subreads，對單分子多測序序列進行自我糾錯處理，獲得環形一致性序列（Circular consensus sequence，CCS）。通過檢測確定CCS 序列包含5′端引物、3′端引物以及poly-A 后進行分類，找出全長非嵌合（Full-length non-chimeric read，FLNC）序列。采用GMAP 軟件［16］，將FLNC 序列對比至黃顙魚基因組上［12］，再根據每條FLNC 序列的比對位置，統計分析基因座（loci）和轉錄本異構體（isoform）。另外，通過冗余轉錄本的去除和低可信度轉錄本的過濾獲得合格的isofrms。將測序得到的loci 和isoform 與參考基因組注釋的loci 和isoform 進行比較，可以確定檢測到已知基因新的isoform 以及鑒定到新基因的isoform。測序得到的基因滿足以下任一條即判定為新基因：①與已注釋基因沒有overlap 或overlap 小于20%；②與已注釋基因overlap 大于20%，但基因方向不一致。將本次測序獲得的轉錄本和參考基因組注釋得到的轉錄本進行比較分析，如果參考基因組注釋的基因轉錄本與三代轉錄組測序分析的isoform不同時為單外顯子，或轉錄組測序分析的isoform 存在1 個以上新的剪切位點，則認為該同源異構體是新的同源異構體。

1.4 IncRNA 預測和基因結構分析

將新基因的isoform、已知基因的新isoform 序列與NR、KOG、KO 庫比對，過濾掉潛在的編碼序列；對于在NR、KO 和KOG 庫中沒有hit 的序列，進一步利用CNCI、CPC2、CPAT 和PLEK 評估序列的編碼潛能，過濾編碼潛能大于設定的cutoff 或長度＜200 bp的序列，取4 個軟件預測結果的交集序列，作為最終的非編碼RNA 預測的結果［17，18］。PacBio 長讀長測序實現了全長轉錄本測序。相對于二代短讀長RNA-Seq 測序識別可變剪接時完全依賴于junction reads 比對的方法，三代全長測序使得直接基于全長isoform 序列相互比較的可變剪接識別成為可能。用ASprofile 軟件對測序得到的isoform 可變剪接事件分別進行分類和統計［19］。融合基因通過以下方式被確定：定位到2個或2個以上的遠距離范圍，定位比對至少占轉錄物的10%，覆蓋率≥99%，每個定位位點必須至少相距100 kb［20］。利用全長轉錄組APA 檢測軟件Tapis進行可變多聚腺苷酸化位點檢測［21］。

1.5 新基因功能注釋

使用以下公共數據庫對新基因進行基因功能注釋：非冗余蛋白數據庫（Non-redundant protein database，NR）、蛋白質真核同源數據庫（Eukaryotic orthologous groups，KOG）、基因本體論數據庫（Gene ontology，GO）、東京基因與基金組百科全書（Kyoto encyclopedia of genes and genomes，KEGG）和蛋白質序列數據庫（Swiss-Prot）。

2 結果與分析

2.1 測序結果與數據分析

全長轉錄組測序獲得9 525 639 個Subreads 和842 830 個環形一致性序列（CCS）Reads，其中，帶有5′端引物的reads 數有755 439個，有3′端引物reads數有772 185個，有Poly-A 的reads 數有760 944 個；全長（Full-length）reads 數共693 262個，全長非嵌合（Full-length non-chimeric read，FLNC）reads數有685 574個，全長非嵌合reads 平均長度為2 736 bp，全長非嵌合N50 為3 067 bp（表1）。采用GMAP 軟件將全長轉錄組FLNC 序列與黃顙魚基因組對比分析，基因組已注釋的loci 和isoform 均為24 552個，而基于PacBio 測得的FLNC 序列篩選后與基因組對比分析，鑒定出26 664 個loci 和72 509 個isoforms，isoforms 平均長度為2 918 bp。PacBio 分析的全部isoform 和基因組已注釋isoform 的長度分布和同一loci內isoform 的個數如圖1a、圖1b 所示。將72 509 個isoforms 與參考基因組比較分析發現，13.01%的isoforms（9 437 個）為已知基因isoforms，69.76%的isoforms（50 580 個）為已知基因新isoforms，17.23%的isoforms（12 492個）為新基因isoforms（圖1c）。

圖1 isoforms 統計和分類

表1 reads分類統計

2.2 LncRNA 分析

利用CNCI、CPC2、PLEK 和CPAT分別預測了6 497、6 482、6 390、3 379 個LncRNA，并將4 個軟件預測的交集部分（3 169 個LncRNA）作為最終的非編碼RNA 預測結果（圖2a）。根據LncRNA 在基因組上相對于蛋白編碼基因的位置分為4 種：正義鏈（sense）有861 個、反義鏈（antisense）658 個、內含子間（intronic）746 個、基因間（intergenic）904 個（圖2b），其所在的位置與其功能有一定的相關性。LncRNA 的長度分布見圖2c，最小長度為213 bp，最大長度為5 458 bp，平均長度為2 220 bp，序列的正態分布表明LncRNA 的序列質量較好。

圖2 黃顙魚長鏈非編碼RNA 分析

2.3 可變剪接事件、可變多聚腺苷酸化位點和融合基因分析

分析的可變剪接類型共有45 873個，分為以下幾種類型：外顯子跳躍（Exon skip，ES）有9 437個，占比20.57%；內含子保留（Retained intron，RI）有10 360個，占比22.58%；可變供體位點（Alternate donor site，AD）有2 182個，占比4.76%；可變受體位點（Alternate acceptor site，AA）有4 140個，占比9.03%；其他類型可變剪接形式有19 754個，占比43.06%（圖3a）。利用Tapis 軟件測到20 774 個polyA 位點來源于11 651 個基因（圖3b），共分析出4 881 個存在可變多聚腺苷酸化位點。此外，共檢測到304 個融合基因，其中，265 個融合基因來源于不同染色體上不同基因的融合，39 個融合基因來源于同一染色體上不同基因的融合（圖3c）。

圖3 全長轉錄本基因結構分析

2.4 新基因功能注釋

已發現的12 492 個新基因與公共數據進行功能注釋，成功注釋了7 233 個isoforms，其中，非冗余蛋白（NR）數據庫、基因本體論（GO）數據庫、東京基因與基金組百科全書（KEGG）數據庫、蛋白質真核同源（KOG）數據庫和蛋白質序列（Swiss-Prot）數據庫分別注釋到7 224、4 291、3 865、2 477、3 641 個新基因，還有5 259 個新基因未注釋上（圖4）。利用GO數據庫對4 291 條新基因進行注釋，并分類到生物學過程、細胞組分和分子功能，其中，細胞過程（Cell process，2 345 個）、細胞組分（Cell，1 979 個）和結合功能（Binding，2 195 個）分別在三大類中數量最多（圖5）。3 865 個新基因參與KEGG 代謝通路，并富集到細胞過程（Cellular processes）、環境信息處理（Environmental information processing）、遺傳信息處理（Genetic information processing）、新陳代謝（Metabolism）和生物體系統（Organismal system）上，其中，運輸和分解代謝（Transport and catabolism，343 個）、信號轉導（Signal transduction，712 個）、折疊分類和降解（Folding，sorting and degradation，194 個）、脂質代謝（Lipid metabolism，138 個）和免疫系統（Immune system，372 個）分別在五大類中數量最多（圖6）。

圖4 全長轉錄本新基因注釋

圖5 GO 注釋分類

此外，涉及黃顙魚生殖與繁殖相關的內分泌系統代謝途徑包括催產素信號通路（Oxytocin signaling pathway，70 個）、雌二醇信號通路（Estrogen signaling pathway，50 個）、孕酮介導的卵母細胞成熟（Progesterone-mediated oocyte maturation，44 個）、促性腺激素釋放激素信號通路（GnRH signaling pathway，30 個）和卵巢類固醇合成（Ovarian steroidogenesis，6 個）（圖7）。

圖7 KEGG 注釋的內分泌系統相關基因

3 小結與討論

在過去的研究中，由于缺乏黃顙魚基因組信息，關于黃顙魚遺傳和生理相關研究受到限制，只能通過操作繁瑣、效率較低的cDNA 末端快速克隆技術（RACE PCR）獲得黃顙魚部分基因轉錄本序列［22，23］。隨著高通量測序技術的快速發展，二代轉錄組測序技術廣泛應用于黃顙魚的研究［11，14，15］。由于缺乏基因組信息和二代轉錄組測序長度短等限制，絕大多數已發表的黃顙魚轉錄組均通過無參分析，導致基因注釋困難，對基因的可變剪切、融合基因和基因家族不能準確地檢測［1，24］。

本研究基于PacBio 平臺的單分子實時測序技術對黃顙魚的肝臟、腎、背部肌肉、腦、脾、心臟、皮膚、血液、鰓和性腺等組織進行混樣測序，共獲得全長非嵌合reads 數685 574個，全長非嵌合reads 平均長度為2 736 bp；相比二代測序技術獲得的黃顙魚轉錄本長度大幅提升，如Chen等［11］基于454 GSFLX 測序平臺獲得的黃顙魚混樣轉錄組unique sequences 平均長度僅601 bp，而Wu等［15］和Zhu等［25］基于Illumina 測序平臺的unigenes 平均長度分別為944 bp 和716 bp。LncRNA 是長度大于200 bp 的長鏈非編碼RNA，在生物體內廣泛存在，并介導許多復雜的生命活動過程［26，27］。在本研究中，利用CNCI、CPC2、PLEK 和CPAT 軟件預測到3 169 個LncRNA。在此基礎上，挖掘LncRNA 與內分泌和生殖過程的相關性，對其展開功能研究，對探索黃顙魚內分泌及生殖過程具有重要意義?；蚪Y構如可變剪接事件（AS）、可變多聚腺苷酸化位點（APA）和基因融合分析可增加轉錄多樣性和基因功能復雜性［28-30］。有些基因的前體mRNA（pre-mRNA）通過不同的剪接方式（選擇不同的剪接位點）產生不同的mRNA 剪接異構體，該過程稱為可變剪接（或選擇性剪接）（Alternative splicing，AS）。Weirather等［31］證明PacBio 測序平臺比二代轉錄組測序技術更有利于AS 事件的鑒定。在本研究中，共挖掘到45 873 個可變剪接事件。Poly-A 位點的改變也是一類重要的RNA 轉錄后調控修飾，產生具有不同UTRs 和編碼序列的mRNAs，其功能與選擇性剪接相似［32］。在本研究中，共挖掘到4 881 個基因存在可變多聚腺苷酸化位點。融合基因是指來源于不同基因的2 個片段被拼接在一起形成的新基因［33］。導致2 個基因發生融合的機制包括基因組結構變異、轉座或基因轉錄后的反式剪接等。在本研究中，共檢測到304 個融合基因，其中，265 個融合基因來源于不同染色體不同基因融合，39 個融合基因來源于同一染色體不同基因融合。

目前，新品種全雄黃顙魚“全雄1 號”和雜交黃顙魚“黃優1 號”已成為中國主流養殖品種。由于全雄黃顙魚均為雄性、雜交黃顙魚性腺退化無法繁殖、大型湖泊禁捕以及“長江十年禁漁計劃”等因素，導致中國黃顙魚母本資源短缺，苗種供應不足。人們利用魚類性逆轉技術結合黃顙魚性別連鎖分子標記，成功將XX 雌性黃顙魚逆轉為XX 雄性黃顙魚，然后XX 雄性和雌性黃顙魚繁殖后獲得黃顙魚全雌配套系［34-36］，黃顙魚全雌配套系還需要進一步提升其繁殖性能。Hu等［37］發現部分雌性黃顙魚存在排卵障礙問題，在常規催產藥物中加入鯉腦垂體提取物可順利排卵。23 號染色體上miR-200 簇敲除的斑馬魚存在輸卵管發育缺陷和排卵障礙問題，注射催產素（OXT 和AVT）+hCG 或鯉腦垂體提取物均可以促使排卵［38］，表明催產素在魚類的排卵過程中起著重要作用。本研究挖掘到的2 477 個新基因通過KEGG 分析富集到269 個已知途徑中，其中，292 個基因與動物內分泌系統相關，70 個基因富集到催產素信號、50 個基因富集到雌二醇信號通路、44 個基因富集到孕酮介導的卵母細胞成熟、30 個基因富集到促性腺激素釋放激素信號通路和6 個基因富集到卵巢類固醇合成。這些數據為以后的黃顙魚生殖和繁殖相關機制研究提供了科學依據。