基于轉(zhuǎn)錄組測序技術(shù)研究三角魴和團(tuán)頭魴抗嗜水氣單胞菌感染基因表達(dá)差異

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

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

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

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

1 材料與方法

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

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

1.2 RNA提取及建庫測序

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

1.4 差異基因分析

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

2 結(jié)果與分析

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

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

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

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

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

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

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

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

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

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

注：不同顏色線表示基于DESeq分析的差異基因在不同時(shí)間點(diǎn)和樣本之間的表達(dá)變化。

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

2.4 差異表達(dá)基因的富集分析

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

3 討論

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

基于DESeq的基因表達(dá)趨勢聚類分析表明，三角魴和團(tuán)頭魴在3 h和24 h時(shí)間節(jié)點(diǎn)上的基因表達(dá)均表現(xiàn)為表達(dá)量無變化、表達(dá)上調(diào)和表達(dá)下調(diào)等3種情況，但其中主要以三角魴表達(dá)量無變化、團(tuán)頭魴表達(dá)下調(diào)(Cluster 4)，團(tuán)頭魴表達(dá)量無變化、三角魴表達(dá)下調(diào)(Cluster 1)為主?？梢姡囚?、團(tuán)頭魴之間差異表達(dá)的基因主要以三角魴或團(tuán)頭魴下調(diào)表達(dá)為主，這與羅非魚(Oreochromisniloticus)的情況類似[13]，但與蛋白質(zhì)組學(xué)分析結(jié)果存在差異[6]。此外，為從總體上評估三角魴和團(tuán)頭魴感染后基因表達(dá)差異，進(jìn)一步基于DESeq和DESeq2方法進(jìn)行了分析，并進(jìn)一步鑒定出上調(diào)表達(dá)基因FN1b和CFH、下調(diào)表達(dá)基因CFHL4和CCL4等，提示以上基因表達(dá)與三角魴和團(tuán)頭魴的感染差異相關(guān)。FN1b是編碼纖維連接蛋白(Fibronectin，F(xiàn)N)的兩種基因亞型之一[14]，纖維連接蛋白則是細(xì)胞外基質(zhì)的重要成分之一，而細(xì)胞外基質(zhì)被認(rèn)為是影響器官形成與修復(fù)的關(guān)鍵因子[15]。CFH(Complement factor H)即補(bǔ)體因子H ，調(diào)節(jié)血漿中補(bǔ)體的替代途徑，抑制補(bǔ)體活性并介導(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]。三角魴相對于團(tuán)頭魴，F(xiàn)N1b和CFH基因表達(dá)顯著上調(diào)，表明三角魴在受到細(xì)菌攻擊時(shí)更傾向于自我修復(fù)，而CCL4基因的下調(diào)，則表明團(tuán)頭魴調(diào)動(dòng)了更多的細(xì)胞因子用于消除細(xì)菌。CFHL4則是補(bǔ)體因子H相關(guān)蛋白(Complement factor H-related protein，F(xiàn)HRs)之一[16]，目前關(guān)于其具體機(jī)制并未明確，但越來越多的數(shù)據(jù)表明，F(xiàn)HRs與CFH具有相反的作用，即它們直接增強(qiáng)補(bǔ)體活性[20]。因此，推測三角魴相對于團(tuán)頭魴C(jī)FHL4基因表達(dá)的下調(diào)也與三角魴細(xì)胞自我修復(fù)基因有關(guān)。

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