基于高通量測(cè)序的黃姑魚轉(zhuǎn)錄組從頭組裝和基因注釋分析

王余菊，婁方瑞,2，水柏年

(1.浙江海洋大學(xué)水產(chǎn)學(xué)院，浙江舟山 316022；2.中國(guó)海洋大學(xué)水產(chǎn)學(xué)院，山東青島 266003)

黃姑魚Nibea albiflora，隸屬于鱸形目Perciformes、石首魚科Sciaenidea、黃姑魚屬Nibea。黃姑魚較為廣泛的分布在中國(guó)近海及日本和朝鮮半島海域，其作為東海較為重要的捕撈和養(yǎng)殖種類，具有重要的經(jīng)濟(jì)價(jià)值[1-2]。然而，近年來(lái)隨著黃姑魚捕撈量的增加，其自然資源量逐步下降，黃姑魚養(yǎng)殖業(yè)迅速發(fā)展，目前全國(guó)黃姑魚養(yǎng)殖年產(chǎn)量已超過(guò)3 萬(wàn)t[3-4]。國(guó)內(nèi)外學(xué)者在黃姑魚的基礎(chǔ)生物學(xué)、人工繁育、養(yǎng)殖生態(tài)和分子生物學(xué)等各方面已開展了較為廣泛研究[4-6]。值得注意的是，關(guān)于黃姑魚基因組學(xué)方面的數(shù)據(jù)信息較少，其基因組組裝結(jié)果尚未發(fā)表，限制了黃姑魚分子育種工作的開展。

轉(zhuǎn)錄組通常是指在生物的某一狀態(tài)下，由特定的組織或細(xì)胞所轉(zhuǎn)錄出的全部RNA 的集合，反應(yīng)了當(dāng)時(shí)基因表達(dá)的情況[7]。轉(zhuǎn)錄組測(cè)序能夠在缺少基因組信息的前提下，揭示基因與生物學(xué)特性和外界環(huán)境之間的關(guān)系，能夠有效地開發(fā)大量基因資源[8]。轉(zhuǎn)錄組測(cè)定是對(duì)于生物基因功能及其結(jié)構(gòu)進(jìn)行分析的首要環(huán)節(jié)，且其相比較于全基因組測(cè)序，具備更為迅速的測(cè)序速度和更低的測(cè)序成本。目前轉(zhuǎn)錄組測(cè)序并已在大黃魚Larimichthys crocea[9]、鮸魚Miichthys miiuy[10]、口蝦蛄Oratosquilla oratoria[11]、大竹蟶Solen grandis[12]和曼氏無(wú)針烏賊Sepiella maindroni[13]等多個(gè)海洋生物中得到廣泛應(yīng)用，ZHAN Wei，et al[14]報(bào)道了黃姑魚轉(zhuǎn)錄組測(cè)序組裝結(jié)果,Unigenes 的N50 為1 279 bp，組裝后轉(zhuǎn)錄本序列較短，這可能影響后續(xù)應(yīng)用。為了進(jìn)一步開發(fā)黃姑魚的基因資源，本研究中以采自舟山近海的黃姑魚自然個(gè)體為材料，利用Illumina Hiseq XTen平臺(tái)獲得黃姑魚轉(zhuǎn)錄組序列遺傳變異等信息，從中挖掘黃姑魚的基因數(shù)據(jù)和和微衛(wèi)星分子標(biāo)記，以期為后續(xù)黃姑魚及其他黃姑魚屬魚類的分子學(xué)探究提供數(shù)據(jù)基礎(chǔ)。

1 材料與方法

1.1 樣品采集

轉(zhuǎn)錄組測(cè)序所采用的黃姑魚樣品于2018 年10 月采自舟山附近海域，在20 ℃的實(shí)驗(yàn)室海水環(huán)境內(nèi)暫養(yǎng)6 d，暫養(yǎng)期間不投喂，采集活體黃姑魚的多種組織(肝、卵巢、鰓、眼和肌肉)分別用液氮速凍并保存在-80 ℃的超低溫環(huán)境中。

1.2 RNA 提取、文庫(kù)構(gòu)建和上機(jī)測(cè)序

