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

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

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

1 材料與方法

1.1 試驗材料

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

1.2 低鹽脅迫試驗

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

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

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

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

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

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

利用Trinity軟件提供的ORF預測流程(http:∥trinityrnaseq.sourceforge.net/analysis/extract_proteins_from_trinity_transcripts.html)對組裝得到的所有轉(zhuǎn)錄本和Unigenes序列進行基因預測，運用HMMER軟件和Pfam數(shù)據(jù)庫(http:∥pfam.sanger.ac.uk/)對照，得到完整的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ù)庫進行比對，獲得Unigenes相應(yīng)的功能注釋信息。

1.6 基因表達量估計

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

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

差異表達基因的校驗公式：FDR<0.05

(1)

差異表達基因的校驗公式：log2|FC|1

(2)

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

2 結(jié) 果

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

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

注：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)計

注：取最長的轉(zhuǎn)錄本作為Unigene；N50.按照長度將組裝轉(zhuǎn)錄本/Unigenes由大到小排序，累加轉(zhuǎn)錄本/Unigenes的長度到總長度的一半時，對應(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ù)庫對Unigenes進行功能注釋，COG、KOG和NOG數(shù)據(jù)庫共注釋Unigenes 4179條，其中COG注釋2709條、KOG注釋3416條和NOG注釋156條(表3)。

表3 金烏賊Unigenes注釋成功率

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

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

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

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

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

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

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

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

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

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

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

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

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

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

圖7 基因表達量分布

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

圖9 差異表達基因GO富集柱狀圖

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

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

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

圖11 差異表達基因KEGG富集柱狀圖

表5 主要差異表達基因及所在信號通路

(續(xù)表5)

3 討 論

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

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

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

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