基于RNA-Seq篩選新型鴨細(xì)小病毒感染鴨胚成纖維細(xì)胞的差異表達(dá)基因

崔 元，李曉軒，孫繼國,2★，袁萬哲,2★

（1.河北農(nóng)業(yè)大學(xué)動物醫(yī)學(xué)院，河北 保定071001；2.河北省獸醫(yī)生物技術(shù)工程技術(shù)研究中心，河北 保定071001）

2014年以來，在中國大陸東南多省陸續(xù)發(fā)生了大規(guī)模的“鴨短喙-侏儒綜合征”（Beak atrophy and dwarfism syndrome，BADS），該病以鴨舌頭外伸腫脹、生長發(fā)育嚴(yán)重遲緩為主要特征，經(jīng)過病毒的分離、測序鑒定、遺傳進(jìn)化分析、動物模型建立等試驗操作，確定引發(fā)該病的病原為鴨細(xì)小病毒（duck parvovirus，DPV），也稱新型鴨細(xì)小病毒、新型鵝細(xì)小病毒等[1-4]。近年來，該病發(fā)病區(qū)域的不斷擴大給養(yǎng)鴨業(yè)造成了巨大的經(jīng)濟(jì)損失。目前，國內(nèi)外研究人員對水禽細(xì)小病毒病的病原展開了多方面的深入研究，并取得了很大進(jìn)展，目前的研究主要集中在病毒的生物學(xué)特性[5]、基因組結(jié)構(gòu)特征[6]、蛋白特性[7]、分子及免疫學(xué)檢測方法的建立[8-9]等方面，而極少有研究涉及到宿主感染病毒前后差異基因表達(dá)水平的變化。

轉(zhuǎn)錄組測序（RNA-Seq）是一種能夠用于研究宿主抗病毒感染應(yīng)答差異表達(dá)基因的技術(shù)，該技術(shù)能夠在單核苷酸水平針對特定物種的整體轉(zhuǎn)錄活動進(jìn)行檢測，從而全面快速地獲得該物種在某一狀態(tài)下的幾乎所有轉(zhuǎn)錄本信息，且測序結(jié)果的準(zhǔn)確度很高[10]，已被廣泛地應(yīng)用于篩選相關(guān)疾病基因。有學(xué)者利用RNA-Seq技術(shù)研究發(fā)現(xiàn)，吸血后的蚊子體內(nèi)有30%的轉(zhuǎn)錄本信息發(fā)生變化，為蚊媒病提供了控制依據(jù)[11]。有研究人員利用RNA-Seq技術(shù)研究了小麥矮縮病毒侵染小麥的復(fù)雜過程，為研究病毒致病的分子機制奠定了基礎(chǔ)[12]。本研究以DEF細(xì)胞為研究對象，采用RNA-Seq技術(shù)篩選DEF細(xì)胞感染DPV前后的差異表達(dá)基因，以期為揭示DPV感染宿主的分子機制，進(jìn)一步闡明DPV與宿主的相互作用，為DPV感染和發(fā)病機制提供重要的理論依據(jù)。

1 材料和方法

1.1 主要試驗材料RNA抽提試劑盒mirVana?miRNA ISOlation Kit購自Thermo Scientific公司，反轉(zhuǎn)錄試劑盒HiScriptⅡQ RT SuperMix for qPCR購自Vazyme公司，熒光定量PCR Quanti-Fast?SYBR?Green PCR Kit購自Qiagen公司。9～12日齡鴨胚購自保定某鴨場。DPV SD（KU516831）分離株鴨胚尿囊液由本試驗室分離保存。

1.2 鴨胚成纖維細(xì)胞的制備 取10～12日齡鴨胚，去除頭、四肢及內(nèi)臟；PBS清洗胚體后剪碎呈肉糜狀；加入適量胰酶置大離心管中；37℃水浴中消化30 min，離心去除上清；用PBS反復(fù)清洗、離心，直至上層液澄清；加入10%DMEM生長培養(yǎng)基，細(xì)胞篩過濾，計數(shù)；置細(xì)胞培養(yǎng)箱中培養(yǎng)，觀察細(xì)胞貼壁情況。

