低鹽脅迫下金烏賊轉(zhuǎn)錄組學(xué)分析

王興強(qiáng)，曹 梅，崔春輝，鄭年昊，沈 曄，劉豪杰，丁雪峰

( 江蘇海洋大學(xué) 海洋生命與水產(chǎn)學(xué)院，江蘇 連云港 222005 )

金烏賊(Sepiaesculenta)又名墨魚、大魷魚，屬軟體動(dòng)物門、頭足綱、鞘亞綱、烏賊目、烏賊科、烏賊屬，屬廣溫性種類，廣泛生活在渤海、黃海、東海、南海和日本列島附近海域，有洄游習(xí)性，對(duì)鹽度要求較為苛刻[1]。金烏賊喜歡生活在幾十米深的沙質(zhì)水域和水質(zhì)良好、流速適中的島嶼附近，偶爾穴居，趨光性強(qiáng)，白天潛底躲避敵害生物，夜晚上浮至海水上層捕食。自然狀況下，金烏賊1年左右即達(dá)到性成熟；人工養(yǎng)殖條件下，若采取升溫、過量投喂、鹽度調(diào)整等措施可加快性腺發(fā)育。金烏賊雌雄異體，一年僅繁殖一次，繁殖時(shí)有較復(fù)雜的求偶行為，產(chǎn)卵前雌性有噴沙習(xí)性，雌性烏賊每次可產(chǎn)卵莢約1500粒，雄性烏賊則排出約800個(gè)精莢，產(chǎn)卵、排精后親體很快死亡[2]。金烏賊幼時(shí)以捕食端足類和橈足類為生，成體則捕食各種甲殼類和魚類為主，在缺少食物時(shí)，偶爾會(huì)發(fā)生同類相殘的現(xiàn)象。近年來，因?yàn)檫^度捕撈、海洋污染和天然產(chǎn)卵場、棲息地的不斷萎縮，金烏賊生物資源量急劇下降。為了恢復(fù)金烏賊種質(zhì)資源，人工養(yǎng)殖迫在眉睫[3]。與傳統(tǒng)的魚蝦蟹類不同，金烏賊養(yǎng)殖過程中只攝食活的餌料，其喜歡噴墨的生物學(xué)習(xí)性為人工養(yǎng)殖增加了很大難度，加之在養(yǎng)殖收獲和活體運(yùn)輸方面沒有天然優(yōu)勢，導(dǎo)致金烏賊的人工養(yǎng)殖停滯不前。鑒于上述原因，研究環(huán)境脅迫對(duì)金烏賊生理機(jī)制的影響具有重要意義。筆者通過比較正常鹽度(30)與低鹽度脅迫(15)處理6 h的金烏賊幼體轉(zhuǎn)錄組學(xué)數(shù)據(jù)，探討金烏賊響應(yīng)低鹽脅迫的分子機(jī)制，獲得低鹽脅迫條件下金烏賊生長、生殖、代謝和免疫方面相關(guān)的基因資源，為以后金烏賊分子標(biāo)記的開發(fā)和關(guān)鍵基因的克隆提供技術(shù)支撐。

1 材料與方法

1.1 試驗(yàn)材料

試驗(yàn)所用金烏賊均為通過人工孵化所得幼體，人工育苗所用金烏賊親體來自海州灣海域，親體暫養(yǎng)及育苗工作均在連云港市連云區(qū)核電南路6號(hào)大板橋試驗(yàn)基地育苗室中完成。

1.2 低鹽脅迫試驗(yàn)

孵化約30 d的金烏賊幼體，體質(zhì)量(1.3±0.3) g，苗種培育鹽度30；據(jù)預(yù)試驗(yàn)得知，金烏賊在鹽度15的條件下可存活6～8 h，因此低鹽脅迫時(shí)間選擇6 h。試驗(yàn)結(jié)束后分別取對(duì)照組和低鹽脅迫組20只金烏賊幼體，直接投入液氮速凍。

1.3 RNA提取、cDNA文庫構(gòu)建及原始測序數(shù)據(jù)獲取

