高亞硝酸鹽脅迫下日本囊對(duì)蝦肝胰腺的轉(zhuǎn)錄組分析

陳亭君,栗志民,袁 樂(lè),劉建勇,梁彩鳳

(廣東海洋大學(xué) 水產(chǎn)學(xué)院,廣東 湛江 524088)

日本囊對(duì)蝦(Marsupenaeus japonicus)是全球最具有經(jīng)濟(jì)價(jià)值的甲殼類動(dòng)物之一,廣泛分布于印度-西太平洋、日本、中國(guó)和澳大利亞等地區(qū)[1-4],養(yǎng)殖規(guī)模呈逐年升高態(tài)勢(shì),其養(yǎng)殖方式包括集約化和半集約化養(yǎng)殖[5]。國(guó)內(nèi)對(duì)該蝦類的養(yǎng)殖方式已逐漸由室外土池養(yǎng)殖轉(zhuǎn)變?yōu)槭覂?nèi)高集約化養(yǎng)殖,以滿足市場(chǎng)日益增長(zhǎng)的需求[6-7]。要提高集約化養(yǎng)殖的生產(chǎn)效率和經(jīng)濟(jì)效益,水質(zhì)控制格外重要[8]。近年來(lái),在集約化養(yǎng)殖中,由于過(guò)多的投餌量和不徹底的養(yǎng)殖池排污,使高含氮(N)排泄物和殘餌在水體中累積,超過(guò)了養(yǎng)殖水體中亞硝化細(xì)菌和硝化細(xì)菌的代謝能力,使得亞硝酸鹽氮的質(zhì)量濃度達(dá)到了較高水平[9-10]。

亞硝酸鹽氮是氨氧化成硝酸鹽的中間產(chǎn)物,它使氧合血藍(lán)蛋白轉(zhuǎn)化為脫氧血藍(lán)蛋白,降低血淋巴對(duì)氧的親和性,降低機(jī)體的輸氧能力,最終造成水產(chǎn)動(dòng)物缺氧甚至死亡[11-13]。大量研究已經(jīng)證實(shí)亞硝酸鹽氮對(duì)多種蝦類具有較強(qiáng)的毒性[14-18],包括損壞器官[19]、生長(zhǎng)和發(fā)育變慢[20-21]、降低存活率[22]和免疫能力[23-24]等。值得一提的是,鹽度下降會(huì)導(dǎo)致蝦類對(duì)亞硝酸鹽的耐受性降低[15,25],pH下降時(shí)亞硝酸鹽應(yīng)激會(huì)導(dǎo)致蝦類氮的排泄、離子調(diào)節(jié)和氣體交換被干擾,并可能導(dǎo)致載氧能力下降[26]。過(guò)去的研究中,發(fā)現(xiàn)亞硝酸鹽脅迫對(duì)日本囊對(duì)蝦的免疫和代謝均能產(chǎn)生影響,在亞硝酸鹽脅迫下,日本囊對(duì)蝦血淋巴中氧血色素苷、蛋白質(zhì),以及氧血色素苷/蛋白質(zhì)的比值減少,氮代謝和酸堿平衡發(fā)生改變,滲透壓降低,尿素增加[27-29]。血淋巴亞硝酸鹽和血淋巴尿素隨環(huán)境亞硝酸鹽和暴露時(shí)間的增加而增加[29]。此外,Zheng等[30]克隆了與凋亡相關(guān)的基因caspase-3和DAD-1,初步探討了亞硝酸鹽脅迫對(duì)免疫相關(guān)基因和凋亡相關(guān)蛋白的遺傳響應(yīng)的分子機(jī)制。

在甲殼類生物學(xué)領(lǐng)域,利用轉(zhuǎn)錄組測(cè)序技術(shù)(RNA-seq)研究基因表達(dá)已廣泛用于比較生理學(xué)、生態(tài)學(xué)、進(jìn)化、環(huán)境監(jiān)測(cè)和商業(yè)化養(yǎng)殖中[31]。近年來(lái),關(guān)于亞硝酸鹽脅迫下甲殼類轉(zhuǎn)錄組的研究?jī)H見于凡納濱對(duì)蝦(Litopenaeus vannamei)[32]和日本沼蝦(Macrobrachium nipponense)[33],研究顯示凡納濱對(duì)蝦和日本沼蝦在亞硝酸鹽脅迫下的免疫相關(guān)通路和凋亡通路非常活躍,已篩選得到許多與免疫反應(yīng)、解毒、凋亡途徑相關(guān)的候選基因。然而,有關(guān)亞硝酸鹽脅迫下日本囊對(duì)蝦的分子機(jī)制仍然知之甚少。甲殼類動(dòng)物的肝胰腺不僅是重要的消化器官,在脂質(zhì)、碳水化合物等代謝過(guò)程中起重要作用,而且是不可或缺的免疫器官,跟解毒和免疫息息相關(guān)[33-34]。因此,本研究通過(guò)轉(zhuǎn)錄組測(cè)序技術(shù)獲得在高亞硝酸鹽脅迫下日本囊對(duì)蝦肝胰腺轉(zhuǎn)錄組信息,為探討高亞硝酸鹽脅迫下的分子機(jī)制、豐富cDNA 數(shù)據(jù)庫(kù)的信息、識(shí)別免疫和凋亡等通路的差異基因提供分子證據(jù)。

