劉 凱，謝 楠，戴瑜來，阮松林，馮曉宇

(杭州市農(nóng)業(yè)科學(xué)研究院，浙江 杭州 310024)

魴屬(Magalobrame)魚類，在分類地位上屬于鯉形目(Cypriniformes)、鯉科(Cyprinidae)和鲌亞科(Cultrinae)，該屬主要包括4種魚類：三角魴(Megalobramaterminalis)、團頭魴(M.amblycephala)、廣東魴(M.hoffmanni)和厚頜魴(M.pellegrini)等。團頭魴，又名“武昌魚”，是我國特有的優(yōu)良草食性魚類，因食性廣、生長快、養(yǎng)殖成本低等優(yōu)點，其還是我國主要的大宗淡水水產(chǎn)養(yǎng)殖品種之一[1]。三角魴，在錢塘江附近地區(qū)被稱為“三角鳊”“塔鳊”，在黑龍江附近地區(qū)被稱為“法羅”。相對于團頭魴，三角魴的養(yǎng)殖區(qū)域較小，但因具有快速生長期長、商品個體大等特點，三角魴在浙江、湖南、山東、江蘇等地也具有一定的養(yǎng)殖規(guī)模，具有良好的養(yǎng)殖前景[2-3]。

近年來，團頭魴病害頻發(fā)，由嗜水氣單胞菌(Aeromonashydrophila)引起的細(xì)菌性敗血癥是團頭魴養(yǎng)殖過程中危害較大的疾病之一，嚴(yán)重影響?zhàn)B殖效益[4]。而養(yǎng)殖實踐發(fā)現(xiàn)，與同屬的團頭魴相比，三角魴的細(xì)菌性敗血癥發(fā)生率較低，即使發(fā)生了該病，使用常規(guī)藥物就可較好地控制，未有大面積發(fā)病的情況，且病死率低。三角魴和團頭魴之間差異小，親緣關(guān)系較近[5]，但2種魚類在感染嗜水氣單胞菌后所表現(xiàn)出的差異，除去環(huán)境因素外，遺傳因素是其中一個重要的原因。此前，方獻平等[6]采用非標(biāo)記定量蛋白質(zhì)組學(xué)技術(shù)分析了在嗜水氣單胞菌侵染脅迫后3 h、10 h和24 h肝組織蛋白質(zhì)組的變化，結(jié)果表明三角魴在戊糖磷酸途徑、脂肪酸代謝、亮氨酸降解等通路與團頭魴存在顯著差異。為在基因組注釋的豐富度以及期望于基因水平上了解三角魴和團頭魴在嗜水氣單胞菌侵染過程中的表達差異，本文基于轉(zhuǎn)錄組測序技術(shù)研究了三角魴和團頭魴在感染嗜水氣單胞菌后轉(zhuǎn)錄組的變化，以期了解三角魴和團頭魴在嗜水氣單胞菌感染上的差異性，為闡明嗜水氣單胞菌致病機制乃至抗病育種提供研究基礎(chǔ)。

1 材料與方法

1.1 實驗材料

實驗用三角魴和團頭魴由杭州市農(nóng)業(yè)科學(xué)研究院水產(chǎn)研究所提供，嗜水氣單胞菌株分離自江蘇武進某養(yǎng)殖場。致病實驗與方獻平等[6]報道實驗相同，選取了3 h和24 h兩個時間點，每個時間點分別取三角魴和團頭魴各5尾，體質(zhì)量為130.5～150.5 g。每5尾魚分別取肌肉、肝臟、血液3個組織并分別浸泡于3份RNAfixer無液氮RNA樣品儲存液中，于-80℃冰箱中保存?zhèn)溆?。測序時，每3個組織等質(zhì)量混合后作為1份實驗材料，進行建庫和測序。各實驗樣品被分別標(biāo)記為SJF3、SJF20、TTF3、TTF20，表示3 h和24 h兩個時間點采集的三角魴和團頭魴樣品。

1.2 RNA提取及建庫測序