按照Trizol法提取總RNA，并用Takara miniBEST試劑盒進(jìn)行純化，運(yùn)用Nanodrop檢測得到的RNA，將達(dá)到要求的RNA樣品液氮包埋，送至上海元莘生物進(jìn)行Illumina高通量測序，得到的原始圖像數(shù)據(jù)經(jīng)過Base Calling轉(zhuǎn)化為序列數(shù)據(jù)，以FASTQ文件格式存儲(chǔ)。

1.4 原始測序數(shù)據(jù)處理及轉(zhuǎn)錄組從頭拼接

對(duì)原始數(shù)據(jù)進(jìn)行過濾，獲得高質(zhì)量測序數(shù)據(jù)。使用SeqPrep(https:∥github.com/jstjohn/SeqPrep)和 Sickle(https:∥github.com/najoshi/sickle)對(duì)質(zhì)量修剪前后的序列進(jìn)行數(shù)據(jù)量統(tǒng)計(jì)。利用Trinity(http:∥trinityrnaseq.sourceforge.net/)對(duì)高質(zhì)量的待分析數(shù)據(jù)進(jìn)行從頭組裝，對(duì)組裝得到的轉(zhuǎn)錄組進(jìn)行各項(xiàng)基礎(chǔ)指標(biāo)的統(tǒng)計(jì)。

1.5 從頭組裝轉(zhuǎn)錄組功能注釋

利用Trinity軟件提供的ORF預(yù)測流程(http:∥trinityrnaseq.sourceforge.net/analysis/extract_proteins_from_trinity_transcripts.html)對(duì)組裝得到的所有轉(zhuǎn)錄本和Unigenes序列進(jìn)行基因預(yù)測，運(yùn)用HMMER軟件和Pfam數(shù)據(jù)庫(http:∥pfam.sanger.ac.uk/)對(duì)照，得到完整的Unigenes蛋白家族注釋信息。將拼接得到的所有核苷酸序列，使用BlastX(版本2.2.25)分別與NR(Non-redundant protein sequences from GenPept, Swissprot, PIR, PDF, PDB, and NCBI RefSeq)、Swiss-Prot(Swiss-Prot protein sequence database)、GO(Gene ontology)、COG(Cluster of orthologous groups of proteins)、KOG(Clusters of orthologous groups for eukaryotic complete genomes)和KEGG(Kyoto encyclopedia of genes and genomes)數(shù)據(jù)庫進(jìn)行比對(duì)，獲得Unigenes相應(yīng)的功能注釋信息。

1.6 基因表達(dá)量估計(jì)

利用Bowtie(http:∥bowtie-bio.sourceforge.net/index.shtml)將對(duì)照和低鹽脅迫處理組高質(zhì)量的測序數(shù)據(jù)分別比對(duì)到組裝完成的轉(zhuǎn)錄組序列上，RSEM(http:∥deweylab.biostat.wisc.edu/rsem/)中利用Bowtie的比對(duì)結(jié)果進(jìn)行表達(dá)量估計(jì)并進(jìn)行FPKM轉(zhuǎn)換，從而得到基因的表達(dá)水平?；虿町惐磉_(dá)分析使用EdgeR(http:∥www.bioconductor.org/packages/2.12/bioc/html/edgeR.html)[4-6]。采用Goatools (https:∥github.com/tanghaibao/GOatools)進(jìn)行差異表達(dá)基因GO富集分析和Fisher精確檢驗(yàn)。采用KOBAS(http:∥kobas.cbi.pku.edu.cn/home.do)進(jìn)行差異表達(dá)基因KEGG通路富集分析和Fisher精確檢驗(yàn)，采用BH (FDR)方法進(jìn)行多重檢驗(yàn)。

顯著差異基因篩選的表達(dá)公式如下：

差異表達(dá)基因的校驗(yàn)公式：FDR<0.05

(1)

差異表達(dá)基因的校驗(yàn)公式：log2|FC|1

(2)

式中，F(xiàn)DR為基因/轉(zhuǎn)錄本在兩樣本中差異顯著性檢驗(yàn)結(jié)果，log2|FC|為該基因/轉(zhuǎn)錄本在兩樣間差異倍數(shù)以2為底的對(duì)數(shù)值，同時(shí)滿足公式(1)與公式(2)視為該基因/轉(zhuǎn)錄本在兩樣本間的差異表達(dá)顯著。

