脊尾白蝦響應(yīng)環(huán)境脅迫的轉(zhuǎn)錄組學(xué)分析

王興強， 盧雪旎， 沈 曄， 曹 梅， 韋壽永， 李慶國， 茆丹婷， 秦傳新

(1.江蘇海洋大學(xué)海洋科學(xué)與水產(chǎn)學(xué)院，江蘇連云港 222005； 2.連云港旺島旅游開發(fā)有限公司，江蘇連云港 222000；3. 中國水產(chǎn)科學(xué)研究院南海水產(chǎn)研究所，廣東廣州 510300)

溫度、鹽度和溶解氧均是甲殼類生活環(huán)境中重要的環(huán)境因子。脊尾白蝦()具有耐低溫、低鹽和低氧等特點，可以在溫度2～35 ℃和鹽度4～35的環(huán)境條件下生存，時刻關(guān)注這些環(huán)境因子的變化在白蝦養(yǎng)殖中非常重要。低溫脅迫可使凡納濱對蝦()血淋巴溶菌和抗菌能力顯著降低，而低鹽和低氧脅迫會降低其抗氧化能力和免疫酶活性。斑節(jié)對蝦()體內(nèi)酚氧化酶活性在低溫脅迫下顯著下降，而C型凝集素表達量在低鹽脅迫下顯著上升。低氧脅迫時日本沼蝦()需要消耗大量糖原和磷酸精氨酸來維持能量需要，而高溫時其幼體也會產(chǎn)生明顯的應(yīng)激反應(yīng)。環(huán)境脅迫后的物種均可通過高通量轉(zhuǎn)錄組測序，經(jīng)生物信息學(xué)分析來發(fā)掘其響應(yīng)環(huán)境變化的調(diào)控基因。運用轉(zhuǎn)錄組測序技術(shù)，孫政等發(fā)現(xiàn)，脊尾白蝦體內(nèi)各種免疫相關(guān)基因表達上調(diào)以響應(yīng)各種病原刺激。梅園通過尼羅羅非魚()轉(zhuǎn)錄組測序，在溫度脅迫產(chǎn)生的差異表達基因中發(fā)現(xiàn)一種共有結(jié)構(gòu)域。本研究通過Illumina測序得到脊尾白蝦低溫、低鹽和低氧脅迫下的轉(zhuǎn)錄組數(shù)據(jù)，統(tǒng)一組裝，然后對不同脅迫條件下的差異表達基因進行注釋和富集分析，以期闡明脊尾白蝦在環(huán)境脅迫下的生理調(diào)控機制。

1 材料與方法

1.1 低鹽試驗

試驗所用脊尾白蝦購自連云港南極路水產(chǎn)品市場，在實驗室暫養(yǎng)7 d以適應(yīng)環(huán)境，暫養(yǎng)期間水體溫度24～26 ℃、鹽度31、pH值8.25以及水中溶解氧含量為7.6～8.1 mg/L。在低鹽試驗開始之際，隨機挑選出120尾脊尾白蝦[平均濕體質(zhì)量(4.12±0.56) g]用于試驗，設(shè)低鹽脅迫(鹽度0.2)、自然海水(鹽度31)2個處理，每個處理3個重復(fù)，共6個水族箱(水體約30 L)，每個水族箱分別投放脊尾白蝦20尾，然后對水體進行人工增氧。在進行低鹽脅迫之后，水族箱中的白蝦約在5.5 h之后蝦體顏色變化大致可描述為從尾部逐漸變白，在6 h之后蝦體整體呈現(xiàn)側(cè)臥狀，在這種狀態(tài)下根據(jù)經(jīng)驗預(yù)測被低鹽脅迫的脊尾白蝦已處于昏厥狀態(tài)。因此為了保證脊尾白蝦的取樣效果，對于低鹽脅迫組的脊尾白蝦，當其整體處于側(cè)臥狀態(tài)下5 min左右時開始進行取樣操作，分別在低鹽脅迫組與自然海水組的每個水族箱中隨機選取脊尾白蝦5尾，并使用液氮將整蝦快速冷凍，轉(zhuǎn)錄組測序備用。

1.2 低溫試驗

