陳亭君,栗志民,袁 樂,劉建勇,梁彩鳳

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

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

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

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

1 材料和方法

1.1 實驗材料

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

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

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

1.2 RNA提取及Illumina測序

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

1.3 測序數(shù)據(jù)過濾和組裝

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

1.4 基因差異表達分析及功能注釋

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

基于非冗余蛋白質(zhì)數(shù)據(jù)庫(Non-Redundant Protein Sequence Database,Nr)、非冗余核苷酸數(shù)據(jù)庫(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ù)庫[41]對基因功能進行注釋。

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

隨機選擇9個DEGs,進行qPCR 驗證,分別是:ATP 結(jié)合盒轉(zhuǎn)運蛋白(ATP-binding cassette transporters,ABC transporters)、二氫嘧啶脫氫酶(dihydropyrimidine dehydrogenase,DPD)、酚氧化酶原b(prophenoloxidase b,proPO b)、長鏈脂酰輔酶A 脫氫酶(long-chain specific acyl-Co A dehydrogenase,LCAD)、細胞色素家族(cytochrome P450,family 2,subfamily J,CYP2J)、質(zhì)子偶聯(lián)氨基酸轉(zhuǎ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),對上述9個DEGs進行qPCR 驗證。利用Primer 5.0軟件設(shè)計特異性引物(表1),送至生工生物工程(上海)股份有限公司合成。

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

2 結(jié)果

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

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

表3 拼接長度分布Table 3 Splicing length distribution

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

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

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

2.2 基因的功能注釋

將拼接得到的48 807條unigenes通過NR、NT、Swiss-Prot、PFAM、KO、KOG 和GO 七大數(shù)據(jù)庫進行基因功能注釋,注釋到的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個數(shù)據(jù)庫中至少注釋到1個數(shù)據(jù)庫的unigenes數(shù)量為34 361條,占總數(shù)的74.2%,而在7個數(shù)據(jù)庫都注釋到的unigenes數(shù)量為1 501條,占unigenes總數(shù)的3.24%(表6)。其中,比對到NR 數(shù)據(jù)庫的數(shù)據(jù)根據(jù)物種來源來分析,注釋到前3位的物種分別是美洲鉤蝦(Hyalella azteca),日本囊對蝦(Marsupenaeus japonicus)和凡納濱對蝦(Litopenaeus vannamei),匹配unigenes數(shù)量分別占總數(shù)的31.6%、14.4%和6.6%(圖1)。

圖1 NR 庫比對物種分布Fig.1 Comparable species distribution in the NR database

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

將全部unigenes進行GO 數(shù)據(jù)庫比對,結(jié)果顯示:一共有25 015條unigenes被成功注釋和分類到生物過程(Biological Process)(25個亞類)、細胞組分(Cellular Component)(20個亞類)和分子功能(Molecular Function)(10個亞類)三大類中。在生物過程中,參與細胞過程(GO:0009987)、代謝過程(GO:0008152)、單細胞組織過程(GO:0044699)最多,分別注釋到14 335(57.31%)、13 080(52.29%)和12 325(49.27%)條unigenes;在細胞組分中,主要與細胞組分(GO:0044464)、細胞(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)計Fig.2 Classification statistics of GO annotation classification statistics

將unigenes與KOG 數(shù)據(jù)庫進行比對,結(jié)果顯示,9 474條unigenes被成功注釋并按26個KOG 進行分類:①注釋到一般功能預(yù)測(General Function Prediction Only)的unigenes最多,為1 364 條,占總數(shù)的14.40%;②翻譯后修飾、蛋白質(zhì)轉(zhuǎn)化和分子伴侶(Posttranslational Modification,Protein Turnover,Chaperones)、信號轉(zhuǎn)導(dǎo)機制(Signal Transduction Mechanisms)、氨基酸運輸和代謝(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)計Fig.3 Classification statistics of KOG annotation

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

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

2.3 差異基因表達分析

針對日本囊對蝦在5組不同處理時間進行了N 組和CG 組兩兩比較分析,識別出亞硝酸鹽脅迫下發(fā)生變化的基因,并繪制火山圖(圖5)。結(jié)果顯示,日本囊對蝦在亞硝酸鹽協(xié)迫下N 組與CG 組兩兩相比共篩選出3 733個差異表達基因(Differentially Expressed Genes,DEGs);與對照組相比,N6、N12、N24、N48和N96組分別識別出593、606、1 089、497和988個DEGs(表7),其中差異基因數(shù)量最多的為N24(1 089個),最少的為N48(497個)?？偟膩碚f,下調(diào)的DEGs數(shù)量(1 943個)比上調(diào)的DEGs數(shù)量(1 830個)多,在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ù)量進行統(tǒng)計,繪成韋恩圖(Venn Diagram),找到在亞硝酸鹽脅迫下5 組不同處理時間的19 個共同DEGs(圖6)。把DEGs分別與NR、NT、KOG/COG、PFAM、Swiss-Prot、KEGG 和GO 七個數(shù)據(jù)庫進行比對和功能注釋,至少注釋到1個數(shù)據(jù)庫的概率分別為89.88%、89.93%、90.36%、89.54%和90.49%(表7)。19 條共同的DEGs注釋和上下調(diào)情況如表8,結(jié)果顯示有8個假定蛋白基因(Hypothetical Protein),其他的多為免疫基因,如ATP結(jié)合盒轉(zhuǎn)動體(ATP-binding cassette transporter)、C 型凝集素(C-type lectin)、磷酸烯醇丙酮酸羧激酶醇(phosphoenolpyruvate carboxykinase)、磷脂氫過氧化物(phospholipid-hydroperoxide glutathione peroxidase)、谷胱甘肽過氧化物酶(prophenoloxidase b)和預(yù)測:細胞色素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 不同脅迫時間下共同差異基因篩選情況Table 8 Screening of common differential genes under different stress times

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

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

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

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

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

2.5 qPCR驗證RNA-seq

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

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

3 討論

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

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

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