液氮研磨并等量混勻黃姑魚組織，而后加入標(biāo)準(zhǔn)Trizol Reagent Kit 試劑獲取混勻后組織的總RNA?；诃傊悄z電泳、Nanodrop 分光光度計(jì)、Qubit 熒光定量?jī)x和Agilent 2100 生物分析儀對(duì)總RNA 的降解程度、純度、濃度和完整性分別進(jìn)行檢測(cè)。檢測(cè)合格后，基于磁珠法富集混勻后組織總RNA 中的mRNA，隨后將其打斷成特定短片段并以此為模板進(jìn)行文庫(kù)構(gòu)建。委托北京諾禾致源公司使用Illumina Hiseq XTen 平臺(tái)對(duì)庫(kù)檢合格的測(cè)序文庫(kù)進(jìn)行雙末端150 bp 測(cè)序。

1.3 原始測(cè)序數(shù)據(jù)過(guò)濾和質(zhì)量評(píng)估

使用Trimmomatic 軟件對(duì)于測(cè)序平臺(tái)中獲取的原始測(cè)序序列(Raw Reads) 進(jìn)行過(guò)濾以獲取高質(zhì)量的Clean Reads，其中，需要對(duì)帶有測(cè)序接頭的Raw Reads、Qphred<=20 的堿基達(dá)到所有堿基數(shù)量50%以上的Raw Reads 以及無(wú)法確定的堿基信息比例高于10%的Raw Reads 進(jìn)行剔除；剩余的高質(zhì)量Clean Reads 經(jīng)FastQC 軟件進(jìn)行質(zhì)量評(píng)估后用于后續(xù)分析。

1.4 轉(zhuǎn)錄本拼接和基因功能注釋

使用Trinity 2.4.0 軟件[15]對(duì)Clean Reads 進(jìn)行從頭組裝獲取轉(zhuǎn)錄本序列(Transcripts)，所使用的參數(shù)為：--min_kmer_cov:3?；趍ap 到轉(zhuǎn)錄本序列上的reads 數(shù)及表達(dá)模式，Corset 1.05 軟件[16]進(jìn)一步對(duì)獲取的轉(zhuǎn)錄本(Transcripts)進(jìn)行層次聚類，從而獲取最長(zhǎng)的聚類序列(Unigenes)，軟件參數(shù)使用默認(rèn)參數(shù)。Unigenes序列用于后續(xù)分析。而后，利用Diamond 和Blast 等比對(duì)軟件基于七大基因注釋數(shù)據(jù)庫(kù)Nr (Non-Redundant Protein Sequence Database)、Nt (Nucleotide Sequence Database)、Pfam (Pfam protein families database)、KOG(Clusters of Orthologous Groups of proteins)、Swiss-prot,KEGG(Kyoto Encyclopedia of Genes and Genomes)和GO (Gene Ontology)對(duì)所有的Unigene 進(jìn)行比對(duì)和注釋，從而獲得黃姑魚的表達(dá)基因功能，基因功能注釋軟件及參數(shù)設(shè)置如表1。Diamond 是一個(gè)用于比對(duì)數(shù)據(jù)庫(kù)蛋白，其蛋白質(zhì)比對(duì)性能為blast+的500～20 000 倍，因此比對(duì)Nr、KOG 和Swiss-prot 等3 個(gè)蛋白數(shù)據(jù)庫(kù)采用Diamond 軟件。

表1 基因功能注釋軟件及參數(shù)Tab.1 The software and parameters associated with gene functional annotation.

1.5 編碼序列(CDS)預(yù)測(cè)與簡(jiǎn)單重復(fù)序列標(biāo)記(SSR)位點(diǎn)搜索

編碼序列(cording sequences，CDS)預(yù)測(cè)分析步驟如下[11]：首先，將聚類獲取的Unigenes 依次比對(duì)到Nr和Swissprot 蛋白庫(kù)并獲取比對(duì)成功的Unigenes 中的編碼框核苷酸序列；第二，應(yīng)用estscan3.0.3 軟件預(yù)測(cè)未比對(duì)成功的Unigenes 中的編碼序列核苷酸信息；最后，按照5'->3'的順序，將編碼閱讀框堿基序列翻譯為氨基酸序列。轉(zhuǎn)錄中的微衛(wèi)星序列(simple sequence repeats，SSR)則主要依靠MISA 1.0 軟件進(jìn)行檢測(cè)[13]，軟件設(shè)置使用默認(rèn)參數(shù)，其中，重復(fù)單位及其最少重復(fù)次數(shù)分別設(shè)定為：?jiǎn)魏塑账嶂貜?fù)10 次、雙核苷酸重復(fù)6 次、三核苷酸重復(fù)5 次、四核苷酸重復(fù)5 次、五核苷酸重復(fù)5 次和六核苷酸重復(fù)5 次。