1.3 DPV感染DEF細(xì)胞模型 把DEF細(xì)胞以每孔為1×106～2×106個細(xì)胞接種到無菌的6孔塑料板上，將F5代DPV尿囊液經(jīng)0.22 μm濾器過濾除菌后，接種至已長成單層的DEF細(xì)胞，接種劑量為104EID50。置細(xì)胞培養(yǎng)箱培養(yǎng)144 h，收獲細(xì)胞培養(yǎng)物，盲傳5代，獲得的細(xì)胞培養(yǎng)物用PCR方法[2]進(jìn)行檢測。分別設(shè)置試驗組和對照組，試驗組是將PCR檢測為陽性的F5代細(xì)胞培養(yǎng)液接種至DEF細(xì)胞，對照組加等量的2%DMEM。兩組均于細(xì)胞培養(yǎng)箱中培養(yǎng)至72 h后，用適量的TRIzol消化細(xì)胞，于液氮或干冰中保存?zhèn)溆谩?/p>

1.4 測序文庫的構(gòu)建與測序 由上海歐易公司利用Illumina HiSeqTM 2500儀器進(jìn)行樣品的測序。主要步驟為：應(yīng)用試劑盒抽提試驗組和對照組樣品中的總RNA，經(jīng)質(zhì)檢合格后，富集純化mRNA；將mRNA片段化，合成一鏈cDNA，然后合成二鏈cDNA，并純化雙鏈cDNA；純化的雙鏈cDNA進(jìn)行末端修復(fù)、加A尾并連接序列接頭；選擇片段并利用PCR方法進(jìn)行DNA片段富集；構(gòu)建好的文庫經(jīng)Agilent 2100Bioanalyzer質(zhì)檢合格后，使用Illumina HiSeqTM 2500測序儀進(jìn)行測序。

1.5 數(shù)據(jù)預(yù)處理與質(zhì)量控制 將測序得到的原始數(shù)據(jù)（raw data）使用NGS QC Toolkit軟件進(jìn)行質(zhì)控（QC）并去除接頭，在此基礎(chǔ)上過濾掉低質(zhì)量堿基以及N堿基，最終得到適用于后續(xù)分析的高質(zhì)量的clean reads。由于鴨品種并無參考基因組序列，因此應(yīng)用Trinity軟件paired-end的拼接方法，將有overlap的clean reads拼接得到轉(zhuǎn)錄本序列，再利用TGICL軟件聚類去冗余延伸得到一套最終的Unigene，以此作為后續(xù)分析的參考序列。

1.6 差異表達(dá)基因的篩選與分析 通過blastx將Unigene序列分別與NR、SWISSPROT、KOG、GO和KEGG數(shù)據(jù)庫進(jìn)行比對，得到的基因注釋信息，以得到的Unigene的序列信息做庫，使用bowtie 2軟件通過mapping得到能比對到參考基因的clean reads的數(shù)量及占比，然后使用FPKM（fragments Per kb per Million reads）法求得各Unigene分別在試驗組和對照組樣本中的表達(dá)豐度，應(yīng)用DESeq軟件計算基因差異表達(dá)量，并進(jìn)行差異顯著性檢驗，根據(jù)log2FC≥1或≤-1且P≤0.05的條件，從中篩選出2組樣品中差異表達(dá)的基因，基于差異表達(dá)基因，進(jìn)行GO富集分析和KEGG信號通路富集分析，并結(jié)合注釋結(jié)果對其進(jìn)行描述。統(tǒng)計每個GO條目和KEGG條目中所包括的差異基因的個數(shù)，并用超幾何分布檢驗方法計算每個條目中差異基因富集的顯著性，當(dāng)P≤0.05時，認(rèn)為差異基因在該條目中顯著富集。

1.7 差異表達(dá)基因的Real-time FQ-PCR驗證 為了證實RNA-Seq數(shù)據(jù)統(tǒng)計的準(zhǔn)確性，本研究中選擇7個上調(diào)基因（IL6、CCL20、CXCL4、CSF3、CCL4、IL8、IL1β）、3個下調(diào)基因（MHC2、FABP4、LDLR）進(jìn)行Real-time FQ-PCR檢測。根據(jù)GenBank上發(fā)表的序列，應(yīng)用軟件Primer 5.0設(shè)計各基因的特異性引物（含內(nèi)參基因引物，見表1），提取試驗組和對照組細(xì)胞中的總RNA樣本，反轉(zhuǎn)錄成為cDNA后進(jìn)行Real-time FQ-PCR檢測。反應(yīng)體系為：2×Quanti-Fast?SYBR?Green PCR Master Mix，5 μL；上游引 物（10 μM），0.2 μL；下游引物（10 μM），0.2μL；cDNA，1μL；Nuclease-free H2O，3.6 μL。PCR擴增程序為：95℃預(yù)變性5 min；95℃變性10 s，60℃延伸30 s，共進(jìn)行40個循環(huán)。每個樣本進(jìn)行3次重復(fù)檢測。最后獲得各組Ct值，按照2－△△Ct法計算目的基因的mRNA表達(dá)量，比較其與RNA-Seq數(shù)據(jù)分析所得的FC值的變化，并應(yīng)用SPSS 22.0統(tǒng)計軟件對兩種方法的檢測結(jié)果進(jìn)行相關(guān)性分析。