1 材料和方法

1.1 實(shí)驗(yàn)材料

由廣東國(guó)聯(lián)水產(chǎn)有限公司提供第3代健康的120日齡混合家系日本囊對(duì)蝦,其平均體長(zhǎng)(49.28±4.79)mm,平均體重(1.39±0.38)g。實(shí)驗(yàn)前,將日本囊對(duì)蝦在水溫(28±0.2)℃、鹽度(29.8±0.2)、溶解氧(Dissolved Oxygen,DO)6.0 mg/L條件下于水泥池中馴化7 d,以減輕應(yīng)激反應(yīng)。

通過(guò)預(yù)實(shí)驗(yàn),確定高亞硝酸鹽脅迫質(zhì)量濃度為80 mg/L(在此質(zhì)量濃度下,脅迫96 h 的存活率約為80%)。實(shí)驗(yàn)分為對(duì)照組(CG)和高亞硝酸鹽組(N)。以新鮮海水作為對(duì)照(亞硝酸鹽質(zhì)量濃度＜0.02 mg/L)。采用分析純(NaNO2)溶于新鮮海水配制質(zhì)量濃度為2 000 mg/L母液,貯藏于陰暗干燥環(huán)境備用,實(shí)驗(yàn)時(shí)通過(guò)稀釋母液配制高亞硝酸鹽(80 mg/L)。

采用1 500 L塑料桶進(jìn)行對(duì)照組和高亞硝酸鹽組實(shí)驗(yàn),其他實(shí)驗(yàn)條件與馴化條件一致。實(shí)驗(yàn)期間不投餌,持續(xù)96 h,每24 h用高亞硝酸鹽試紙測(cè)定1次質(zhì)量濃度,以保持恒定水平。實(shí)驗(yàn)設(shè)置3個(gè)重復(fù)組。于6、12、24、48和96 h分別從高亞硝酸鹽組(N6、N12、N24、N48和N96)和對(duì)照組(CG6、CG12、CG24、CG48和CG96)取樣,每個(gè)時(shí)間點(diǎn)各取9尾蝦采集肝胰臟,保存于含1 m L RNAhold的離心管中。樣品在4℃下保存過(guò)夜,然后在-20 ℃下保存,直到提取RNA。

1.2 RNA提取及Illumina測(cè)序

利用TRIzol試劑(Invitrogen,US)從肝胰臟中提取總RNA,用1%瓊脂糖凝膠電泳檢測(cè)RNA 降解和污染狀況。分別通過(guò)NanoPhotometer?分光光度計(jì)(Implen,CA,USA)和Qubit?RNA Assay Kit and Qubit?2.0 熒光計(jì)(Life Technologies,CA,USA)檢查RNA 純度和濃度。采用生物分析儀2100系統(tǒng)(Agilent Technologies,CA,USA)中的RNA Nano 6000 檢測(cè)試劑盒評(píng)估RNA 完整性。本研究使用Illumina?的NEBNext?UltraTMRNA Library Prep Kit(NEB,USA)生成,共構(gòu)建了10個(gè)文庫(kù)。首先,利用帶有Oligo(d T)的磁珠從1μg總RNA 中富集有poly-A 尾的m RNA;然后,加入Fragmentation Buffer,將m RNA 隨機(jī)斷裂成200 bp 左右的小片段;第3 步,采用SuperScript Double-Stranded cDNA Synthesis Kit(Invitrogen)試劑盒,加入六堿基隨機(jī)引物(Illumina),以m RNA 為模板反轉(zhuǎn)錄合成第1鏈cDNA,進(jìn)行第2鏈合成,形成穩(wěn)定的雙鏈結(jié)構(gòu);第4步,雙鏈的cDNA 結(jié)構(gòu)為黏性末端,加入End Repair Mix將其補(bǔ)成平末端,隨后在3'末端加上1個(gè)A 堿基,用于連接Y 字形的接頭,具體步驟參見試劑盒說(shuō)明書;最后,Agencourt AMPure XP(Beckman Coulter,Brea,CA,USA)對(duì)PCR 產(chǎn)物進(jìn)行純化,并在Agilent 2100生物分析儀系統(tǒng)上對(duì)文庫(kù)質(zhì)量進(jìn)行評(píng)估。庫(kù)檢合格后,把不同文庫(kù)按照有效濃度及目標(biāo)下機(jī)數(shù)據(jù)量的需求混合(pooling)后,使用NovaSeq6000儀器進(jìn)行(Illumina,美國(guó))測(cè)序。