2 結(jié) 果

對(duì)正常鹽度(30，對(duì)照組)和低鹽脅迫(15)兩組處理進(jìn)行轉(zhuǎn)錄組測序，測序共獲得87 326 026條序列，質(zhì)量剪切后對(duì)照組序列43 680 773條，低鹽脅迫組40 494 301條，Q30超過91.14%，GC含量超過38.59%(表1)。

表1 測序數(shù)據(jù)統(tǒng)計(jì)

注：Q20、Q30. Phred值>20、30的堿基占總體堿基的百分比.

Note: Q20. Percentage of bases whose phred value is greater than or equal to 20; Q30. Percentage of bases whose phred value is greater than or equal to 30.

經(jīng)過原始測序數(shù)據(jù)剪切和Trinity軟件從頭拼接獲得575 171條轉(zhuǎn)錄本，平均長度349.09 bp，N50為326 bp；獲得513 053條Unigenes，平均長度332.25 bp，N50為310 bp(表2)。

表2 組裝結(jié)果統(tǒng)計(jì)

注：取最長的轉(zhuǎn)錄本作為Unigene；N50.按照長度將組裝轉(zhuǎn)錄本/Unigenes由大到小排序，累加轉(zhuǎn)錄本/Unigenes的長度到總長度的一半時(shí)，對(duì)應(yīng)轉(zhuǎn)錄本/Unigenes的長度.

Note:take the longest transcript as Unigene; N50. order of sequence of the assembled transcripts/Unigenes from large to small according to their length, and add up the length of the transcripts/Unigenes to half of the total length, corresponding to the length of the transcripts or Unigenes.

分別在NR、Swiss-Prot、KEGG、String和Pfam數(shù)據(jù)庫對(duì)Unigenes進(jìn)行功能注釋，COG、KOG和NOG數(shù)據(jù)庫共注釋Unigenes 4179條，其中COG注釋2709條、KOG注釋3416條和NOG注釋156條(表3)。

表3 金烏賊Unigenes注釋成功率

利用GO數(shù)據(jù)庫，7232條Unigenes按照參與生物過程、分子功能和細(xì)胞組分3個(gè)方面進(jìn)行分類注釋。生物過程包括25個(gè)分支，細(xì)胞組分包括20個(gè)分支和分子功能包括16個(gè)分支(圖1)。生物過程得到的GO注釋信息數(shù)量最多(19 730條)，其次是細(xì)胞組分(10 525條)，注釋信息最少的是分子功能(8645條)。生物過程主要包含代謝過程(4205條)、細(xì)胞過程(4126條)、單生物過程(3116條)和生物調(diào)節(jié)(1473條)相關(guān)的GO注釋。分子功能主要包含結(jié)合(3732條)、催化活性(3514條)、轉(zhuǎn)運(yùn)活性(524條)和分子傳感器活性(207條)相關(guān)的GO注釋。

隨著層數(shù)的增加，對(duì)分子功能、生物過程和細(xì)胞組分的描述也就更加詳細(xì)，基因功能注釋更加清晰(圖2)。

COG數(shù)據(jù)庫將Unigenes分為25類，其中，通用功能預(yù)測含有的Unigenes數(shù)量最多(459條)；接下來依次是復(fù)制、重組和修復(fù)(175條)、轉(zhuǎn)錄(143條)、能源生產(chǎn)和轉(zhuǎn)化(126條)和氨基酸運(yùn)輸和代謝(106條)，沒有序列注釋到[W] 細(xì)胞外結(jié)構(gòu)上(圖3)。

本次在KEGG數(shù)據(jù)庫中獲得有注釋的Unigenes 13 835條(2.7%)，圖4列出含有基因數(shù)量最多的前20個(gè)通路，包含Unigenes數(shù)目較多的KEGG通路有嘌呤、嘧啶和碳代謝，PI3K-AKT、cAMP和Rap1信號(hào)通路，內(nèi)吞作用，RNA轉(zhuǎn)運(yùn)，局灶性黏附，賴氨酸降解和泛素介導(dǎo)的蛋白水解等。