2 結(jié)果

2.1 黃姑魚轉(zhuǎn)錄組測(cè)序質(zhì)量評(píng)價(jià)和序列組裝

轉(zhuǎn)錄組測(cè)序共獲取43 385 432 條Raw Reads，經(jīng)過(guò)去除低品質(zhì)的原始序列，得到42 809 266 條高質(zhì)量的Clean Reads。質(zhì)量評(píng)估結(jié)果表明，Clean Reads 的堿基錯(cuò)誤率、Q20、Q30 和GC 含量值分別為0.02%、97.15%、92.28%和49.27%。Clean Reads 拼接為57 654 條Transcripts，Transcripts 的平均長(zhǎng)度和N50 值分別為1 950 bp 和3 255 bp；所有的transcripts 進(jìn)一步聚類后獲得32 623 條Unigenes，Unigenes 的平均長(zhǎng)度和N50 值分別為1 646 bp 和2 777 bp。Trianscripts 和Unigenes 在不同長(zhǎng)度區(qū)間內(nèi)的頻數(shù)統(tǒng)計(jì)結(jié)果如圖1 所示。

圖1 黃姑魚轉(zhuǎn)錄組Transcripts 和Unigenes 長(zhǎng)度分布Fig.1 Distribution of the length of transcripts and unigene in the N.albiflora transcriptome.

2.2 基因功能注釋

2.2.1 基因功能注釋成功率統(tǒng)計(jì)

基于7 大數(shù)據(jù)庫(kù)對(duì)黃姑魚轉(zhuǎn)錄組的Unigenes 序列進(jìn)行基因功能分析，注釋成功率如表2 所示。結(jié)果表明，28 645 條Unigenes(87.81%)匹配到至少一個(gè)數(shù)據(jù)庫(kù)中，占總Unigene 的；7 023 條Unigenes(21.52%)可以匹配到所有數(shù)據(jù)庫(kù)中。

表2 注釋結(jié)果統(tǒng)計(jì)Tab.2 The annotation results.

2.2.2 Unigene 序列相似性分析

為了獲取黃姑魚基因序列及其功能信息，并進(jìn)一步確定黃姑魚與其近緣種基因序列的相似度，將轉(zhuǎn)錄組中的Unigenes 與NCBI 中的Nr 數(shù)據(jù)庫(kù)進(jìn)行比對(duì)。結(jié)果表明，黃姑魚轉(zhuǎn)錄組中16 291、1 106、665、370 和346 條Unigenes 分別與大黃魚、尖吻鱸Lates calcarifer、高體鰤Seriola dumerili、貝氏隆頭魚Labrus bergylta 和眼斑雙鋸魚Amphiprion ocellaris的基因序列具有較高的相似性，剩余的3 085 條Unigenes 與132 種其他物種的基因序列具有相似性(圖2)。

2.2.3 Unigene 功能分類

9 605 條在KOG 數(shù)據(jù)庫(kù)中注釋成功黃姑魚轉(zhuǎn)錄組的Unigene 被歸屬于26 個(gè)KOG 功能大類中(圖3)。其中，數(shù)量最多的功能類為信號(hào)轉(zhuǎn)導(dǎo)機(jī)制(1 843 條)，而后依次為一般功能預(yù)測(cè)(1 507 條)、轉(zhuǎn)錄后修飾、蛋白折疊和分子伴侶(946 條)等。

圖2 Nr 數(shù)據(jù)庫(kù)中Unigenes 注釋信息Fig.2 The annotation information of all unigenes blast in Nr database

圖3 黃姑魚轉(zhuǎn)錄組Unigenes 的KOG 分類Fig.3 KOG classification of unigenes in the N.albiflora transcriptome