1.3 測(cè)序數(shù)據(jù)過(guò)濾和組裝

為了得到高質(zhì)量的測(cè)序數(shù)據(jù),必須將測(cè)序得到的原始測(cè)序序列(Sequenced Reads)或粗讀本(Raw Reads)過(guò)濾為凈讀本(Clean Reads):①去掉含測(cè)序接頭(Adapter)的讀數(shù);②去掉N(N 代表無(wú)法確定堿基信息)的比例＞10%的讀數(shù);③去除低質(zhì)量讀數(shù),即堿基質(zhì)量(Phred score,Qphred≤20的堿基數(shù)占整個(gè)堿基的50%以上的讀數(shù))。同時(shí),計(jì)算clean reads的Q＞20、30的堿基,以及G 和C的數(shù)量總和占總的堿基數(shù)量的百分比(Q20、Q30和GC含量)。所得到的高質(zhì)量clean reads用于后續(xù)分析,并采用Trinity軟件對(duì)claan reads進(jìn)行組裝[35]。

1.4 基因差異表達(dá)分析及功能注釋

首先,利用bowtie2軟件[36]將凈讀本比對(duì)到轉(zhuǎn)錄組序列上;然后,使用RSEM[37](http:∥deweylab.biostat.wisc.edu/rsem/)對(duì)bowtie2軟件的比對(duì)結(jié)果進(jìn)行統(tǒng)計(jì),進(jìn)一步得到每個(gè)樣品比對(duì)到每個(gè)基因上的read count數(shù)目,并對(duì)其進(jìn)行(Fragments Per Kilobase Million,FPKM)轉(zhuǎn)換[38]。先使用edgeR v.3.0.8軟件和1個(gè)尺度歸一化因子調(diào)整每個(gè)序列庫(kù)的讀計(jì)數(shù),而后采用TMM 對(duì)read count數(shù)據(jù)進(jìn)行標(biāo)準(zhǔn)化處理,再使用DEGseq對(duì)N 組和GC組之間的DEGs進(jìn)行差異分析,需q值(q-value)結(jié)合差異倍數(shù)(Fold Change,FC)來(lái)篩選,q≤0.05且|log2FC(Sample2/Sample1)|(基因表達(dá)量差異倍數(shù)是以2為底數(shù)的對(duì)數(shù)值)≥1,則該基因?yàn)轱@著差異表達(dá)基因[39]。

基于非冗余蛋白質(zhì)數(shù)據(jù)庫(kù)(Non-Redundant Protein Sequence Database,Nr)、非冗余核苷酸數(shù)據(jù)庫(kù)(Nucleotide Sequence Database,Nt)、蛋白質(zhì)和真核生物的同源群(eu Karyotic Ortholog Groups/Clusters of Orthologous Groups of Proteins,KOG/COG)、蛋白質(zhì)家族(Protein family,Pfam)、京都基因與基因組百科全書(Kyoto Encyclopedia of Genes and Genomes,KEGG)[40]和GO(Gene Ontology,GO)數(shù)據(jù)庫(kù)[41]對(duì)基因功能進(jìn)行注釋。

1.5 熒光定量驗(yàn)證轉(zhuǎn)錄組數(shù)據(jù)

隨機(jī)選擇9個(gè)DEGs,進(jìn)行qPCR 驗(yàn)證,分別是:ATP 結(jié)合盒轉(zhuǎn)運(yùn)蛋白(ATP-binding cassette transporters,ABC transporters)、二氫嘧啶脫氫酶(dihydropyrimidine dehydrogenase,DPD)、酚氧化酶原b(prophenoloxidase b,proPO b)、長(zhǎng)鏈脂酰輔酶A 脫氫酶(long-chain specific acyl-Co A dehydrogenase,LCAD)、細(xì)胞色素家族(cytochrome P450,family 2,subfamily J,CYP2J)、質(zhì)子偶聯(lián)氨基酸轉(zhuǎn)運(yùn)蛋白4(Proton-coupled Amino acid Transporter 4,PAT4)、甜菜堿同型半胱氨酸S-甲基轉(zhuǎn)移酶(betaine-homocysteine S-methyltransferase,BHMT)、C型凝集素(C-type lectin,CLEC)和磷酸烯醇丙酮酸羧基激酶(phosphoenolpyruvate carboxykinase,PEPCK),對(duì)上述9個(gè)DEGs進(jìn)行qPCR 驗(yàn)證。利用Primer 5.0軟件設(shè)計(jì)特異性引物(表1),送至生工生物工程(上海)股份有限公司合成。

