感染GCRV的草魚miRNA測序與比較分析

肖調(diào)義 龍 哲 李 彬 張秋實(shí) 秦 玲 熊舒婷

(湖南農(nóng)業(yè)大學(xué)動(dòng)物科學(xué)技術(shù)學(xué)院， 長沙 410128)

草魚(Ctenopharyngodon idellus)是我國傳統(tǒng)的淡水養(yǎng)殖魚類之一， 2019年總產(chǎn)量達(dá)553.3×107kg，產(chǎn)量居全國第一[1]。但人工養(yǎng)殖條件下的草魚易感染草魚呼腸孤病毒(Grass carp reovirus， GCRV)引起草魚出血病， 嚴(yán)重制約了草魚產(chǎn)業(yè)的健康發(fā)展。GCRV疫苗是防治草魚出血病的重要手段之一， 但GCRV的多亞型與變異特性限制了疫苗的預(yù)防成效[2，3]。近年來逐步興起的抗病育種探索有了新的研究進(jìn)展[4]， 利用分子標(biāo)記輔助選育， 可大大提高選擇效率， 縮短育種年限， 已有一些水生動(dòng)物的抗病相關(guān)分子被標(biāo)記， 并利用該標(biāo)記培育成優(yōu)良品系[5，6]。生物機(jī)體的免疫水平與遺傳因素密切相關(guān)[7]， 汪亞平團(tuán)隊(duì)發(fā)現(xiàn)中國長江、珠江及黑龍江等不同水系草魚的GCRV抗病能力存在差異[8]， 因此基于草魚機(jī)體自身的GCRV抗性水平差異的特點(diǎn)開展高抗性草魚的選育是實(shí)現(xiàn)草魚抗病育種的又一重要策略。不僅是表達(dá)基因， miRNA也可列入草魚的主效抗性分子資源的篩選范圍[9，10]， 可望為高抗性草魚的分子輔助育種提供理論參考。

miRNA是一類長約22nt的內(nèi)源性的、進(jìn)化上高度保守的非編碼小RNA單鏈分子， 通過與靶基因的3′端非翻譯區(qū)(Untranslated regions， UTR)結(jié)合來抑制靶基因的翻譯或降解靶基因， 調(diào)控生物體的各項(xiàng)生命活動(dòng)[11]。miRNA調(diào)控生命體的大部分基因[12]，參與了生物體重要的生命活動(dòng)， 如個(gè)體發(fā)育、細(xì)胞增殖、分化、凋亡和免疫等[13—15]。病原體的入侵也會影響機(jī)體細(xì)胞內(nèi)miRNAs的表達(dá)水平， 借助高通量測序分析鑒定功能miRNA是目前比較普遍的方法。近年來， 一系列組學(xué)研究表明miRNAs是魚類抗病相關(guān)研究的重要候選功能分子， 參與硬骨魚類的抗病毒免疫反應(yīng)[16—18]。Zhang等[19]綜述了在一系列感染虹彩病毒rock bream iridovirus (RBIV)-C1的牙鲆(Paralichthysolivaceus)脾臟中篩選出了121個(gè)差異表達(dá)的miRNAs， 其預(yù)測靶基因主要富集在免疫應(yīng)答、信號轉(zhuǎn)導(dǎo)和細(xì)胞凋亡等信號通路中，揭示miRNA在病毒感染過程中參與了天然免疫過程。

草魚出血病是目前草魚產(chǎn)業(yè)的極大威脅， 草魚抗病力水平的高低與遺傳因素是緊密相關(guān)。因此，本研究借助高通量測序分析草魚感染前后的頭腎中的差異表達(dá)的miRNAs， 篩選出了一些候選功能miRNA分子， 并了解其在各免疫器官中的表達(dá)模式， 旨在為草魚的抗病毒分子機(jī)制的研究提供理論補(bǔ)充， 為高抗性草魚的選育及分子設(shè)計(jì)育種工作提供可參考的信息。

1 材料與方法

1.1 實(shí)驗(yàn)動(dòng)物