圖1 Unigenes GO二級(jí)分類統(tǒng)計(jì)

圖2 Unigenes GO二、三、四級(jí)分類統(tǒng)計(jì)

圖3 Unigenes COG分類統(tǒng)計(jì)

圖4 包含Unigenes數(shù)目最多的前20個(gè)通路

根據(jù)涉及的KEGG代謝通路，將基因分為代謝(A)、遺傳信息處理(B)、環(huán)境信息處理(C)、細(xì)胞過程(D)和有機(jī)系統(tǒng)(E)5個(gè)分支(圖5)。其中，代謝分支中的基因數(shù)占比最大，達(dá)2545條。在這分支內(nèi)部，概述網(wǎng)包含的基因最多(1914條)，其余依次是氨基酸代謝(835條)、碳水化合物代謝(500條)、脂質(zhì)代謝(419條)和能量代謝(385條)。排名第二的是遺傳信息處理，共有1989條。排名第三的是有機(jī)系統(tǒng)，共有1858條，在其內(nèi)部內(nèi)分泌系統(tǒng)包含的基因最多(779條)，其余依次為免疫系統(tǒng)(504條)和消化系統(tǒng)(424條)。排名第四的是環(huán)境信息處理，共有1783條信息，其中信號(hào)轉(zhuǎn)導(dǎo)包含的基因數(shù)目最多(1462條)。數(shù)量最少的是細(xì)胞過程，共有1502條，其中運(yùn)輸和分解代謝(750條)、細(xì)胞生長和死亡(575條)數(shù)量較多。

圖5 Unigenes KEGG注釋統(tǒng)計(jì)

對(duì)比到NR數(shù)據(jù)庫中的19 932條Unigenes中，其中注釋到長牡蠣(Crassostreagigas)2886條，注釋到霸王蓮花青螺(Lottlagigantea)2732條，注釋到海蝸牛(Aplysiacalifornica)1111條，注釋到海蠕蟲(Capitellateleta)742條，注釋到豬毛尾線蟲(Trichurissuis)536條，注釋到夏威夷四盤耳烏賊(Euprymnascolopes)157條，注釋到其余物種6378條，注釋到未知物種3043條(圖6)。

圖6 Unigenes物種分類

8 080 012條對(duì)照組和6 169 906條低鹽脅迫組的高質(zhì)量測序數(shù)據(jù)分別比對(duì)到組裝得到的轉(zhuǎn)錄組序列上，定位百分比分別為24%和20%(表4)。

表4 與轉(zhuǎn)錄組比對(duì)統(tǒng)計(jì)

低鹽脅迫產(chǎn)生1923條差異表達(dá)基因(圖7、圖8)，其中1114條顯著上調(diào)，809條顯著下調(diào)。

GO功能富集分析顯示，一些可能與低鹽脅迫相關(guān)的生物學(xué)過程如α-氨基酸代謝、羧酸代謝、氧乙酸代謝、有機(jī)酸代謝、RNA代謝和發(fā)育等過程得到了顯著富集(圖9)。

圖7 基因表達(dá)量分布

圖8 上下調(diào)基因GO注釋柱形

圖9 差異表達(dá)基因GO富集柱狀圖

GO可視化分析發(fā)現(xiàn)，低鹽脅迫對(duì)α-氨基酸代謝、水解酶活性、金屬離子、陰離子及核苷酸結(jié)合等過程影響顯著(圖10)。

KEGG通路富集分析顯示，低鹽脅迫6 h后的差異表達(dá)基因主要富集到雌激素信號(hào)通路、心肌細(xì)胞腎上腺素信號(hào)通路、類固醇生物合成、抗原加工與表達(dá)、脂肪消化與吸收、蛋白質(zhì)消化與吸收、甘油脂質(zhì)代謝、花生四烯酸代謝和酪氨酸代謝等信號(hào)通路上(圖11)，主要差異表達(dá)基因及所在通路見表5。

圖10 顯著性GO有向無環(huán)圖