qPCR 使用TransStart Tip Green Super Mix(北京全式金生物科技有限公司)試劑,以延伸因子1α(EF1α)為參考基因,通過(guò)羅氏LightCycler480 II實(shí)時(shí)熒光定量PCR 系統(tǒng)進(jìn)行擴(kuò)增。擴(kuò)增在384孔板上進(jìn)行,反應(yīng)總體積為10μL,包括:1μL cDNA、每個(gè)基因特異性正向和反向引物各0.2μL、5μL TransStart Tip Green qPCR Super Mix和3.6μL無(wú)酶水。qPCR 步驟為:94 ℃30 min;94 ℃5 s,60 ℃15 s,72 ℃10 s(45個(gè)循環(huán));95 ℃10 s,65 ℃60 s,95 ℃1 s。相對(duì)表達(dá)量采用2-ΔΔCT法計(jì)算,數(shù)據(jù)為平均值±標(biāo)準(zhǔn)差(Means±SD),通過(guò)SPSS19.0軟件中的單因素方差分析(one way ANOVA)進(jìn)行統(tǒng)計(jì)學(xué)檢驗(yàn),差異顯著性為P＜0.05,差異極顯著為P＜0.01。

2 結(jié)果

2.1 轉(zhuǎn)錄組的測(cè)序和從頭組裝

轉(zhuǎn)錄組測(cè)序后,從CG 組和NG 組構(gòu)建的10個(gè)文庫(kù)中共產(chǎn)生961 590 184個(gè)raw reads,除去包含適配器序列或poly-N 序列的讀取和原始讀取中的低質(zhì)量讀取,共獲得920 785 608個(gè)clean reads。在所有的文庫(kù)里,堿基質(zhì)量及組成分析顯示,各樣品Q30均≥93%,GC 含量為50.29%～52.95%(表2)。利用Trinity軟件對(duì)獲得的clean reads進(jìn)行組裝,去除冗余之后分別獲得74 856條轉(zhuǎn)錄本(Transcripts)和46 308條單基因簇(Universal Gene,Unigenes)。轉(zhuǎn)錄本的N50長(zhǎng)度(將轉(zhuǎn)錄本按照長(zhǎng)度從長(zhǎng)到短排序,依次累加轉(zhuǎn)錄本的長(zhǎng)度,當(dāng)累計(jì)轉(zhuǎn)錄本長(zhǎng)度達(dá)到轉(zhuǎn)錄本總長(zhǎng)的50%時(shí),拼接的轉(zhuǎn)錄本的長(zhǎng)度為N50,可用于評(píng)估拼接效果)和N90的長(zhǎng)度分別為2 408和470 bp,unigenes的N50和N90的長(zhǎng)度分別為1 833和435 bp。對(duì)組裝出來(lái)的unigenes進(jìn)行長(zhǎng)度分布統(tǒng)計(jì),其最小長(zhǎng)度為301 bp,分布在300～500 bp的有18 954條,占總數(shù)的40.93%,數(shù)量最多;大于2 000 bp只有6 323條,只占總數(shù)的13.65%,平均長(zhǎng)度為1 300 bp(表3和表4)。將Trinity軟件拼接得到的轉(zhuǎn)錄本序列,作為后續(xù)分析的參考序列。

表3 拼接長(zhǎng)度分布Table 3 Splicing length distribution

表4 長(zhǎng)度頻數(shù)分布Table 4 Length frequency distribution

采用單拷貝直向同源數(shù)據(jù)庫(kù)BUSCO 對(duì)拼接得到的unigenes進(jìn)行拼接質(zhì)量評(píng)估,結(jié)果顯示:有978個(gè)BUSCO 被完全覆蓋,完全匹配到的單拷貝(Complete and Single-Copy)的unigenes為902 條,占總數(shù)的92.2%;多拷貝(Complete and Duplicated)、部分片段匹配(Fragmented)和沒(méi)有匹配(Missing)的unigenes分別為28、26和22條,分別占總數(shù)的2.9%、2.7%和2.2%(表5)。

表5 拼接轉(zhuǎn)錄本BUSCO 評(píng)估Table 5 Busco evaluation of splicing transcripts

2.2 基因的功能注釋

將拼接得到的48 807條unigenes通過(guò)NR、NT、Swiss-Prot、PFAM、KO、KOG 和GO 七大數(shù)據(jù)庫(kù)進(jìn)行基因功能注釋,注釋到的unigenes 數(shù)量(比例)分別為:25 833(55.78%)、21 342(46.08%)、22 628(48.86%)、25 015(54.01%)、3 615(7.8%)、9 474(20.45%)和25 015(54.01%)條。在7個(gè)數(shù)據(jù)庫(kù)中至少注釋到1個(gè)數(shù)據(jù)庫(kù)的unigenes數(shù)量為34 361條,占總數(shù)的74.2%,而在7個(gè)數(shù)據(jù)庫(kù)都注釋到的unigenes數(shù)量為1 501條,占unigenes總數(shù)的3.24%(表6)。其中,比對(duì)到NR 數(shù)據(jù)庫(kù)的數(shù)據(jù)根據(jù)物種來(lái)源來(lái)分析,注釋到前3位的物種分別是美洲鉤蝦(Hyalella azteca),日本囊對(duì)蝦(Marsupenaeus japonicus)和凡納濱對(duì)蝦(Litopenaeus vannamei),匹配unigenes數(shù)量分別占總數(shù)的31.6%、14.4%和6.6%(圖1)。