此外，與GO 數(shù)據(jù)庫(kù)進(jìn)行比對(duì)結(jié)果表明，17 812 條匹配成功的Unigenes 被富集在1 590 個(gè)GO 功能詞條中，如圖4。其中，在3 個(gè)GO 大類生物過(guò)程(Biological process)、細(xì)胞成分(Cellular component)和分子功能(Molecular Function)中分別有54 612、33 293 和22 539 條Unigenes 獲得注釋信息。在生物過(guò)程大類中，富集在細(xì)胞過(guò)程(GO：0008152)最多Unigenes 數(shù)量較多，達(dá)10 595 條，其次是代謝過(guò)程(GO：0009987)詞條為8 527 條；在細(xì)胞成分大類中，富集在細(xì)胞(GO：0005623)和細(xì)胞組分(GO：0044464)詞條中的Unigenes 數(shù)量較多，均為6 042 條；在分子功能大類中，富集在結(jié)合(GO：0005488)最多為11 082 條，其次是催化活性(GO：0003824)中詞條為6 928 條。

圖4 黃姑魚轉(zhuǎn)錄組Unigenes 的GO 分類Fig.4 GO classification of unigenes in the N.albiflora transcriptome.

2.2.4 Unigene 代謝途徑分析

參與同一代謝通路的黃姑魚的Unigene 可以通過(guò)KEGG 分析富集在一起。黃姑魚轉(zhuǎn)錄組數(shù)據(jù)的代謝途徑分析結(jié)果表明，所有的Unigenes 與KEGG 數(shù)據(jù)后進(jìn)行比對(duì)后，14 754 條Unigene 顯著富集在231 種代謝通路中。其中，各大類中富集Unigene 數(shù)量較多的通路如圖5 所示，涉及信號(hào)轉(zhuǎn)導(dǎo)、內(nèi)分泌系統(tǒng)、免疫系統(tǒng)等。

圖5 黃姑魚轉(zhuǎn)錄組Unigenes 的KEGG 分類Fig.5 The KEGG pathways of unigenes in the N.albiflora transcriptome (A,Cellular Processes; B,Environmental Information Processing; C,Genetic Information Processing; D,Metabolism; E,Organismal Systems).

2.3 CDS 預(yù)測(cè)和SSR 位點(diǎn)分析

通過(guò)Nr 和Swissprot 數(shù)據(jù)庫(kù)，16 803 個(gè)CDS 從比對(duì)成功的Unigenes 中被提??；使用estscan 對(duì)未比對(duì)成功的Unigene 進(jìn)一步分析后獲取13 579 個(gè)預(yù)測(cè)CDS。同時(shí)，我們共獲取了29 022個(gè)SSR 標(biāo)記，這些標(biāo)記位于13 508 條Unigenes 中，如圖6 所示。結(jié)果表明，SSR 標(biāo)記的數(shù)量與SSR 標(biāo)記的重復(fù)單位數(shù)量成負(fù)相關(guān)，單堿基重復(fù)類型SSR 標(biāo)記數(shù)量最多，為14 702 個(gè)；六堿基重復(fù)類型SSR 標(biāo)記數(shù)量最少，僅為40 個(gè)。從SSR 重復(fù)基序角度來(lái)說(shuō)，AC 重復(fù)是黃姑魚最普遍的SSR 標(biāo)記，此外，最普遍的三堿基重復(fù)是AGG 重復(fù)，AAAC 重復(fù)則是最普遍的四堿基重復(fù)。

圖6 黃姑魚轉(zhuǎn)錄組中SSR 分析統(tǒng)計(jì)結(jié)果Fig.6 Number of SSRs in the N.albiflora transcriptome

3 討論

3.1 黃姑魚轉(zhuǎn)錄組測(cè)序與組裝