RNA提取及建庫測序委托杭州英睿生物科技有限公司完成。取總RNA樣品，利用試劑盒合成cDNA，通過末端修復(fù)、連接接頭和純化，獲得三角魴和團頭魴的cDNA文庫，建好的文庫通過Illumina HiseqTM 2000平臺進行測序。測序樣本的原始數(shù)據(jù)(Raw Reads)過濾與組裝等步驟參照劉凱等[7]的方法進行，過濾數(shù)據(jù)時，采用Cutadapt去除3′端帶接頭的序列以及去除平均質(zhì)量分?jǐn)?shù)低于Q20的Reads，使用Trinity軟件對所有樣本的Clean Reads進行拼接，最終組裝獲得Unigenes，作為參考序列。

1.3 功能注釋、分類及代謝途徑分析

通過blastx、Blast2GO和KOBAS等程序?qū)nigenes比對到數(shù)據(jù)庫NR(NCBI non-redundant protein sequences)、COG(Cluster of orthologous groups of proteins)、GO(Gene ontology)和KEGG(Kyoto encyclopedia of genes and genomes)等，進而得到Unigenes的功能注釋信息。

1.4 差異基因分析

使用轉(zhuǎn)錄組表達定量軟件RSEM，以組裝獲得的Unigenes作為參考序列，分別將每個樣品的Clean Reads比對到Unigenes序列上，統(tǒng)計每個樣品比對到每一個基因上的Reads數(shù)，并計算每個基因的FPKM值，最終獲得基因表達譜。分別使用DESeq[8]和DESeq2[9]的R包對比分析三角魴和團頭魴在感染嗜水氣單胞菌后的差異基因表達，將P<0.05并且log2|(Fold change)|>1的基因定義為差異基因，分析示意圖見圖1。使用clusterProfiler[10]的R包通過GO和KEGG數(shù)據(jù)庫對主要的差異表達基因進行富集分析。

2 結(jié)果與分析

2.1 轉(zhuǎn)錄組數(shù)據(jù)與組裝

三角魴和團頭魴的樣本經(jīng)建庫、測序后，每個樣本獲得了超過4 G的數(shù)據(jù)量，共獲得Raw Reads 114 827 718條。去除接頭序列、低質(zhì)量序列后，共獲得Clean Reads 114 824 154條，占總Reads的99.996 8%(表1)。序列經(jīng)拼接和組裝后，共獲得62 780個Unigenes，平均長度達到531.883 bp，Unigenes的N50為652 bp，組裝完整性較高(表2)。

表1 測序數(shù)據(jù)質(zhì)量評估

表2 測序數(shù)據(jù)的統(tǒng)計匯總

2.2 轉(zhuǎn)錄組功能注釋

在拼接獲得Unigenes后，與公共數(shù)據(jù)庫進行比對注釋，62 780個Unigenes中有35 443個可注釋到NR數(shù)據(jù)庫，比例最高，為56.46%；僅有8 632個可注釋到COG數(shù)據(jù)庫，比例最低，為13.75%(表3)。在NR數(shù)據(jù)庫注釋中，注釋來源最多物種為斑馬魚(Daniorerio)，其次為鰱魚(Hypophthalmichthysmolitrix)、鳙魚(Hypophthalmichthysnobilis)等(圖2)。在GO數(shù)據(jù)庫注釋中，11 622個Unigenes被注釋到38個GO Term分類上，其中生物學(xué)過程GO Term分類的相關(guān)基因較多，如GO：0009987 細(xì)胞過程GO Term分類上注釋了4 599個Unigenes；分子功能GO Term分類中，GO：0005488 細(xì)胞結(jié)合GO Term分類相關(guān)基因最多，注釋了7 413個Unigenes。在COG數(shù)據(jù)庫注釋中，一般功能預(yù)測分類注釋的Unigenes最多，為1 979個；其次為翻譯后修飾、蛋白質(zhì)周轉(zhuǎn)、伴侶的分類，注釋Unigenes 1 353個；之后為轉(zhuǎn)錄分類，注釋Unigenes 684個。在KEGG數(shù)據(jù)庫注釋中，注釋Unigenes 最多的前3個通路分別是補體和凝血級聯(lián)、金黃色葡萄球菌感染、吞噬體，注釋Unigenes 達到9 898、3 721和1 653個。

表3 基因注釋結(jié)果統(tǒng)計

2.3 基因表達的差異分析