圖11 差異表達(dá)基因KEGG富集柱狀圖

表5 主要差異表達(dá)基因及所在信號(hào)通路

(續(xù)表5)

3 討 論

近年來，高通量測序技術(shù)廣泛應(yīng)用在眾多動(dòng)植物的轉(zhuǎn)錄組分析方面，成果豐富。在對(duì)植物的研究中，玉米、水稻和大豆等糧食作物占重要地位，學(xué)者們對(duì)魚、蝦、蟹、貝等水產(chǎn)動(dòng)物也做了大量探索[7-8]。金烏賊屬狹鹽性的養(yǎng)殖品種，不能適應(yīng)低鹽度的海水，鹽度一旦低于24，會(huì)出現(xiàn)大量死亡的現(xiàn)象[3]。低鹽脅迫下，金烏賊生長速度減緩、免疫力下降，脅迫嚴(yán)重時(shí)會(huì)導(dǎo)致大規(guī)模死亡。因此，在金烏賊養(yǎng)殖過程中，尤其是汛期發(fā)生時(shí)，要嚴(yán)格控制水體鹽度的變化。當(dāng)前宜研究鹽度對(duì)金烏賊影響的生理機(jī)制，挖掘響應(yīng)鹽度脅迫的敏感指標(biāo)。本研究運(yùn)用高通量測序技術(shù)，對(duì)金烏賊幼體進(jìn)行測序，共獲得87 326 026條序列，經(jīng)過原始測序數(shù)據(jù)質(zhì)量剪切和從頭拼接得到575 171條轉(zhuǎn)錄本和513 053條Unigenes(表1、表2)。Pfam數(shù)據(jù)庫構(gòu)建了多個(gè)不同族序的氨基酸序列的統(tǒng)計(jì)模型，蛋白質(zhì)由一個(gè)或多個(gè)結(jié)構(gòu)域構(gòu)成，不同的結(jié)構(gòu)域賦予蛋白質(zhì)不同的功能，通過識(shí)別蛋白質(zhì)的結(jié)構(gòu)域序列，以此來預(yù)測蛋白質(zhì)的作用[9]。分別在NR、Swiss-Prot、KEGG、String和Pfam數(shù)據(jù)庫對(duì)Unigenes進(jìn)行功能注釋，但絕大部分Unigenes未獲得注釋，定位不清(表3)，這可能是由于金烏賊有墨囊，在一定程度上影響RNA提取及轉(zhuǎn)錄的效果；而金烏賊為海產(chǎn)動(dòng)物，雜合度較高，導(dǎo)致拼裝的Unigenes序列過短；Unigenes的N50僅為310，較一般轉(zhuǎn)錄組的組裝水平差。通過與NR庫的比對(duì)，由注釋的情況看，金烏賊與貝類品種長牡蠣、霸王蓮花青螺和海蝸牛親緣關(guān)系較近(圖6)。由于頭足類的研究剛剛起步，還需要更多的轉(zhuǎn)錄組信息豐富到NR數(shù)據(jù)庫。與String數(shù)據(jù)庫比對(duì)結(jié)果，可獲得基因?qū)?yīng)的COG數(shù)量，根據(jù)COG數(shù)量對(duì)所有轉(zhuǎn)錄本進(jìn)行功能歸類[10-11]。注釋基因在經(jīng)過代謝路徑分析與功能分類后，COG功能注釋分為25類，合計(jì)有2709條Unigenes(圖3)。