圖1 NR 庫(kù)比對(duì)物種分布Fig.1 Comparable species distribution in the NR database

表6 基因注釋成功率統(tǒng)計(jì)Table 6 Statistics of success rate of gene annotation

將全部unigenes進(jìn)行GO 數(shù)據(jù)庫(kù)比對(duì),結(jié)果顯示:一共有25 015條unigenes被成功注釋和分類到生物過(guò)程(Biological Process)(25個(gè)亞類)、細(xì)胞組分(Cellular Component)(20個(gè)亞類)和分子功能(Molecular Function)(10個(gè)亞類)三大類中。在生物過(guò)程中,參與細(xì)胞過(guò)程(GO:0009987)、代謝過(guò)程(GO:0008152)、單細(xì)胞組織過(guò)程(GO:0044699)最多,分別注釋到14 335(57.31%)、13 080(52.29%)和12 325(49.27%)條unigenes;在細(xì)胞組分中,主要與細(xì)胞組分(GO:0044464)、細(xì)胞(GO:0005623)和膜(GO:0016020)有關(guān),分別注釋到6 724(26.88%)、6 724(26.88%)和5 159(20.62%)條unigenes;在分子功能中,主要與結(jié)合功能和催化活性有關(guān),分別注釋到11 937(47.72%)和11 516(46.04%)條unigenes(圖2)。

圖2 GO 注釋分類統(tǒng)計(jì)Fig.2 Classification statistics of GO annotation classification statistics

將unigenes與KOG 數(shù)據(jù)庫(kù)進(jìn)行比對(duì),結(jié)果顯示,9 474條unigenes被成功注釋并按26個(gè)KOG 進(jìn)行分類:①注釋到一般功能預(yù)測(cè)(General Function Prediction Only)的unigenes最多,為1 364 條,占總數(shù)的14.40%;②翻譯后修飾、蛋白質(zhì)轉(zhuǎn)化和分子伴侶(Posttranslational Modification,Protein Turnover,Chaperones)、信號(hào)轉(zhuǎn)導(dǎo)機(jī)制(Signal Transduction Mechanisms)、氨基酸運(yùn)輸和代謝(Amino Acid Transport and Metabolism)及翻譯、核糖體結(jié)構(gòu)與生物發(fā)生(Translation,Ribosomal Structure and Biogenesis)注釋到的unigenes分別有957(10.10%)、882(9.31%)、801(8.45%)和718(7.58%)條;③未知蛋白(Unamed Protein)注釋到的unigenes最少,僅占總數(shù)的0.02%(圖3)。

圖3 KOG 注釋分類統(tǒng)計(jì)Fig.3 Classification statistics of KOG annotation

對(duì)3 615條unigenes進(jìn)行KO 成功注釋后,根據(jù)unigenes參與的KEGG 代謝通路進(jìn)行分析,將其分為細(xì)胞過(guò)程(A)、環(huán)境信息處理(B)、遺傳信息處理(C)、代謝(D)和有機(jī)系統(tǒng)(E)五個(gè)分支,分別占總unigenes的13.67%、12.72%、16.74%、53.72%和23.82%。代謝通路的過(guò)程很多,其中富集的前3條通路分別為碳水化合物代謝(Carbohydrate Metabolism)、信號(hào)轉(zhuǎn)導(dǎo)(Signal Transduction)和氨基酸代謝(Amino Acid Metabolism),富集的unigenes數(shù)量分別為364(10.07%)、362(10.01%)和332(9.18%)條。另外,信號(hào)轉(zhuǎn)導(dǎo)(Signal Transduction)和免疫系統(tǒng)(Immune System)通路也被顯著富集(圖4)。

圖4 KEGG 代謝通路分類統(tǒng)計(jì)Fig.4 Classification statistics of KEGG metabolic pathway

2.3 差異基因表達(dá)分析

針對(duì)日本囊對(duì)蝦在5組不同處理時(shí)間進(jìn)行了N 組和CG 組兩兩比較分析,識(shí)別出亞硝酸鹽脅迫下發(fā)生變化的基因,并繪制火山圖(圖5)。結(jié)果顯示,日本囊對(duì)蝦在亞硝酸鹽協(xié)迫下N 組與CG 組兩兩相比共篩選出3 733個(gè)差異表達(dá)基因(Differentially Expressed Genes,DEGs);與對(duì)照組相比,N6、N12、N24、N48和N96組分別識(shí)別出593、606、1 089、497和988個(gè)DEGs(表7),其中差異基因數(shù)量最多的為N24(1 089個(gè)),最少的為N48(497個(gè))。總的來(lái)說(shuō),下調(diào)的DEGs數(shù)量(1 943個(gè))比上調(diào)的DEGs數(shù)量(1 830個(gè))多,在N12、N48和N96組,下調(diào)的DEGs比上調(diào)的多,而在N6和N24組則相反。