62 780個Unigenes中，有33 693個Unigenes在三角魴中表達，有32 936個Unigenes在團頭魴中表達，其中有24 834個Unigenes在三角魴和團頭魴中共表達?；贒ESeq，分別比較了三角魴、團頭魴在3 h和24 h時間點上基因表達的差異變化(圖3)。對TTF20-TTF3的差異分析表明，團頭魴感染24 h后與感染3 h后相比，顯著上調(diào)表達的基因有4個，分別為FN1b、PAH、CBLN14和INIH3b.2，顯著下調(diào)表達基因的前10個主要有ApolipoproteinA-I、C-typelectin14以及UBA、MHCclassⅠantigin和MHCclassⅡalphaantigen等。對SJF20-SJF3的差異分析表明，三角魴感染24 h后與感染3 h后相比，顯著上調(diào)表達的基因的前10個主要有Ovary-specificC1q-likefactor、ComplementC3以及bain-DAA*0101、Hymo-UA和Meam-UBA等主要組織相容性復(fù)合物(Major histocompatibility complexes，MHC)相關(guān)基因，顯著下調(diào)表達基因的前10個主要有CBLN6LOC793037、Phosphotransferase、TBL2和Cyca-UA1*01等。對SJF3-TTF3的差異分析表明，在感染3 h后，三角魴相對于團頭魴，顯著上調(diào)表達的基因的前10個主要有Preprohepcidin、FN1b、TBL2、Phosphotransferase和Cyca-UA1*01等，顯著下調(diào)表達基因的前10個主要有MHCclassⅡalphaantigen、HBB、HBBLOC113082649和Ba1globin等。對SJF20-TTF20的差異分析表明，在感染24 h后，三角魴相對于團頭魴，顯著上調(diào)表達的基因的前10個主要有ApolipoproteinA-I、Ovary-specificC1q-likefactor、Serotransferrin、Hymo-UA和bain-DAA*0101等，顯著下調(diào)表達基因的前10個主要有CFHL4、BHMTwu：fb53h01、G6PC1a.2和CBLN10RP71-1G18.9等。將以上顯著上調(diào)或下調(diào)表達的基因進行表達趨勢聚類分析(圖4)，三角魴、團頭魴在3 h和24 h時間點上基因表達趨勢可劃分為6類，其中聚類模式4所包含的差異表達基因最多，為16個；聚類模式1所包含的差異表達基因數(shù)量為15個；聚類模式6包含基因12個；聚類模式3包含基因10個；聚類模式2、5分別包含基因5個、4個。由此可見，三角魴、團頭魴基因表達趨勢聚類主要模式為4、1、6和3。從以上主要聚類模式可見，三角魴、團頭魴之間差異表達的基因主要以三角魴或團頭魴下調(diào)表達的較多，而上調(diào)表達的基因較少。

注：紅色點表示顯著上調(diào)的基因，其中顯著上調(diào)的前10個基因標(biāo)注了基因名。藍色點表示顯著下調(diào)的基因，其中顯著下調(diào)的前10個基因標(biāo)注了基因名。圖5同此。

注：不同顏色線表示基于DESeq分析的差異基因在不同時間點和樣本之間的表達變化。

基于DESeq和DESeq2，從總體上比較了三角魴相對于團頭魴基因表達的差異變化(圖5)?；贒ESeq對SJF-TTF的差異分析表明，顯著上調(diào)表達的基因的前10個主要有FN1b、CXCL12a、CFH、Transferrin(Fragment)和Serotransferrin等，顯著下調(diào)表達基因的前10個主要有CCL4、CFHL4、G6PC1a.2和CTSL等?；贒ESeq2對SJF-TTF的差異分析表明，顯著上調(diào)表達基因的前10個主要有Meam-UAA、TBL2、CFH、PAH、FN1b和Preprohepcidin等，顯著下調(diào)表達基因的前10個主要有Cyca-DAB1、CFHL4、C-typelectin14、G6PC1a.2和CCL4等。比較DESeq和DESeq2分析結(jié)果發(fā)現(xiàn)，有7個差異表達的基因在DESeq和DESeq2中均被鑒別出來，分別是上調(diào)表達基因FN1b和CFH，下調(diào)表達基因CCL4、CFHL4、G6PC1a.2，以及2個未注釋的基因g.Unigene024805和g.Unigene001582。

2.4 差異表達基因的富集分析