2 結(jié)果與分析

2.1 總RNA樣本質(zhì)檢及測序質(zhì)量的評估 由表2可見，抽提的總RNA樣本的質(zhì)量滿足試驗要求（A類：RIN≥7且28S/18S≥0.7），RNA純度高、完整性好，可進(jìn)行后續(xù)試驗。提取樣品的總RNA后，對試驗組和對照組樣本進(jìn)行高通量RNA-Seq，分別獲得52 479 122和52 622 754條raw data，經(jīng)過過濾去除低質(zhì)量reads，獲得適合于后續(xù)數(shù)據(jù)分析的高質(zhì)量reads分別為50 475 368和50 881 866條，有效堿基百分比分別為96.1%和96.6%，說明測序得到的序列質(zhì)量較高。

表1 Real-time FQ-PCR引物序列表Table 1 Primers used in Real-time FQ-PCR

表2 總RNA樣本質(zhì)量檢測結(jié)果Table 2 Quality test results of total RNA samples

2.2 各樣本reads與Un igene的比對結(jié)果 使用bowtie 2軟件將預(yù)處理后的clean reads與得到的Unigene參考基因進(jìn)行比對，結(jié)果顯示，試驗組和對照組中能比對到參考基因的序列分別為44 239 788條（87.64%）和45 170 806條（88.77%）。

2.3 差異基因的篩選 使用FPKM法計算2組樣本中的基因表達(dá)量，以log2FC≥1或≤-1且P≤0.05為條件，篩選2組樣品的差異表達(dá)基因，結(jié)果顯示，差異表達(dá)基因共計4 073個，其中上調(diào)基因1 904個，下調(diào)基因2 169個。部分差異表達(dá)基因見表3。

表3 篩選的部分差異表達(dá)基因Table 3 Partial differently expressed genes

2.4 GO功能富集分析 對4 073個差異表達(dá)基因的生物學(xué)過程（BP）、細(xì)胞組分（CC）以及分子功能（MF）進(jìn)行GO功能富集分析。分析結(jié)果顯示，共有743個基因得到了共計3 275個GO功能注釋，其中有1 908個GO條目顯著富集差異基因（P≤0.05）。GO功能相關(guān)的差異表達(dá)基因中，試驗組與對照組相比共有375個基因上調(diào)表達(dá)，368個基因下調(diào)表達(dá)。其中有2 223個（67.88%）歸為BP，362個（11.05%）為CC，690個（21.07%）為MF，由此可見，大部分差異表達(dá)基因在BP中起到了關(guān)鍵作用。差異表達(dá)基因參與的BP主要包括炎癥反應(yīng)、膽固醇生物合成過程、正調(diào)節(jié)淋巴細(xì)胞趨化、趨化因子介導(dǎo)的信號通路、免疫反應(yīng)、對IFN-?的細(xì)胞應(yīng)答等；CC包括細(xì)胞外隙、肌鈣蛋白復(fù)合體、胞外區(qū)、細(xì)胞外基質(zhì)等；參與的MF主要表現(xiàn)在趨化因子活動、細(xì)胞因子活動、低密度脂蛋白粒子結(jié)合、神經(jīng)生長因子結(jié)合、CCR趨化因子受體結(jié)合等。

2.5 KEGG信號通路富集分析 對4 073個差異表達(dá)基因進(jìn)行KEGG信號通路富集分析，結(jié)果顯示，有368個差異基因得到了KEGG數(shù)據(jù)庫的功能注釋，共得到301個信號通路的功能注釋，其中在KEGG數(shù)據(jù)庫中顯著富集（P≤0.05）差異基因的信號通路條目有92個。富集到信號通路條目下的這368個差異表達(dá)基因中，試驗組與對照組相比共有200個基因上調(diào)表達(dá)，168個基因下調(diào)表達(dá)。在這301個KEGG注釋的信號通路中，試驗組與對照組相比上調(diào)的信號通路共243個，下調(diào)的共261個，有203個信號通路中既存在上調(diào)差異表達(dá)基因又存在下調(diào)基因。其中，差異基因顯著富集（P≤0.05）的上調(diào)信號通路條目共86個，主要包括細(xì)胞因子受體相互作用、IL-17信號通路、TNF信號通路、NF-κB信號通路、Jak-STAT信號通路、Toll樣受體信號通路等；差異基因顯著富集（P≤0.05）的下調(diào)信號通路條目共62個，主要包括類固醇生物合成、萜類化合物骨架生物合成、脂肪酸新陳代謝、補體和凝血級聯(lián)反應(yīng)、DNA復(fù)制、甘氨酸、絲氨酸和蘇氨酸代謝等。