脊尾白蝦來源和暫養(yǎng)同“1.1”節(jié)，試驗開始時隨機挑選180尾脊尾白蝦，平均濕體質(zhì)量(1.27±0.15) g，設(shè)低溫脅迫(0 ℃)、自然水溫(25 ℃)2個溫度梯度處理，2個處理各設(shè)置3個重復(fù)，總計水族箱6個(水體30 L左右)，其中每個水族箱分別放置脊尾白蝦30尾，人工增氧，低溫驟變處理4 h。對照組、低溫脅迫組同時整蝦取樣，速凍于液氮中，轉(zhuǎn)錄組測序備用。

1.3 低氧試驗

脊尾白蝦來源和暫養(yǎng)同“1.1”節(jié)，試驗開始時隨機挑選60尾脊尾白蝦，平均濕體質(zhì)量(3.15±0.26) g，設(shè)低氧脅迫、自然溶解氧[溶解氧(7.59±0.66) mg/L]2個溶解氧梯度處理，2個處理各設(shè)置3個重復(fù)，總計水族箱6個(水體30 L左右)，每個水族箱分別放置脊尾白蝦10尾。用2 mm厚液體石蠟封閉低氧脅迫組水族箱，脅迫約4 h后白蝦處于昏厥狀態(tài)，溶解氧降低到(1.13±0.26) mg/L，整蝦取樣，速凍于液氮中，對照組同樣取樣保存，轉(zhuǎn)錄組測序備用。

1.4 RNA提取與文庫構(gòu)建

用液氮將凍存的樣品進行研磨，采取常規(guī)Trizol試劑法來提取總RNA，同時對它的純度、濃度和完整性進行檢測。樣品達到要求后，富集mRNA并以其為模板合成雙鏈cDNA鏈，最后進行PCR擴增，構(gòu)建測序文庫。

1.5 轉(zhuǎn)錄組測序與生物信息學(xué)分析

將樣品數(shù)據(jù)送至上海元莘生物醫(yī)藥科技有限公司進行高通量測序，評估原始測序數(shù)據(jù)和文庫構(gòu)建質(zhì)量，通過去除接頭序列、剪切低質(zhì)量讀段、舍棄N率(測序中識別不出的堿基率)達到10%讀數(shù)及長度小于25 bp片段的數(shù)據(jù)，最終得到高質(zhì)量測序數(shù)據(jù)(Clean data)。通過Trinity對所有高質(zhì)量數(shù)據(jù)進行從頭組裝，最終得到轉(zhuǎn)錄本和unigene。緊接著進行功能注釋，運用BlastX分別與NR、STRING、Swiss-Prot和KEGG數(shù)據(jù)庫進行比對，對基因的生物學(xué)功能進行進一步的探究。通過RSEM軟件對基因的表達量進行評估，通過FPKM表示基因的表達水平，并運用edgeR軟件進行基因的差異表達分析，主要富集到的功能和通路通過差異表達基因的GO富集分析和KEGG富集分析得到。

2 結(jié)果與分析

2.1 測序數(shù)據(jù)質(zhì)量分析

經(jīng)過濾后得到高質(zhì)量測序數(shù)據(jù)，低鹽、低溫和溶解氧處理組和對照組均產(chǎn)出>35M條高質(zhì)量序列，堿基數(shù)均>4 Gb，GC含量均>43%。測序質(zhì)量的好壞是用質(zhì)量值(Q，quality value)來評估的，Q30表示堿基測序出錯的概率為0.001，值越高出錯率越低。本次測序樣品的Q30堿基百分比均大于92%(表1)。

表1 不同處理數(shù)據(jù)評估統(tǒng)計

使用Trinity對高質(zhì)量數(shù)據(jù)進行序列組裝，將Contigs合并分組生成Graph后得到轉(zhuǎn)錄本，最終組裝得到轉(zhuǎn)錄本436 023條和unigene 305 682條，轉(zhuǎn)錄本和unigene的N50分別為814 bp和552 bp，其中序列長度為200～400 bp的unigene最多，占總數(shù)的67.79%，長度在1 kb以上的占8.11%(表2)。

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

2.2 unigene功能注釋

unigene序列是通過BLAST在NR、Swiss-Prot、GO、COG、KEGG各大數(shù)據(jù)庫中進行比對的，大部分unigene(56 898條)注釋到NR數(shù)據(jù)庫，31 695條注釋到KEGG數(shù)據(jù)庫， 28919 條注釋到GO數(shù)據(jù)庫。