圖5 差異基因的火山圖Fig.5 Volcano map of differential genes

將N6 vs CG6、N12 vs CG12、N24 vs CG24、N48 vs CG48和N96 vs CG96各組的DEGs數(shù)量進(jìn)行統(tǒng)計(jì),繪成韋恩圖(Venn Diagram),找到在亞硝酸鹽脅迫下5 組不同處理時(shí)間的19 個(gè)共同DEGs(圖6)。把DEGs分別與NR、NT、KOG/COG、PFAM、Swiss-Prot、KEGG 和GO 七個(gè)數(shù)據(jù)庫(kù)進(jìn)行比對(duì)和功能注釋,至少注釋到1個(gè)數(shù)據(jù)庫(kù)的概率分別為89.88%、89.93%、90.36%、89.54%和90.49%(表7)。19 條共同的DEGs注釋和上下調(diào)情況如表8,結(jié)果顯示有8個(gè)假定蛋白基因(Hypothetical Protein),其他的多為免疫基因,如ATP結(jié)合盒轉(zhuǎn)動(dòng)體(ATP-binding cassette transporter)、C 型凝集素(C-type lectin)、磷酸烯醇丙酮酸羧激酶醇(phosphoenolpyruvate carboxykinase)、磷脂氫過(guò)氧化物(phospholipid-hydroperoxide glutathione peroxidase)、谷胱甘肽過(guò)氧化物酶(prophenoloxidase b)和預(yù)測(cè):細(xì)胞色素P450 9e2 亞型X1(cytochrome P450 9e2isoform X)等。

圖6 差異基因韋恩圖Fig.6 Venn diagram of DEGs

表7 DEGS的數(shù)目及注釋率Table 7 Quantity and annotation rate of DEGs

表8 不同脅迫時(shí)間下共同差異基因篩選情況Table 8 Screening of common differential genes under different stress times

2.4 差異基因的GO 功能分類和KEGG 富集分析

對(duì)N6 vs CG6、N12 vs CG12、N24 vs CG24、N48 vs CG48和N96 vs CG96五組比較所得的DEGs分別進(jìn)行GO 富集分析,共富集到1 758～2 399個(gè)GO terms type。DEGs富集到生物過(guò)程中的數(shù)量最多,占總數(shù)的58.76%～61.06%,分子功能次之(表9)。對(duì)5組顯著富集(P＜0.05)的GO 相關(guān)的上調(diào)基因(紅色)和下調(diào)基因(藍(lán)色)的分類統(tǒng)計(jì)做成柱狀圖(圖7),由圖7可見,N24 vs CG24組DEGs顯著富集到生物過(guò)程、細(xì)胞組分和分子功能,其他組富集到生物過(guò)程和分子功能;在生物過(guò)程中,幾丁質(zhì)代謝過(guò)程、含氨基葡萄糖的復(fù)合代謝過(guò)程和氨基糖代謝過(guò)程富集到所有組中;另外,氧化還原過(guò)程在N6 vs CG6、N12 vs CG12和N48 vs CG48組中顯著富集且富集到DEGs最多;在分子功能中,氧化還原酶活性、鐵離子結(jié)合、血紅素結(jié)合和幾丁質(zhì)結(jié)合富集于所有組中;在N24 vs CG24中富集到最多DEGs的3個(gè)細(xì)胞組分是細(xì)胞外區(qū)、線粒體和線粒體部分。

表9 差異基因GO 富集分析Table 9 GO enrichment analysis of differential genes

對(duì)N6 vs CG6、N12 vs CG12、N24 vs CG24、N48 vs CG48 和N96 vs CG96 的DEGs數(shù)量分別進(jìn)行KEGG 通路富集分析,結(jié)果顯示,分別有82、100、174、73和166個(gè)通路發(fā)生激活或抑制。對(duì)每組的前20個(gè)(top20)通路富集做散點(diǎn)圖(圖8)。如圖8所示,發(fā)生顯著富集的KEGG 通路有溶酶體(ko04142)、TGF-β信號(hào)通路(ko04350)、AMPK 信號(hào)通路(ko04152)、PI3K-Akt信號(hào)通路(ko04151)、p53信號(hào)通路(ko04115)、過(guò)氧化物酶體(ko04146)、吞噬體(ko04145)、細(xì)胞凋亡(ko04210)、PPAR 信號(hào)通路(ko03320)、膽堿代謝(ko05231)、亞油酸代謝(ko00591)、精氨酸和脯氨酸代謝(ko00330)等。