關(guān)于鹽度對(duì)水產(chǎn)動(dòng)物基因表達(dá)的影響，在紅耳龜(Trachemysscriptaelegans)、尼羅羅非魚(Oreochromisniloticus)、凡納濱對(duì)蝦(Litopenaeusvannamei)、中華絨螯蟹(Eriocheirsinensis)、香魚(Plecoglossusaltivelis)、條紋鲇(Pangasianodonhypophthalmus)、長牡蠣和文蛤(Meretrixmeretrix)等已有諸多研究[12-19]。本研究中，對(duì)照組和低鹽脅迫組的高質(zhì)量測序數(shù)據(jù)比對(duì)到組裝得到的轉(zhuǎn)錄組序列上的定位百分比分別為24%和20%，可能與拼裝的Unigenes序列過短有關(guān)，導(dǎo)致比對(duì)率較低。GO數(shù)據(jù)庫有助于理解基因背后所代表的生物學(xué)意義，GO功能顯著性富集分析可以說明差異基因的功能富集情況，在基因功能水平闡明樣本間的差異[20]。由圖1可見，7232條Unigenes按照其參與的生物過程、分子功能及細(xì)胞組分3個(gè)方面進(jìn)行了分類注釋。根據(jù)圖9和圖10推測，低鹽脅迫對(duì)金烏賊α-氨基酸代謝、羧酸代謝、氧乙酸代謝、有機(jī)酸代謝、水解酶活性和金屬離子、陰離子及核苷酸結(jié)合等影響顯著，可以進(jìn)一步挖掘感興趣的分子標(biāo)志(圖9和圖10)。具體到Unigenes，低鹽脅迫后，金烏賊O-晶體蛋白log2FC為-8.94，鈣調(diào)蛋白、膠原酶3、鈣依賴性蛋白激酶14、瞬時(shí)受體電位蛋白、Apolipophorins、促腎上腺皮質(zhì)激素釋放因子結(jié)合蛋白、含有Arrestin結(jié)構(gòu)域的odr-4蛋白、Ran結(jié)合蛋白9、XdhC、SFI1和磺基轉(zhuǎn)移酶等基因也得到了顯著表達(dá)。

KEGG數(shù)據(jù)庫中豐富的通路信息將有助于從系統(tǒng)水平去了解基因的生物學(xué)功能，針對(duì)兩兩分組的差異表達(dá)基因，可將差異基因顯示在KEGG 通路圖上[21]。本研究中，13 835條(2.7%)Unigenes注釋到KEGG數(shù)據(jù)庫中，PI3K-AKT信號(hào)通路、賴氨酸降解、cAMP信號(hào)通路、碳代謝和Rap1信號(hào)通路等值得進(jìn)一步挖掘和探討。根據(jù)KEGG代謝通路分類系統(tǒng)，氨基酸代謝、碳水化合物代謝、脂質(zhì)代謝、能量代謝、內(nèi)分泌系統(tǒng)、免疫系統(tǒng)、細(xì)胞生長和死亡相關(guān)基因具有重要的地位(圖5)。由圖11可見，低鹽脅迫對(duì)金烏賊類固醇生物合成、抗原加工與表達(dá)、脂肪消化與吸收、甘油脂質(zhì)代謝、花生四烯酸代謝、蛋白質(zhì)消化與吸收、酪氨酸代謝、雌激素信號(hào)通路和心肌細(xì)胞腎上腺素信號(hào)等通路影響顯著，可以進(jìn)一步做不同脅迫時(shí)間下以上通路的具體變化過程，增進(jìn)金烏賊對(duì)低鹽脅迫響應(yīng)的分子機(jī)制的了解。如對(duì)心肌細(xì)胞腎上腺素信號(hào)通路挖掘發(fā)現(xiàn)，低鹽度脅迫條件下，金烏賊肌鈣蛋白c、鈣調(diào)素和磷脂酶c基因顯著下調(diào)，而鈉/鉀運(yùn)輸酶亞單元β-2和幽閉緊結(jié)基因顯著上調(diào)；對(duì)IL-17信號(hào)通路挖掘發(fā)現(xiàn)，低鹽度脅迫條件下金烏賊腫瘤壞死因子α誘導(dǎo)蛋白3、膠原酶3和絲氨酸/蘇氨酸蛋白激酶基因顯著下調(diào)，而MMP3基因顯著上調(diào)。

本次試驗(yàn)得到了金烏賊轉(zhuǎn)錄組數(shù)據(jù)庫，獲取了與金烏賊脅迫有關(guān)的基因資源，豐富了金烏賊轉(zhuǎn)錄組數(shù)據(jù)庫，可為今后進(jìn)行金烏賊低鹽脅迫生理機(jī)制的探討、分子標(biāo)記的開發(fā)和關(guān)鍵基因的克隆等提供技術(shù)支撐。