COG功能分類結(jié)果見圖1，主要注釋到翻譯、核糖體結(jié)構(gòu)與生物生成(translation，ribosomal structure and biogenesis)，翻譯后修飾、蛋白質(zhì)周轉(zhuǎn)與伴侶蛋白(posttranslational modification，protein turnover，chaperones)，僅通用功能預(yù)測(general function prediction only)和信號轉(zhuǎn)導(dǎo)機制(signal transduction mechanisms)等功能中。

2.3 差異表達基因分析

差異表達基因中FDR<0.05并且log|FC|≥1的為差異顯著。低溫脅迫產(chǎn)生的差異表達基因中，具有顯著性表達差異的有264條，包括215條上調(diào)基因和49條下調(diào)基因，其中159條得到注釋。低鹽脅迫產(chǎn)生差異表達顯著的基因僅有1條，為顯著下調(diào)的α-葡萄糖苷酶(α-glucosidase)，低氧脅迫沒有產(chǎn)生差異表達顯著的基因。

2.3.1 GO注釋 低溫脅迫產(chǎn)生的264條顯著差異表達基因中，共有66條受到功能分類，GO注釋結(jié)果見圖2，其中，結(jié)構(gòu)分子活性中富集到52個上調(diào)基因。

2.3.2 GO富集分析 使用Goatools軟件對低溫脅迫產(chǎn)生的差異表達基因進行GO功能顯著性富集分析，為控制計算的假陽性率，對值進行校正，當校正后的<0.05時，GO功能顯著富集。<0.001的有表皮結(jié)構(gòu)組成(structural constituent of cuticle)、脂質(zhì)轉(zhuǎn)運活性(lipid transporter activity)和結(jié)構(gòu)分子活性(structural molecule activity)等分子功能，其中表皮結(jié)構(gòu)組成的顯著性較高，所富集到的差異表達基因也較多，為51條。<0.01的有肽基二肽酶活性(peptidyl dipeptidase activity)和有機物運輸(organic substance transport)等，<0.05的有單生物體運輸(single organism transport)、單有機體定位(single organism localization)和底物特異性轉(zhuǎn)運體活性(substrate specific transporter activity)等。同時對得到的功能節(jié)點標出編號和功能，以值進行顯著性排序，最顯著的10個節(jié)點信息見表3。

表3 GO富集(低溫脅迫)

2.3.3 KEGG注釋 使用KOBAS進行KEGG通路富集分析，經(jīng)過校正的值以0.05為閾值，定義差異表達基因中顯著富集的KEGG通路。低溫脅迫產(chǎn)生的差異表達基因顯著富集到腎素-血管緊張素系統(tǒng)(renin-angiotensin system，RAS)、美國錐蟲病(american trypanosomiasis)、阿米巴原蟲病(amoebiasis)、腎素分泌(renin secretion)、肥厚性心肌病(hypertrophic cardiomyopathy)、氨基糖和核苷酸糖代謝(amino and nucleotide sugar metabolism)、霍亂弧菌感染(vibrio cholerae)和補體與凝血級聯(lián)(complement cascades with clotting)等通路中(表4)。

表4 KEGG富集部分結(jié)果(低溫脅迫)

差異表達的基因在該通路的比例與所有注釋基因在該通路的比例的比值越大，校正后的值越小，通路越具有參考價值。由腎素-血管緊張素系統(tǒng)信號通路(圖3)可以看出，綠色邊框有3個，表示3個下調(diào)基因。由表5展示信息發(fā)現(xiàn)，3個下調(diào)基因均為血管緊張素轉(zhuǎn)換酶(angiotensin converting enzyme，簡稱ACE)基因，差異倍數(shù)對數(shù)值logFC分別為 -8.80、-6.65和-7.81。

表5 腎素-血管緊張素系統(tǒng)信號通路包含序列

3 討論與結(jié)論