2.6 Real-time FQ-PCR驗證差異表達(dá)基因 選取10個差異表達(dá)基因，并以β-actin為內(nèi)參基因，應(yīng)用Real-time FQ-PCR驗證RNA-Seq結(jié)果的準(zhǔn)確性。2組樣本中每個基因做3次重復(fù)檢測，將檢測結(jié)果進(jìn)行單因素方差分析，結(jié)果顯示，2組樣本中10個基因表達(dá)量水平的差異均極顯著（P≤0.01，見表4）。將2種方法檢測得到的基因表達(dá)水平FC值進(jìn)行相關(guān)性分析，得到Pearson相關(guān)系數(shù)r=0.751，P=0.012。驗證結(jié)果表明，每個基因的Real-time FQ-PCR檢測結(jié)果與RNA-Seq結(jié)果均具有相同的方向性，通過Real-time FQ-PCR檢測的10個基因表達(dá)水平的變化情況與通過RNA-Seq分析預(yù)測的變化（上調(diào)或下調(diào)）一致（圖1）。

表4 Real-time FQ-PCR驗證差異表達(dá)基因Table 4 Verifying of differently expressed genes using real-time FQ-PCR

3 討論與結(jié)論

有研究表明，感染DPV的病鴨體內(nèi)不同的器官可出現(xiàn)不同程度的水腫、充血、出血和壞死等病理變化，病理切片觀察顯示，病鴨的心、肝、腎、小腸可出現(xiàn)不同程度的炎性細(xì)胞浸潤，體內(nèi)各臟器（心、肝、脾、肺、腎、胸腺、胰腺、舌等）均可用PCR方法檢測出目的條帶，表明DPV對鴨的上述器官均可產(chǎn)生不同程度的致病作用和病理損傷。因此，一旦DPV感染機體，宿主將迅速做出反應(yīng)，開啟防御機制引發(fā)機體的免疫應(yīng)答，這必然會引起某些基因表達(dá)水平的變化，而分析研究這些差異基因?qū)τ谔骄考膊“l(fā)生發(fā)展規(guī)律及宿主抗病機制具有重要作用。目前篩選處理差異表達(dá)基因的技術(shù)主要有消減雜交法、基因芯片、轉(zhuǎn)錄組測序、定量或半定量PCR等。其中較為傳統(tǒng)的基因芯片技術(shù)雖然以其高通量、快速等優(yōu)點被廣泛應(yīng)用，但是其只能檢測已知的序列信息。本研究中涉及到的物種鴨目前并無可匹配的基因芯片，而RNA-Seq技術(shù)兼具高通量、檢測范圍廣、定量準(zhǔn)確、可重復(fù)性高等優(yōu)點，還能應(yīng)用于未知轉(zhuǎn)錄本發(fā)現(xiàn)、單核苷酸多態(tài)性等領(lǐng)域的研究，是基因水平研究領(lǐng)域的一個強大工具。本研究采用RNA-Seq技術(shù)篩選分析DEF細(xì)胞感染DPV后差異表達(dá)基因的數(shù)據(jù)結(jié)果，與Real-time FQ-PCR的驗證結(jié)果具有較高的符合率，說明RNA-Seq的測序結(jié)果能夠良好的反應(yīng)DEF細(xì)胞感染DPV后基因表達(dá)水平的總體變化情況。

圖1 Real-time FQ-PCR驗證差異表達(dá)基因結(jié)果Fig.1 Verifying of differently expressed genes using real-time FQ-PCR