圖8 差異基因KEGG 通路富集散點(diǎn)圖Fig.8 Enrichment scatter plot of differential genes in KEGG pathway

2.5 qPCR驗(yàn)證RNA-seq

采用qPCR 檢測(cè)在高亞硝酸鹽脅迫下9 個(gè)差異基因(DPD、ABCH、Pro POb、ACADL、CYP2J、PAT4、BH MT、CLEC和PEPCK)的表達(dá)情況。qPCR 結(jié)果顯示,基因在不同時(shí)間的高亞酸鹽脅迫下,呈現(xiàn)顯著上調(diào)或下調(diào)(P＜0.05或P＜0.01),且與RNA-seq趨勢(shì)一致。這一結(jié)果進(jìn)一步驗(yàn)證了RNA-seq的可靠性和準(zhǔn)確性(圖9)。

圖9 9個(gè)差異表達(dá)基因的qRT-PCR 與轉(zhuǎn)錄組的比較Fig.9 Comparison of 9 DEGs by qRT-PCR and transcriptome

3 討論

日本囊對(duì)蝦集約化的養(yǎng)殖模式容易引起亞硝酸鹽的積累,從而影響該蝦類的健康養(yǎng)殖[14]。為了在分子水平上更好地了解蝦類對(duì)亞硝酸鹽脅迫反應(yīng),采用轉(zhuǎn)錄組測(cè)序研究高亞硝酸鹽脅迫下日本囊對(duì)蝦調(diào)控機(jī)制和差異基因表達(dá)。研究利用Illumina測(cè)序平臺(tái),獲得了不同時(shí)間(6、12、24、48和96 h)高亞硝酸鹽脅迫下日本囊對(duì)蝦肝胰腺的轉(zhuǎn)錄組數(shù)據(jù)。通過(guò)Trinity軟件對(duì)獲得的凈讀本進(jìn)行組裝,去除冗余之后獲得46 308條unigenes,N50和N90分別為1 833 bp和435 bp,平均長(zhǎng)度為1 098 bp。墨吉明對(duì)蝦(Fenneropenaeus mer-guiensis)的轉(zhuǎn)錄組獲得41 877條unigenes,N50為1 533 bp[42]。凡納濱對(duì)蝦的轉(zhuǎn)錄組獲得52 073條unigenes,平均長(zhǎng)度為520 bp,N50為745 bp[43]。日本囊對(duì)蝦組裝的unigenes數(shù)量比凡納濱對(duì)蝦少,而比墨吉明對(duì)蝦多,但組裝長(zhǎng)度大于凡納濱對(duì)蝦和墨吉明對(duì)蝦。采用BUSCO 對(duì)拼接得到的unigenes進(jìn)行質(zhì)量評(píng)估,結(jié)果顯示有978個(gè)BUSCO 被完全覆蓋,完全匹配到的unigenes為902條,占總數(shù)的95.1%,部分片段匹配和沒(méi)有匹配的unigenes分別占總數(shù)的2.7%和2.2%。因此,認(rèn)為該轉(zhuǎn)錄組具有高質(zhì)量的組裝拼接數(shù)據(jù)。

甲殼類動(dòng)物缺乏適應(yīng)性免疫系統(tǒng),而完全依賴于先天免疫系統(tǒng)抵制入侵的病原體或響應(yīng)環(huán)境脅迫[44-45]。甲殼類動(dòng)物免疫學(xué)研究主要集中在識(shí)別感染過(guò)程中激活的防御機(jī)制和生化途徑,如凝集素、酚氧化酶原系統(tǒng)(proPO)、吞噬和包圍等[46-47]。例如,中國(guó)對(duì)蝦被WSSV感染后,吞噬體、補(bǔ)體和凝血級(jí)聯(lián)反應(yīng)等與免疫應(yīng)答有關(guān)的通路以及免疫基因可被激活[48]。本研究對(duì)差異基因進(jìn)行KEGG 注釋,發(fā)現(xiàn)吞噬體(ko04145)通路在6、24和48 h時(shí)被富集,且該通路的差異基因顯著上調(diào),這表明在高亞硝酸鹽脅迫下通過(guò)激活吞噬體通路來(lái)參與免疫調(diào)節(jié)應(yīng)答,這與中國(guó)對(duì)蝦相似。另外,還發(fā)現(xiàn)大量與免疫相關(guān)其他的通路,如溶酶體(ko04142)、TGF-β信號(hào)通路(ko04350)、PI3K-Akt信號(hào)通路(ko04151)、p53信號(hào)通路(ko04115)、吞噬體(ko04145)、細(xì)胞凋亡(ko04210)、PPAR信號(hào)通路(ko03320)等。Guo等[32]研究了高亞硝酸鹽脅迫凡納濱對(duì)蝦的轉(zhuǎn)錄組同樣發(fā)現(xiàn)了凋亡信號(hào)通路、p53信號(hào)通路、PPAR 信號(hào)通路、MAPK 信號(hào)通路以及吞噬作用通路,這與本研究一致。而溶酶體(ko04142)、TGF-β信號(hào)通路(ko04350)在亞硝酸鹽脅迫日本沼蝦的轉(zhuǎn)錄組中也被富集到[33]。這些研究結(jié)果表明,亞硝酸鹽脅迫下,對(duì)蝦的這些免疫通路起著至關(guān)重要的作用,具體的機(jī)制還有待進(jìn)一步研究。