王余菊獲取日本蟳()鰓組織轉(zhuǎn)錄組69.22 Gb高質(zhì)量測序讀數(shù)，Q30達93.46%以上，組裝得到90 973條轉(zhuǎn)錄本和52 972條unigene，N50值為2 450、1 775 bp，注釋率37.84%。董麗君等對凡納濱對蝦進行高通量測序，獲得50 921條unigene。N50為1 589 bp，注釋率37.28%。李喜蓮等進行紅螯螯蝦()轉(zhuǎn)錄組測序，獲得147 915 744條高質(zhì)量讀數(shù)，Q30達94.76%，組裝得到67 369條unigene，N50長度為 1 376 bp，注釋率30.83%。本次測序各樣品均產(chǎn)出>35×10條序列數(shù)，堿基數(shù)均>4 Gb，GC含量均>43%，Q30堿基百分比均>92%，因而測序堿基識別出錯率很低，測序質(zhì)量較好。組裝得到轉(zhuǎn)錄本436 023條和unigene 305 682 條，轉(zhuǎn)錄本和unigene的N50分別是814、552 bp。N50較短可能是本次測序數(shù)據(jù)量龐大，且200～400 bp長度范圍占比較多所導(dǎo)致。56 898條unigene注釋到NR數(shù)據(jù)庫，31 695條注釋到KEGG數(shù)據(jù)庫，28 919條注釋到GO數(shù)據(jù)庫，注釋率低的原因可能是在數(shù)據(jù)庫中近源物種基因注釋信息較少。

-葡萄糖苷酶又稱為葡萄糖基轉(zhuǎn)移酶(EC.3.2.1.20)，是糖苷水解酶家族的一員，水解碳水化合物中的-糖苷鍵生成葡萄糖。-葡萄糖苷酶具有雙重作用，即水解和轉(zhuǎn)糖苷，可以對-葡萄糖苷酶的活性進行抑制，從而降低生物體內(nèi)葡萄糖的生成速率，同時可以對生物體內(nèi)的脂肪水平進行調(diào)節(jié)。低鹽脅迫僅產(chǎn)生1個顯著下調(diào)的基因，為α-葡萄糖苷酶。低鹽脅迫下α-葡萄糖苷酶表達量顯著降低，可能導(dǎo)致脊尾白蝦機體難以維持正常的能量需求。

腎素-血管緊張素系統(tǒng)(RAS)是重要的血壓和水電解質(zhì)平衡調(diào)節(jié)系統(tǒng)，而血管緊張素轉(zhuǎn)移酶(ACE)是RAS中的一種重要的酶，又稱為肽酰二肽水解酶，ACE分為ACE 1(原稱ACE)和ACE 2，都是膜相關(guān)的鋅肽家族成員，可以使肽鏈的C-末端二肽殘基水解。血管緊張素Ⅱ受體1型(AT 1R)響應(yīng)循環(huán)中的肽類激素血管緊張素Ⅱ(Ang Ⅱ)的系統(tǒng)性變化，ACE負責(zé)將Ang Ⅰ轉(zhuǎn)換為AngⅡ，并通過調(diào)節(jié)血壓和電解質(zhì)穩(wěn)態(tài)來使緩激肽失活。Ang Ⅱ在RAS中起到關(guān)鍵作用，能增強血管平滑肌張力，升高血壓，還會造成組織損傷。AT 1R和ACE的特異性抑制劑對心血管系統(tǒng)有重要影響，血管緊張素轉(zhuǎn)換酶抑制劑(ACEⅠ)具有使血壓下降的作用。ACE 2水解AngⅡ生成血管緊張素 1-7[Ang(1-7)]，有舒張血管和抗炎的作用，抵消ACE、AngⅡ和AT 1R對腦神經(jīng)元的損傷作用。低溫脅迫產(chǎn)生的差異表達基因中，具有顯著性表達差異的有264條，包括215條上調(diào)基因和49條下調(diào)基因，159條得到注釋，其中66條受到GO功能注釋，結(jié)構(gòu)分子活性中富集到52個上調(diào)基因，推測與低溫脅迫相關(guān)。GO富集分析發(fā)現(xiàn)，低溫脅迫下最顯著的功能節(jié)點是脂質(zhì)轉(zhuǎn)運活性和肽基二肽酶活性，而肽基二肽酶活性中包含的基因為ACE。同樣，KEGG通路富集分析顯示，差異表達基因顯著富集到RAS中，且該通路富集顯著性最高，具有參考價值。由信號通路圖可以看出，該通路共包含3個表達下調(diào)的基因，均為ACE，推測低溫通過抑制血管緊張素轉(zhuǎn)換酶的表達，降低了血管緊張素Ⅱ的生成，從而降低血管平滑肌張力，抑制血壓的升高和組織損傷；同時動員脂肪和蛋白質(zhì)等能源物質(zhì)，以響應(yīng)低溫脅迫。