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

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

（1.河北農業大學動物醫學院，河北 保定071001；2.河北省獸醫生物技術工程技術研究中心，河北 保定071001）

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

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

1 材料和方法

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

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

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

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

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

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

1.7 差異表達基因的Real-time FQ-PCR驗證 為了證實RNA-Seq數據統計的準確性，本研究中選擇7個上調基因（IL6、CCL20、CXCL4、CSF3、CCL4、IL8、IL1β）、3個下調基因（MHC2、FABP4、LDLR）進行Real-time FQ-PCR檢測。根據GenBank上發表的序列，應用軟件Primer 5.0設計各基因的特異性引物（含內參基因引物，見表1），提取試驗組和對照組細胞中的總RNA樣本，反轉錄成為cDNA后進行Real-time FQ-PCR檢測。反應體系為：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℃預變性5 min；95℃變性10 s，60℃延伸30 s，共進行40個循環。每個樣本進行3次重復檢測。最后獲得各組Ct值，按照2－△△Ct法計算目的基因的mRNA表達量，比較其與RNA-Seq數據分析所得的FC值的變化，并應用SPSS 22.0統計軟件對兩種方法的檢測結果進行相關性分析。

2 結果與分析

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

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

表2 總RNA樣本質量檢測結果Table 2 Quality test results of total RNA samples

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

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

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

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

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

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

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

3 討論與結論

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

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

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

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

[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].中國預防獸醫學報，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]宮曉華，李琦，李傳峰，等.鴨源新型鵝細小病毒DS15株生物學特性研究[J].中國禽業，2017，9（4）：52-55.

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

[7]李書光.小鵝瘟病毒保護性抗原的原核表達及免疫原性研究[D].北京：中國農業大學，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]王燕.新型鵝細小病毒VP3抗體特異性ELISA檢測方法的建立[D].保定：河北大學，2017.

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

[11]張清瑞.基于RNA-Seq的DHAV-1感染雛鴨的組織轉錄組分析[D].武漢：華中農業大學，2013.

[12]王亮.基于RNA-Seq技術分析小麥矮縮病毒侵染后小麥的基因表達差異[D].北京：中國農業科學院，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.