低鹽脅迫下金烏賊轉錄組學分析

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

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

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

1 材料與方法

1.1 試驗材料

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

1.2 低鹽脅迫試驗

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

1.3 RNA提取、cDNA文庫構建及原始測序數據獲取

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

1.4 原始測序數據處理及轉錄組從頭拼接

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

1.5 從頭組裝轉錄組功能注釋

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

1.6 基因表達量估計

利用Bowtie(http:∥bowtie-bio.sourceforge.net/index.shtml)將對照和低鹽脅迫處理組高質量的測序數據分別比對到組裝完成的轉錄組序列上，RSEM(http:∥deweylab.biostat.wisc.edu/rsem/)中利用Bowtie的比對結果進行表達量估計并進行FPKM轉換，從而得到基因的表達水平。基因差異表達分析使用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)

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

2 結 果

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

表1 測序數據統計

注：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.

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

表2 組裝結果統計

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

表3 金烏賊Unigenes注釋成功率

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

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

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

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

圖1 Unigenes GO二級分類統計

圖2 Unigenes GO二、三、四級分類統計

圖3 Unigenes COG分類統計

圖4 包含Unigenes數目最多的前20個通路

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

圖5 Unigenes KEGG注釋統計

對比到NR數據庫中的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條低鹽脅迫組的高質量測序數據分別比對到組裝得到的轉錄組序列上，定位百分比分別為24%和20%(表4)。

表4 與轉錄組比對統計

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

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

圖7 基因表達量分布

圖8 上下調基因GO注釋柱形

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

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

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

圖10 顯著性GO有向無環圖

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

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

(續表5)

3 討 論

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

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

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

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