RNA-Seq分析結(jié)果顯示，DEF細(xì)胞感染DPV后，與免疫、炎癥反應(yīng)相關(guān)的上調(diào)基因過表達(dá)，其中16個差異表達(dá)基因上調(diào)超過5倍，9個差異表達(dá)基因上調(diào)超過10倍，表明這些基因在DPV感染過程中發(fā)揮著重要作用。與白介素相關(guān)的基因如IL-1β、IL6和IL8等均有明顯的上調(diào)表達(dá)。IL-1β作為參與多種通路的關(guān)鍵的促炎細(xì)胞因子，通過募集免疫細(xì)胞參與炎癥反應(yīng)，與許多炎性疾病的病理生理學(xué)相關(guān)[13]。本研究中，趨化因子相關(guān)基因如CCL20、CXCL4和CCL4上調(diào)表達(dá)均超過10倍。其中，CXCL4上調(diào)超過13倍，CXCL4及其受體可參與免疫調(diào)節(jié)，能夠促進(jìn)單核巨噬細(xì)胞分泌大量趨化因子和細(xì)胞因子，可反饋作用于巨噬細(xì)胞抑制其凋亡并增強其吞噬能力[14]。此外，在與免疫和炎癥有關(guān)的顯著富集上調(diào)表達(dá)基因的信號通路中，參與細(xì)胞因子間相互作用的信號通路有33個基因，IL-17信號通路有17個基因，TNF信號通路有17個基因，NF-κB信號通路有16個基因，Jak-STAT信號通路有13個基因。

綜上所述，本研究首次應(yīng)用RNA-Seq技術(shù)評估DPV感染DEF細(xì)胞的基因表達(dá)變化，利用獲得的測序數(shù)據(jù)對DPV感染宿主后的免疫反應(yīng)、炎癥反應(yīng)等途徑的總體模式進(jìn)行了分析，為進(jìn)一步研究DPV的致病作用及宿主抗病機制奠定了基礎(chǔ)。但深入明確DPV的致病機制及相關(guān)基因的具體作用還需要進(jìn)一步的探索和試驗。

[1]Chen H,Dou Y,Tang Y,et al.Isolation and Genomic Characterization of a Duck-Origin GPV-Related Parvovirus from Cherry Valley Ducklings in China[J].PLoS One,2015,10(10):e0140284.

[2]崔元，王建昌，李玉保，等.鴨短喙-侏儒綜合征及其病原的初步研究[J].中國預(yù)防獸醫(yī)學(xué)報，2016，38（5）：348-351.

[3]Wang,S,Cheng,X,Chen,S,et al.Identification of a novel goose parvovirus(GPV)recombinant associated with short beak and dwarfism syndrome in mainland China,2015[J].Infect Genet Evol,2016,41:289-291.

[4]Li C,Li Q,Chen Z,et al.Novel duck parvovirus identified in Cherry Valley ducks(Anas platyrhynchos domesticus),China[J].Infection Genetics&Evolution,2016,44:278-280.

[5]宮曉華，李琦，李傳峰，等.鴨源新型鵝細(xì)小病毒DS15株生物學(xué)特性研究[J].中國禽業(yè)，2017，9（4）：52-55.

[6]傅秋玲，程龍飛，萬春和，等.半番鴨胚中鵝細(xì)小病毒的分離鑒定及其基因組分子特征分析[J].畜牧獸醫(yī)學(xué)報，2017，48（11）：2148-2156.

[7]李書光.小鵝瘟病毒保護(hù)性抗原的原核表達(dá)及免疫原性研究[D].北京：中國農(nóng)業(yè)大學(xué)，2016.

[8]Niu X,Chen H,Yang J,et al.Development of a TaqMan-based Real-Time PCR Assay for the Detection of Novel GPV[J].Journal of Virological Methods,2016,237:32-37.

[9]王燕.新型鵝細(xì)小病毒VP3抗體特異性ELISA檢測方法的建立[D].保定：河北大學(xué)，2017.

[10]張春蘭，秦孜娟，王桂芝，等.轉(zhuǎn)錄組與RNA-Seq技術(shù)[J].生物技術(shù)通報，2012，12：51-56.

[11]張清瑞.基于RNA-Seq的DHAV-1感染雛鴨的組織轉(zhuǎn)錄組分析[D].武漢：華中農(nóng)業(yè)大學(xué)，2013.

[12]王亮.基于RNA-Seq技術(shù)分析小麥矮縮病毒侵染后小麥的基因表達(dá)差異[D].北京：中國農(nóng)業(yè)科學(xué)院，2015.

[13]Arend W P,Palmer G,Gabay C.IL-1,IL-18,and IL-33 families of cytokines[J].Immunological Reviews,2008,223(1):20-38.

[14]Kasper B,Brandt E,Brandau S,et al.Platelet factor 4(CXC chemokine ligand 4)differentially regulates respiratory burst,survival,and cytokine expression of human monocytes by using distinct signaling pathways[J].Journal of Immunology,2007,179(4):2584.