對基于DESeq和DESeq2方法分別獲得的前20個表達差異最顯著的基因進行GO富集分析，通過Rich factor、Pvalue值和富集到GO Term上的基因個數(shù)來衡量富集的程度(圖6)。其中，免疫反應(yīng)(GO：0006955)條目下富集的差異表達基因最多，30個GO注釋的差異基因中有5個差異基因在此條目下；其次為抗原遞呈(GO：0019882)和細(xì)胞膜(GO：0016020)2個條目，以及MHCⅡ類蛋白復(fù)合物(GO：0042613)和碳水化合物結(jié)合(GO：0030246)等條目。以上結(jié)果表明，三角魴和團頭魴之間的抗感染差異與抗原遞呈相關(guān)基因關(guān)系密切。進一步對以上獲得的差異表達基因進行KEGG富集分析，富集程度的表征與GO富集分析類似(圖7)。其中，自身免疫性甲狀腺疾病(ko05320)、Ⅰ型糖尿病(ko04940)、病毒性心肌炎(ko05416)、細(xì)胞黏附分子(ko04514)和抗原遞呈(ko04612)5個KEGG通路是富集程度最高的，18個KEGG注釋的差異基因中，均有5個基因富集在這些通路下。剩下的通路還包括IgA腸道免疫網(wǎng)絡(luò)(ko04672)、弓形蟲病(ko05145)和類風(fēng)濕性關(guān)節(jié)炎(ko05323)等。綜上，依然可見抗原遞呈相關(guān)基因與三角魴和團頭魴之間的抗感染差異有關(guān)，此外細(xì)胞黏附通路(ko04514)也與細(xì)菌感染密切相關(guān)。

3 討論

由于轉(zhuǎn)錄組測序無需事先以已知序列為目標(biāo)，因此可以對任意物種進行測序而無需知道物種的基因組背景，這為研究物種間基因表達差異提供了方便[11-12]。因此，本研究基于轉(zhuǎn)錄組測序分析三角魴和團頭魴抗嗜水氣單胞菌感染的基因表達的差異。將三角魴和團頭魴樣本的測序數(shù)據(jù)經(jīng)過濾后組裝Unigenes作為參考序列，使用RSEM軟件比對Unigenes序列確定表達量，以團頭魴為對照，分析三角魴Unigenes表達量的相對變化。此前，方獻平等[6]采用非標(biāo)記定量蛋白質(zhì)組學(xué)技術(shù)，分析在嗜水氣單胞菌侵染脅迫后3 h、10 h和24 h時三角魴和團頭魴肝組織蛋白質(zhì)組的變化。基于方獻平等[6]的蛋白應(yīng)答聚類分析發(fā)現(xiàn)，10 h和24 h的應(yīng)答蛋白質(zhì)優(yōu)選聚類后再與3 h的應(yīng)答蛋白質(zhì)聚類，表明10 h和24 h蛋白質(zhì)應(yīng)答具有相似性?；诖?，本研究選擇3 h和24 h作為時間節(jié)點開展轉(zhuǎn)錄組分析。由于方獻平等[6]的分析僅基于肝臟組織，不能完全反映整個機體對嗜水氣單胞菌侵染的應(yīng)答情況。因此，本研究進一步分析了肌肉、肝臟和血液的混合組織樣，以更全面了解三角魴和團頭魴在抗嗜水氣單胞菌侵染過程中基因表達的差異性。此外，由于DESeq2可以消除DESeq的邊界效應(yīng)，因此在進行整體分析中采用了DESeq2方法與DESeq方法共同進行了分析，減少單方法分析的誤差。