近年來(lái)，轉(zhuǎn)錄組測(cè)序的優(yōu)勢(shì)性已經(jīng)促使其在許多缺乏基因組信息的海洋生物中得到廣泛應(yīng)用。本次研究對(duì)黃姑魚的多種組織進(jìn)行混勻后首次測(cè)定了其轉(zhuǎn)錄組信息，測(cè)序結(jié)果豐富和完善了黃姑魚的基因數(shù)據(jù)庫(kù)資源。對(duì)原始序列中的低質(zhì)量序列經(jīng)過(guò)剔除后，約有98.67%的高質(zhì)量序列。此外，對(duì)獲取的高質(zhì)量序列進(jìn)行質(zhì)控發(fā)現(xiàn)，序列的堿基錯(cuò)誤率較低且具有相對(duì)較高的質(zhì)量值。這進(jìn)一步證明了本次黃姑魚的轉(zhuǎn)錄組測(cè)序信息的準(zhǔn)確性相對(duì)較高[22]。黃姑魚轉(zhuǎn)錄組拼接出的Transcripts 和Unigenes 的N50 值分別為3 255 bp 和2 777 bp，同之前的黃姑魚轉(zhuǎn)錄組組裝結(jié)果相比有明顯提高(UnigenesN50 值為1 279 bp)[14]，表明黃姑魚的轉(zhuǎn)錄本組裝質(zhì)量較好，適用于后續(xù)的基因功能分析。

3.2 功能注釋

在轉(zhuǎn)錄組數(shù)據(jù)挖掘中，基因功能注釋是極為關(guān)鍵的環(huán)節(jié)之一。本次研究中，約有87.81%的Unigenes 在至少一個(gè)數(shù)據(jù)庫(kù)中成功獲得功能注釋，較高的注釋率再次驗(yàn)證了本次轉(zhuǎn)錄組測(cè)序的大數(shù)據(jù)量和轉(zhuǎn)錄本拼接的有效性。此外，蛋白數(shù)據(jù)庫(kù)中石首魚科蛋白質(zhì)信息較為豐富也是黃姑魚轉(zhuǎn)錄組高注釋率的原因。Nr 注釋結(jié)果表明，大量的Unigenes 與大黃魚的基因序列得到了匹配，這是因?yàn)榇簏S魚的基因組測(cè)序早已完成，且其與黃姑魚均屬于石首魚科[2]。另外，功能聚類和通路富集的分析也為我們快速挖掘與黃姑魚生物學(xué)特性相關(guān)的基因奠定了基礎(chǔ)。

3.3 分子標(biāo)記篩選

在物種鑒定、群體遺傳學(xué)及遺傳圖譜構(gòu)建等研究中，分子標(biāo)記開發(fā)可以提供較為有價(jià)值的基礎(chǔ)信息。有研究認(rèn)為，基因編碼區(qū)域中可能存在著廣泛的分子標(biāo)記，而轉(zhuǎn)錄組測(cè)序獲取的大的信息量實(shí)際上可能涵蓋了成千上萬(wàn)的表達(dá)基因，因此，轉(zhuǎn)錄組測(cè)序技術(shù)可能在分子標(biāo)記開發(fā)方面具有較大的優(yōu)勢(shì)性[7]。本研究中，MISA 軟件共搜索到29 022 個(gè)SSR 位點(diǎn)，從而進(jìn)一步豐富了黃姑魚的分子標(biāo)記資源，但所有的SSR 的有效性仍需進(jìn)一步驗(yàn)證。

4 結(jié)論

本研究初次對(duì)黃姑魚的轉(zhuǎn)錄組信息進(jìn)行了測(cè)序，捕獲的轉(zhuǎn)錄組數(shù)據(jù)有效性和可信度較高，從而豐富了黃姑魚的基因數(shù)據(jù)庫(kù)資源。但約有13%的Unigenes 沒有得到注釋，這些可能是黃姑魚特有的基因轉(zhuǎn)錄本或是由于組裝錯(cuò)誤形成的轉(zhuǎn)錄本，需要后續(xù)三代測(cè)序數(shù)據(jù)進(jìn)行進(jìn)一步驗(yàn)證。此外，轉(zhuǎn)錄本的注釋結(jié)果也為黃姑魚生物學(xué)調(diào)控機(jī)制的基因功能、分子事件和信號(hào)通路奠定了基礎(chǔ)。大量的SSR 位點(diǎn)的獲取可為后續(xù)黃姑魚分子標(biāo)記的開發(fā)提供了參考。黃姑魚轉(zhuǎn)錄本數(shù)據(jù)獲得同時(shí)也為石首魚科物種分化研究、種群適應(yīng)進(jìn)化研究、遺傳多樣性研究和系統(tǒng)發(fā)育學(xué)研究提供了數(shù)據(jù)保障，為黃姑魚屬種質(zhì)資源的開發(fā)利用奠定了基礎(chǔ)。