本實(shí)驗(yàn)的實(shí)驗(yàn)草魚均來自瀏陽烏龍漁場， 選取規(guī)格7—8 cm左右的6月齡、無病無傷、未感染過GCRV的當(dāng)年魚種進(jìn)行攻毒， 于室內(nèi)養(yǎng)殖系統(tǒng)中進(jìn)行攻毒， 養(yǎng)殖水溫為恒定28℃。

1.2 病毒毒株及GCRV侵染

草魚體內(nèi)攻毒實(shí)驗(yàn)所用病毒為草魚呼腸孤病毒Ⅱ型(GCRV-Ⅱ)， 長江所曾令兵實(shí)驗(yàn)室饋贈。取30尾草魚至養(yǎng)殖系統(tǒng)暫養(yǎng)2周(分2組， 對照組10尾，實(shí)驗(yàn)組20尾)。攻毒時(shí)， 實(shí)驗(yàn)組中每100 g草魚注射毒株1 mL， 對照組中的每尾草魚注射同等劑量的PBS緩沖液。

1.3 樣本采集及小RNA文庫構(gòu)建和測序

待攻毒后的草魚出現(xiàn)出血癥狀、瀕臨死亡時(shí)，分別采集對照組及處理組出血癥狀的草魚頭腎， 每組各采8份組織樣品， 迅速置于液氮速凍， 而后放置于-80℃。分別提取對照組及攻毒組的頭腎組織總RNA， 將每組中的8份RNA樣品等量混合， 分別構(gòu)建4個(gè)重復(fù)的miRNA測序文庫， 小RNA測序文庫制備采用TruSeq Small RNA Sample Prep Kits (Illumina公司， 美國)， 而后使用 Illumina Hiseq2000/2500對構(gòu)建好的文庫進(jìn)行測序， 測序讀長為單端 1X50 bp。隨后小RNA文庫的Solexa測序的工作由杭州聯(lián)川生物技術(shù)股份有限公司完成。

1.4 生物信息學(xué)分析