基于DESeq的基因表達趨勢聚類分析表明，三角魴和團頭魴在3 h和24 h時間節(jié)點上的基因表達均表現(xiàn)為表達量無變化、表達上調(diào)和表達下調(diào)等3種情況，但其中主要以三角魴表達量無變化、團頭魴表達下調(diào)(Cluster 4)，團頭魴表達量無變化、三角魴表達下調(diào)(Cluster 1)為主。可見，三角魴、團頭魴之間差異表達的基因主要以三角魴或團頭魴下調(diào)表達為主，這與羅非魚(Oreochromisniloticus)的情況類似[13]，但與蛋白質(zhì)組學(xué)分析結(jié)果存在差異[6]。此外，為從總體上評估三角魴和團頭魴感染后基因表達差異，進一步基于DESeq和DESeq2方法進行了分析，并進一步鑒定出上調(diào)表達基因FN1b和CFH、下調(diào)表達基因CFHL4和CCL4等，提示以上基因表達與三角魴和團頭魴的感染差異相關(guān)。FN1b是編碼纖維連接蛋白(Fibronectin，F(xiàn)N)的兩種基因亞型之一[14]，纖維連接蛋白則是細(xì)胞外基質(zhì)的重要成分之一，而細(xì)胞外基質(zhì)被認(rèn)為是影響器官形成與修復(fù)的關(guān)鍵因子[15]。CFH(Complement factor H)即補體因子H ，調(diào)節(jié)血漿中補體的替代途徑，抑制補體活性并介導(dǎo)細(xì)胞表面對替代途徑激活劑和非激活劑的區(qū)分[16]。CFH蛋白產(chǎn)生缺陷后可致替代途徑異常激活，導(dǎo)致細(xì)胞損傷[17]。CCL4[Chemokine(C-C motif)ligands 4 ]即CC型趨化因子配體4，是CC型趨化因子家族的一員，因其從脂多糖激活的巨噬細(xì)胞中被分離出來，亦稱MIP-1β[18]，它在免疫細(xì)胞遷移、活化、分化和增殖等方面具有重要作用[19]。三角魴相對于團頭魴，F(xiàn)N1b和CFH基因表達顯著上調(diào)，表明三角魴在受到細(xì)菌攻擊時更傾向于自我修復(fù)，而CCL4基因的下調(diào)，則表明團頭魴調(diào)動了更多的細(xì)胞因子用于消除細(xì)菌。CFHL4則是補體因子H相關(guān)蛋白(Complement factor H-related protein，F(xiàn)HRs)之一[16]，目前關(guān)于其具體機制并未明確，但越來越多的數(shù)據(jù)表明，F(xiàn)HRs與CFH具有相反的作用，即它們直接增強補體活性[20]。因此，推測三角魴相對于團頭魴CFHL4基因表達的下調(diào)也與三角魴細(xì)胞自我修復(fù)基因有關(guān)。

基于GO和KEGG富集分析，對DESeq和DESeq2方法鑒定出的主要差異基因進行了富集分析，均表明三角魴和團頭魴的感染差異與抗原遞呈相關(guān)基因密切相關(guān)?；谵D(zhuǎn)錄組學(xué)對羅非魚感染后的基因表達分析表明，吞噬作用在羅非魚的抗感染免疫中起到重要作用，其中主要組織相容性復(fù)合物基因就是主要的吞噬相關(guān)基因[13]。本文結(jié)果與此類似，也確定了MHC基因在感染后的重要作用，并表明其與三角魴和團頭魴抗感染差異相關(guān)?；趫F頭魴MHCⅡα基因的研究也表明，MHC基因與團頭魴抗病性密切相關(guān)[21]。在DESeq和DESeq2共同鑒定出的差異基因中，僅FN1b被KEGG注釋，其富集在ECM-受體相互作用(ko04512)通路(圖7中未顯示)，表明ECM即細(xì)胞外基質(zhì)與三角魴和團頭魴抗感染差異密切相關(guān)。此外，在GO和KEGG富集分析中均富集到CXCL12a基因，分別富集于免疫反應(yīng)(GO：0006955)條目和IgA腸道免疫網(wǎng)絡(luò)(ko04672)通路下。CXCL12[Chemokine(C-X-C motif)ligand 12]即CXC型趨化因子配體12，是一種小分子的細(xì)胞因子，與CCL4同屬于趨化因子家族，對淋巴細(xì)胞有強趨化作用[22]。而CXCL12基因在DESeq分析中，三角魴相對于團頭魴顯著上調(diào)，與團頭魴CCL4基因顯著上調(diào)不同，表明三角魴和團頭魴抗感染差異可能與招募淋巴細(xì)胞、促進細(xì)胞吞噬的機制不同相關(guān)。方獻平等[6]基于蛋白質(zhì)組學(xué)技術(shù)發(fā)現(xiàn)，相對于團頭魴，三角魴β-免疫球蛋白感染后顯著上調(diào)表達。但在本研究中，盡管檢測到了β-免疫球蛋白基因表達，但在三角魴和團頭魴間差異不顯著，這可能是由于蛋白應(yīng)答與基因表達的不同步性[23]，或存在其他原因，都有待進一步研究。