根據(jù)差異基因韋恩圖以及注釋結(jié)果,我們選擇了9個(gè)DEGs進(jìn)行qPCR分析,結(jié)果顯示高亞硝酸鹽脅迫后,這些參與免疫應(yīng)答的基因表達(dá)水平均發(fā)生顯著變化。Wei等[49]和Wang等[50]研究表明,C型凝集素在先天性免疫中起著重要的作用,能有效識(shí)別和消滅病原體。當(dāng)細(xì)菌感染凡納濱對(duì)蝦[48]和中國(guó)對(duì)蝦[50]時(shí),C型凝集素基因的表達(dá)水平會(huì)升高。在本研究中,C型凝集素基因的表達(dá)水平呈現(xiàn)上升-下降-上升的趨勢(shì),在早期(6～24 h)顯著上調(diào),到48 h時(shí)受到抑制,而96 h時(shí)又顯著上調(diào)。墨吉明對(duì)蝦的C型凝集素基因表達(dá)量在感染弧菌后的早期(12 h前)顯著上升,而在24和48 h時(shí)下降[51],與本研究類似。因此,C型凝集素基因在對(duì)蝦環(huán)境脅迫下起著重要的調(diào)節(jié)作用。甲殼類動(dòng)物的酚氧化酶原(pro PO)系統(tǒng)在非自我識(shí)別與宿主免疫反應(yīng)中起重要作用[52-54]。Prophenoloxidase b有助于甲殼類動(dòng)物血漿的黑色素化,為先天免疫系統(tǒng)的一個(gè)主要組成部分[55]。在本研究中發(fā)現(xiàn),Prophenoloxidase b的表達(dá)水平在亞硝酸鹽脅迫下早期顯著下調(diào),12 h之后顯著上調(diào),推測(cè)酚氧化酶原b(ProPO b)基因主要通過(guò)促進(jìn)表達(dá)來(lái)參與免疫應(yīng)答。同時(shí),在共同差異基因中還發(fā)現(xiàn)了與免疫相關(guān)的磷脂氫谷胱甘肽過(guò)氧化物酶(Phospholipid Hydroperoxide Glutathione Peroxidase,PHGPx)和細(xì)胞色素P450(cytochromeP450)異構(gòu)體基因。谷胱甘肽過(guò)氧化物酶(Glutathione Peroxidase,GPx)在脂質(zhì)和過(guò)氧化氫的解毒過(guò)程中起重要作用,脂質(zhì)和過(guò)氧化氫在吞噬或生理代謝過(guò)程中,隨著谷胱甘肽的氧化而迅速形成[56]。PHGPx是谷胱甘肽過(guò)氧化酶(GPx)家族中的一種抗氧化酶,它能降低磷脂的過(guò)氧化氫,維持生物膜的完整性[57]。已有研究表明,在LPS應(yīng)激下,擬穴青蟹(Scylla paramamosain)肝胰腺中的PHGPx表達(dá)水平分別在6和12 h顯著上調(diào),之后表達(dá)量逐漸下調(diào)至正常水平[58]。在本研究中發(fā)現(xiàn)PHGPx在亞硝酸鹽脅迫下表達(dá)量呈現(xiàn)上調(diào)和下調(diào)交替的現(xiàn)象,與擬穴青蟹中PHGPx的表達(dá)情況類似,推測(cè)PHGPx在日本囊對(duì)蝦的免疫調(diào)節(jié)中起著重要作用。細(xì)胞色素P450(CytochromeP450)是一種重要的解毒酶,在甲殼類動(dòng)物體內(nèi)外源性和內(nèi)源性化合物的生物轉(zhuǎn)化中起著重要作用[59-60]。在共同差異基因中發(fā)現(xiàn)一個(gè)細(xì)胞色素P450的異構(gòu)體基因,該基因的表達(dá)水平在高亞硝酸鹽脅迫6～96 h內(nèi)顯著上升。在凡納濱對(duì)蝦轉(zhuǎn)錄組中同樣發(fā)現(xiàn)了很多與P450相關(guān)的基因,這些發(fā)現(xiàn)為深入研究無(wú)脊椎動(dòng)物的解毒反應(yīng)提供了豐富的信息[31]。因此,推測(cè)這些免疫基因可能參與了亞硝酸鹽脅迫的免疫應(yīng)答,具體功能仍需要進(jìn)一步研究。