測序得到的原始數(shù)據(jù)為圖像文件經(jīng)過base calling， 得到以fastq格式存儲的結(jié)果文件， 稱為raw reads， 包含reads的序列及測序質(zhì)量信息， 為了保證信息分析的質(zhì)量， 必須對 raw reads 進(jìn)行質(zhì)控處理，使用FastQC軟件去除原始數(shù)據(jù)中有帶接頭的、低質(zhì)量的垃圾 reads， 得到clean reads， 其中5′接頭序列為: 5′-GTTCAGAGTTCTACAGTCCGACGATC-3′，3′接頭序列為: 5′-TGGAATTCTCGGGTGCCAA GG-3′; 利用ACGT101-miR軟件 (LC Sciences， Houston， Texas， USA)進(jìn)行后續(xù)的一系列分析， 先對clean reads進(jìn)行長度過濾， 篩選保留堿基長度在18—26 nt的序列(一般來說動(dòng)物的miRNA長度區(qū)間為18—26 nt)， 將篩選得到的數(shù)據(jù)與miRBase22.0數(shù)據(jù)庫進(jìn)行比對， 再與草魚的基因組(http://bioinfo.ihb.ac.cn/gcgd/php/index.php?page=download)進(jìn)行比對， 用來鑒定已知的miRNAs， 保守型的miRNA命名參照miRBase數(shù)據(jù)庫; 將沒有比對上已知miRNA的剩余序列與mRNA數(shù)據(jù)庫、Rfam數(shù)據(jù)庫(包含rRNA、tRNA、snRNA、snoRNA、miRNA等非編碼RNA)和Repbase數(shù)據(jù)庫進(jìn)行比對， 盡可能地發(fā)現(xiàn)并去除其中可能的rRNA、snoRNA、snRNA、tRNA及可能的repeat associate sRNA， 最終獲得的序列數(shù)據(jù)即有效的數(shù)據(jù); 將這些有效的數(shù)據(jù)再跟草魚基因組進(jìn)行比對， 過濾掉非基因組的序列， 而后將剩下的miRNA序列用RNAfold軟件進(jìn)行miRNA二級結(jié)構(gòu)的繪制和預(yù)測， 留下符合二級發(fā)卡結(jié)構(gòu)的miRNA序列， 定義為novel miRNA， 并這些novel miRNA用PC(Predicted Candidate)標(biāo)記， 注明5p 或 3p臂端位置; 接下來用ACGT101-miR軟件將所有鑒定到的miRNAs的表達(dá)量進(jìn)行歸一化處理(Norm值)[20]， 同時(shí)使用t檢驗(yàn)(t-test)的算法對這些miRNA進(jìn)行表達(dá)量的差異分析， 檢驗(yàn)篩選顯著性的參數(shù)為P_value，我們將P的閾值設(shè)定為0.01， 即當(dāng)Pvalue＜=0.01時(shí)，則表示該miRNA在兩組樣品間的表達(dá)存在顯著差異。根據(jù)差異miRNAs檢測結(jié)果， 使用R軟件中的pheatmap函數(shù)進(jìn)行層次聚類分析， 熱圖用顏色變化來反應(yīng)miRNA表達(dá)量的變化， 可根據(jù)顏色的變化，直觀地看出不同樣本中miRNA表達(dá)趨勢的變化。與此同時(shí)， 為了檢驗(yàn)對照組和處理組中4個(gè)樣品之間的重復(fù)性， 采用Pearson相關(guān)性分析將歸一化分析的數(shù)據(jù)進(jìn)行兩組樣品之間的重復(fù)性分析。將篩選出來的顯著差異的miRNAs(DE miRNAs)進(jìn)行靶基因預(yù)測， 使用的預(yù)測軟件為TargetScan和miRanda，并取2個(gè)軟件預(yù)測的交集。將預(yù)測得到的DE miRNAs的靶基因分別進(jìn)行GO (Gene Ontology)和KEGG (Kyoto Encyclopedia of Genes and Genomes)的富集分析。

1.5 miRNA靶基因功能分析

為挖掘miRNA的下游靶基因， 我們利用TargetScan和miRanda對差異顯著的miRNA分子進(jìn)行了靶基因預(yù)測。選取2款軟件預(yù)測出來的靶基因的交集結(jié)果， 作為差異miRNA的最終的預(yù)測靶基因，而后根據(jù)miRNA與靶基因之間的對應(yīng)關(guān)系， 首先找出了所有差異表達(dá)miRNA的靶基因mRNA與注釋庫中GO(Gene Ontology)對應(yīng)的數(shù)量， 計(jì)算得到的Pvalue≤0.05位閾值， 滿足此條件的功能定義為miRNA-mRNA關(guān)系對顯著富集的功能， 同理， 對靶基因進(jìn)行KEGG(Kyoto Encyclopedia of Genes and Genimics)的功能富集分析。

1.6 RNA提取、逆轉(zhuǎn)錄及熒光定量

RNA提取使用Trizol(Invitrogen公司， 美國)法，部分頭腎組織RNA用于高通量測序， 另一部分用于后續(xù)的實(shí)驗(yàn)驗(yàn)證， 逆轉(zhuǎn)錄使用miRNA 1ststrand cDNA synthesis kit(by stem loop; Vazyme公司， 中國)， miRNA逆轉(zhuǎn)錄采用Stem-loop RT引物方法， 同時(shí)用snRNA U6作為內(nèi)參; 熒光定量使用SYBR Green(Vazyme公司， 中國)作為染料進(jìn)行熒光定量PCR(quantitative real-time PCR)。利用qPCR對miRNA的表達(dá)進(jìn)行相對定量分析， 并利用2-ΔΔCt方法進(jìn)行miRNA的數(shù)據(jù)處理， qPCR的體系、程序方法及運(yùn)算參照之前的操作[13]。

1.7 數(shù)據(jù)的統(tǒng)計(jì)學(xué)分析

為了減小后續(xù)的候選功能分子的驗(yàn)證工作的難度， 進(jìn)行定量驗(yàn)證時(shí)， 我們縮小了差異基因的篩選范圍， 將t-test的檢驗(yàn)篩選顯著性的參數(shù)設(shè)置為P_value≤0.01， Foldchange＜-1和Foldchange＞1.5 log2(Foldchange)＜-1且log2(Foldchange)＞0.5|log2(Foldchange)|≥1且歸一化均值表達(dá)量(Norm值)≥500。所有的熒光定量的實(shí)驗(yàn)數(shù)據(jù)處理用graphpad7.1處理， 作圖和顯著性分析也均在此軟件中完成。定量相關(guān)實(shí)驗(yàn)結(jié)果的顯示均采用均數(shù)±標(biāo)準(zhǔn)差(mean±SD)，顯著性分析采用的方法為t檢驗(yàn)，*P＜0.05為差異顯著。

2 結(jié)果

2.1 miRNA測序數(shù)據(jù)

為挖掘GCRV應(yīng)答過程中的候選功能miRNAs，本實(shí)驗(yàn)分析了GCRV處理前后草魚頭腎中的miRNAs表達(dá)差異， 實(shí)驗(yàn)分兩組， 每組4個(gè)重復(fù)。利用 Hiseq 2000 技術(shù)平臺進(jìn)行了Solexa測序， 分別獲得了11165309和9946325條原始序列讀數(shù)(Raw Reads)。去除低質(zhì)量序列、接頭序列和短序列后，分別從兩個(gè)文庫中篩選出883407和733032條高質(zhì)量序列進(jìn)行下一步的分析。我們在對過濾后的Valid數(shù)據(jù)的總數(shù)(Total)及種數(shù)(Unique)進(jìn)行了長度分布統(tǒng)計(jì)， 結(jié)果顯示， 本次平均每個(gè)樣品獲得了長度約為18—26 nt的短鏈分子6918762條， 將其長度分布進(jìn)行統(tǒng)計(jì)分析后發(fā)現(xiàn)大部分的序列分子都分布在20—24 nt(占比91.76%)， 符合Dicer酶切割的典型特征， 其中， 22 nt的小RNA最為豐富， 占比達(dá)38.97%。在所有的短鏈RNA序列中， 22個(gè)核苷酸的小RNA最為豐富， 占比達(dá)38.97%(圖1A)， 這也說明了樣品測序結(jié)果的可靠性。同時(shí)為了檢測每個(gè)組內(nèi)的4個(gè)樣品的重復(fù)性， 我們檢測了8個(gè)樣品之間的組間關(guān)系， 組間關(guān)系結(jié)果顯示， 對照組和處理組的關(guān)系界線明確， 同組間的樣品相互之間的關(guān)系均大于99.4%， 說明了組間的重復(fù)好， 可信度高(圖1B)。

2.2 差異miRNA的篩選、鑒定及聚類分析

將miRNA與參考基因組比對， 從對照組和處理組中分別鑒定出703和739個(gè)成熟的miRNAs， 其中118個(gè)miRNAs僅在處理組中可以檢測到， 有82個(gè)miRNAs只在對照組中檢測到(圖1C)。通過數(shù)據(jù)庫比對， 我們檢測發(fā)現(xiàn)671個(gè)miRNA為novel miRNA。差異性顯著分析結(jié)果顯示， 當(dāng)P＜=0.01時(shí)， 獲得46個(gè)極顯著上調(diào)miRNA和88個(gè)極顯著下調(diào)miRNA(圖1D)。這表明在攻毒過程中， 一些miRNA的表達(dá)啟動(dòng)， 一些miRNA的表達(dá)被抑制， 預(yù)示miRNA在病毒侵染過程中發(fā)揮著各自不同的生物學(xué)功能。通過歸一化分析差異miRNA的表達(dá)量數(shù)據(jù)。根據(jù)樣品miRNA表達(dá)譜的相近程度， 我們將差異表達(dá)極顯著的miRNA進(jìn)行了聚類分析， 并采用lg (Norm值)進(jìn)行miRNA表達(dá)展示來直觀反映聚類表達(dá)模式。聚類分析顯示， 2個(gè)測序組中的4個(gè)平行樣品的重復(fù)性非常好(圖2)， 證實(shí)我們的采樣規(guī)范， 測序結(jié)果可信度高。

圖1 攻毒前后的草魚miRNA的測序結(jié)果分析Fig. 1 Comparing analysis of miRNA sequencing after GCRV challenge in grass carp

圖2 攻毒前后的草魚DE miRNA熱圖Fig. 2 Heat map of DE miRNA of after GCRV challenge in grass carp

本實(shí)驗(yàn)旨在篩選草魚抗GCRV的主效功能分子， 因此需盡可能篩選出差異倍數(shù)大、表達(dá)量相對高的miRNA。于是我們將攻毒前后的miRNA表達(dá)量的倍數(shù)變化值(FoldChange)和miRNA的歸一化均值表達(dá)量(Norm值)納入差異miRNA的篩選參數(shù)， 設(shè)置篩選參數(shù)為|log2(FoldChange)|＞=1、Norm值(CON組)＞500， 除去成熟序列和表達(dá)完全一致的重復(fù)miRNA后， 最終獲得3個(gè)極顯著上調(diào)miRNA和29個(gè)極顯著下調(diào)miRNA(表1)。

表1 差異表達(dá)miRNA列表Tab. 1 List of differentially expressed miRNAs (DE miRNA)

2.3 差異miRNAs的靶基因GO及Pathway富集性分析

在生物體內(nèi)， miRNA主要通過靶向抑制mRNA發(fā)揮其生物學(xué)功能， 我們將篩選得到的32個(gè)DE miRNAs進(jìn)行了靶基因預(yù)測， 同時(shí)將預(yù)測得到的所有靶基因進(jìn)行GO和KEGG的富集性分析。GO功能注釋分3大類， 分別是注釋基因的分子功能MF(Molecular function)、所處的細(xì)胞位置CC(Cellular component)及其參與到的生物過程BP(Biological process)。通過GO注釋的結(jié)果了解DE miRNA靶基因的功能分布。圖3A列出了排名前20的差異miRNA靶基因的GO富集， 富集的基因分別分布在分子功能MF中的轉(zhuǎn)移酶活性“Transferase activity”、共轉(zhuǎn)運(yùn)體活性“Symporter activity”、鈉通道調(diào)節(jié)器活性“Sodium channel regulator activity”、DNA結(jié)合轉(zhuǎn)錄因子活性“DNA-binding transcription factor activity”、激酶活性“Kinase activity”、 蛋白激酶抑制因子活性“Protein kinase inhibitor activity”、C-XC趨化因子受體活性“C-X-C chemokine receptor activity”， 生物過程BP中的體節(jié)發(fā)生過程“Somitogenesis”、多細(xì)胞生物發(fā)育過程“Multicellular organism development”、DNA轉(zhuǎn)錄調(diào)控過程“Regulation of transcription， DNA-templated”、凋亡調(diào)節(jié)過程“Regulation of apoptotic process”、磷酸化過程“Phosphorylation”、通過JAK-STAT途徑負(fù)調(diào)控受體信號通路“Negative regulation of receptor signaling pathway via JAK-STAT”、細(xì)胞因子調(diào)控的信號通路“Cytokine-mediated signaling pathway”及細(xì)胞位置CC大類中的線粒體呼吸鏈復(fù)合物Ⅱ， 琥珀酸脫氫酶復(fù)合物(泛素) “Mitochondrial respiratory chain complex Ⅱ， succinate dehydrogenase complex(ubiquinone)”、細(xì)胞膜“Membrane”、 閏盤“Intercalated disc”、細(xì)胞膜完整成分“Integral component of membrane”及內(nèi)質(zhì)網(wǎng)膜“Endoplasmic reticulum membrane”等。

圖3B顯示了排名前20的差異miRNAs靶基因富集通路， 包括纈氨酸、亮氨酸和異亮氨酸降解(泛素介導(dǎo)的蛋白水解， Valine， leucine and isoleucine degradation; Ubiquitin mediated proteolysis)、TGF-β信號通路(TGF-beta signaling pathway)、?；撬岷偷团；撬岽x(Taurine and hypotaurine metabolism)、RIG-I受體信號通路(RIG-I-like receptor signaling pathway)、丙酮酸鹽代謝(Pyruvate metabolism)、丙酮酸代謝(Propanoate metabolism)、PPAR信號通路(PPAR signaling pathway)、p53信號通路(p53 signaling pathway)、糖基磷脂酰肌醇通路(Glycosylphosphatidylinositol)、糖鞘脂生物合成(Glycosphingolipid biosynthesis-lacto and neolact)、糖酵解/糖異生(Glycolysis/Gluconeogenesis)、FoxO信號通路(FoxO signaling pathway)、葉酸生物合成(Folate biosynthesis)、丁酸甲酯代謝(Butanoate metabolism)、不飽和脂肪酸的生物合成(Biosynthesis of unsaturated fatty acids)、β-丙氨酸代謝(beta-Alanine metabolism)、上皮細(xì)胞的細(xì)菌入侵(Bacterial invasion of epithelial cells)、黃曲霉毒素生物合成(Aflatoxin biosynthesis)及脂肪細(xì)胞因子信號通路(Adipocytokine signaling pathway)。

圖3 差異miRNA靶基因的GO和KEGG富集Fig. 3 GO and KEGG enrichment of DE miRNA target genes

2.4 差異miRNAs的驗(yàn)證

為驗(yàn)證測序結(jié)果的準(zhǔn)確性， 我們根據(jù)對照組(control， CON)的歸一化數(shù)據(jù)Norm值的高低， 將這些DE miRNA分為低、中、高表達(dá)的三組， 分別在自行界定的低(500＜Norm＜600)、中(600＜Norm＜1000)和高(Norm ＞1000)3個(gè)表達(dá)量梯度中選擇了12個(gè)miRNA用于qPCR驗(yàn)證[miR-96、miR-194ap3_1ss23CT為低表達(dá)， miR-155、miR-722、miR-200家族成員(miR-200a/b/c、miR-141)為中表達(dá)量組， miR-122-5p、miR-194a、miR-100-5p和miR-21為高表達(dá)]。結(jié)果顯示這12個(gè)miRNA的表達(dá)趨勢與miRNA測序的結(jié)果是非常吻合的， miR-96、miR-194a-p3_1ss23CT、miR-722、miR-200家族成員(miR-200a/b/c、miR-141)、miR-122、miR-100-5p和miR-194a等在GCRV處理后顯著下調(diào)， miR-21、miR-155在GCRV處理后顯著上調(diào)(圖4)， 與測序結(jié)果非常吻合， 再次確認(rèn)miRNA測序數(shù)據(jù)的高準(zhǔn)確性和高可信度。值得注意的是， miR-21在GCRV感染組中的頭腎中表達(dá)量最高， miR-194a在兩個(gè)比較組中的差異倍數(shù)最大， 達(dá)到48倍之高。

2.5 DE miRNAs的組織表達(dá)分析

為篩選潛在的免疫相關(guān)的功能miRNA， 我們篩選了一些候選DE miRNA繼續(xù)比較分析其在各組織間的表達(dá)分布， 旨在探索候選miRNA在肝臟、脾臟、腎臟、頭腎、皮膚和鰓等免疫器官中的表達(dá)情況。我們挑選的差異顯著miRNA如下: miR-200家族成員在垂體中表達(dá)量最高， 其次在免疫組織腎臟、皮膚和鰓等組織中有較高表達(dá); miR-100在肌肉中表達(dá)量最高; miR-21在脾臟中表達(dá)量最高; miR-722在肝臟中表達(dá)量最高; miR-194a在腸道中表達(dá)量最高， 其次是肝臟和腎臟等組織中; miR-122在肝臟中表達(dá)量最高， 其次是脾臟和皮膚等(圖5)。

圖5 DE miRNA的組織分布Fig. 5 Expression pattern of DE miRNAs in tissues

3 討論

miRNA在硬骨魚類的病毒感染過程中發(fā)揮潛在重要功能， 許多miRNA會通過靶向作用免疫分子參與調(diào)控水產(chǎn)動(dòng)物的免疫過程[16]， 一些功能miRNA在病毒感染過程中可通過調(diào)控魚類的病毒識別受體及下游信號轉(zhuǎn)導(dǎo)中參與的免疫相關(guān)基因來及時(shí)地調(diào)控魚類的抗病毒天然免疫反應(yīng)。Najib等[21]在感染病毒性出血性敗血癥病毒(Viral hemorrhagic septicemia virus， VHSV)的牙鲆頭腎組織中鑒定了63個(gè)差異miRNA并預(yù)測了310個(gè)靶基因， 包括IL1B、IRF3、IRF5、IRF7、IL8R和Mx等天然免疫相關(guān)分子; 牙鲆在感染巨細(xì)胞病毒后， 被誘導(dǎo)上調(diào)的 miR-731可通過抑制 IRF7 的表達(dá)而減弱Ⅰ型 IFN 反應(yīng)， 從而有利于病毒自身的復(fù)制[22];在棒狀病毒感染鮸巨噬細(xì)胞過程中， 宿主細(xì)胞miR-3570 的表達(dá)顯著上調(diào)， 可靶向調(diào)節(jié)MAVS介導(dǎo)的 NF-κB 和 IRF3 信號通路， 抑制Ⅰ型 IFN 和抗病毒因子的產(chǎn)生從而促進(jìn)病毒快速增殖[23]; miR-144可靶向作用traf6直接破壞IRF7的轉(zhuǎn)錄， 抑制IFN反應(yīng)促進(jìn)病毒的復(fù)制[24]。為挖掘草魚GCRV侵染中的潛在miRNA-mRNA調(diào)控網(wǎng)絡(luò)， 我們對DE miRNA的預(yù)測靶基因進(jìn)行GO和KEGG的富集性分析。GO功能注釋顯示在分子功能MF中的DNA結(jié)合轉(zhuǎn)錄因子活性、激酶活性、蛋白激酶抑制因子活性及C-X-C趨化因子受體活性， 生物過程BP中的DNA轉(zhuǎn)錄調(diào)控過程、凋亡調(diào)節(jié)過程、磷酸化過程、通過JAK-STAT途徑負(fù)調(diào)控受體信號通路及細(xì)胞因子調(diào)控的信號通路等均與免疫調(diào)控過程緊密相關(guān)。與此同時(shí)， KEGG的功能富集結(jié)果發(fā)現(xiàn)包括纈氨酸、亮氨酸和異亮氨酸降解(泛素介導(dǎo)的蛋白水解)、TGF-β信號通路、RIG-I受體信號通路、PPAR信號通路、p53信號通路(p53 signaling pathway)及FoxO信號通路等免疫相關(guān)信號通路均發(fā)生明顯富集。因此， 下一步的研究中可將這些免疫相關(guān)通路中的靶基因進(jìn)行簡單的篩選和功能驗(yàn)證。

通過縮小遴選范圍， 我們獲得了29個(gè)下調(diào)及3個(gè)上調(diào)共32個(gè)顯著DE miRNA， qPCR結(jié)果顯示12個(gè)挑選的miRNA的表達(dá)均與測序結(jié)果相吻合(圖4)，不僅從實(shí)驗(yàn)上證實(shí)了測序結(jié)果的可靠性， 同時(shí)也預(yù)示了測序結(jié)果中篩選出來的這32個(gè)miRNA可能參與了GCRV應(yīng)答的免疫調(diào)控， 是潛在的功能候選miRNA， 提示我們在頭腎中表達(dá)量高和變化倍數(shù)差異比較大的DE miRNA均具有繼續(xù)研究的價(jià)值。在牙鲆、石斑魚(Epinephelusspp)、烏鱧(Channa argus)、虹鱒(Oncorhynchus mykiss)、大西洋鮭(Salmo salar)、大黃魚(Larimichthys crocea)、鮸(Miichthys miiuy)等不同魚類的不同病毒侵染表達(dá)譜分析中發(fā)現(xiàn)miRNA-21、miRNA-155和miRNA-100的表達(dá)均發(fā)生了顯著變化， 為保守的DE miRNAs， 同時(shí)它們均參與了高等脊椎動(dòng)物的免疫應(yīng)答調(diào)控[25]; 在細(xì)菌感染過程中， 鮸miR-21通過IRAK4介導(dǎo)的NF-κB信號通路調(diào)節(jié)炎癥反應(yīng)并且其在TLR信號通路調(diào)節(jié)中起重要作用[26]; Hua等[27]在草魚腎臟細(xì)胞系和斑馬魚(Danio rerio)成魚進(jìn)行miRNA的過表達(dá)， 提示了miR-155可通過靶向作用mIg的適應(yīng)性免疫中發(fā)揮作用; 蝦的miR-100通過靶向胰蛋白酶基因mRNA而發(fā)揮抗凋亡作用[15]; 在本次草魚GCRV感染下的差異分析中， 這3個(gè)miRNA都發(fā)生了顯著差異。綜上， miR-21、miR-100及miR-155可作為下一步研究對象的重要候選分子。

為了進(jìn)一步篩選免疫分子， 我們對部分DE miRNA進(jìn)行了組織表達(dá)分析， 旨在分析各DE miRNA在免疫相關(guān)組織中的表達(dá)情況。結(jié)果顯示miR-722和miR-122均在肝臟中表達(dá)量最高， 這與大西洋鮭中的結(jié)果是一致的[28]; miR-194a在腸道中表達(dá)量最高， 且miR-194a在頭腎中攻毒前后的表達(dá)差異達(dá)48倍之高， 牙鲆miR-194a可調(diào)控IRF7及IFNI等來調(diào)控牙鲆的免疫應(yīng)答[29]， 提示miR-194a也是一個(gè)潛在的免疫分子; miR-200家族成員在垂體中表達(dá)量最高， 與此同時(shí)其在免疫組織腎臟、皮膚和鰓等組織中有較高表達(dá)， 提示了miRNA家族成員在體內(nèi)生物學(xué)功能中的保守性; 草魚出血病的主要表現(xiàn)為鰓、皮膚、肌肉及內(nèi)臟充血， 而miR-200家族在這些組織中也有高表達(dá)， 并且均是GCRV應(yīng)答中的顯著差異分子， 預(yù)示著miR-200在草魚的GCRV應(yīng)答中發(fā)揮潛在重要功能。半滑舌鰨(Cynoglossus semilaevis)中的結(jié)果也顯示miR-200a和miR-200b可能在免疫調(diào)控中發(fā)揮功能[30]， 揭示了miRNA在種間的功能保守性。值得注意的是， 組織定量結(jié)果顯示miR-200家族均在垂體中有著超高表達(dá)， 這與此前我們在同是鯉科魚類的斑馬魚中的研究結(jié)果類似:miR-200家族成員(miR-200a/b/miR-429a)被證實(shí)在斑馬魚雌魚的生殖發(fā)育過程中發(fā)揮至關(guān)重要的功能[31，32]， 這些結(jié)果說明了miR-200表達(dá)特性在種間具有一定的保守性， 同時(shí)還提示我們miR-200不僅可作為免疫方向的研究內(nèi)容， 同時(shí)miR-200在草魚的生殖發(fā)育方面的功能可作為一個(gè)新的研究點(diǎn)。

綜上， 本研究利用高通量測序?qū)CRV處理前后的草魚頭腎組織進(jìn)行了miRNAs的差異分析， 并獲得了29個(gè)下調(diào)及3個(gè)上調(diào)共32個(gè)基因極顯著DE miRNAs， 通過進(jìn)一步的分析， 挖掘出了miR-21、miR-100、miR-155、miR-722、miR-122、miR-194a及miR-200等7個(gè)miRNA在免疫調(diào)控中發(fā)揮潛在重要功能， 為草魚抗病相關(guān)分子的研究及抗病草魚的選育工作提供研究思路。