隋維愷， 高 艷，楊茜茜，冷思竹，嚴(yán)慧娟，屈 羽，陳維剛*，杜小剛*

（1.四川農(nóng)業(yè)大學(xué)生命科學(xué)學(xué)院，四川雅安 625014；2.成都動(dòng)物園成都市野生動(dòng)物研究所，成都 610081）

豚鹿（Axis porcinus）隸屬于哺乳綱（Mammalia）偶蹄目（Artiodactyla）鹿科（Cervidae）斑鹿屬（Axis）。主要分布于云南西部靠近中緬邊境的耿馬、西盟兩縣的南丁河沿岸。多棲息于草地平原、沿海草地、森林、山地等地帶。1975年被《瀕危野生動(dòng)植物種國(guó)際貿(mào)易公約（CITES)》列入附錄Ⅰ[1]。2008 年，世界自然保護(hù)聯(lián)盟（IUCN）將豚鹿收錄到受威脅物種紅色名錄并列為瀕危（EN）等級(jí)[2]。目前國(guó)內(nèi)僅有北京動(dòng)物園、上海動(dòng)物園和成都動(dòng)物園三家機(jī)構(gòu)有圈養(yǎng)的豚鹿，其中，成都動(dòng)物園飼養(yǎng)數(shù)量占總數(shù)的80%以上[3]。豚鹿種群復(fù)壯在野生動(dòng)物保護(hù)中起到重要作用，但人工繁育的豚鹿幼崽存活率很低，主要原因是病原微生物感染引起肺炎和心肺功能衰竭而導(dǎo)致死亡[4]。疾病已經(jīng)嚴(yán)重影響豚鹿幼仔的存活，進(jìn)而嚴(yán)重影響豚鹿種群的發(fā)展[5]。隨著生物學(xué)技術(shù)的迅猛發(fā)展，運(yùn)用新的技術(shù)手段探究野生動(dòng)物的抗病毒機(jī)制已經(jīng)成為未來(lái)的發(fā)展趨勢(shì)[6]，為野生動(dòng)物疾病預(yù)防提供了新的研究思路。

雙鏈 RNA（double-stranded RNA，dsRNA）病毒是RNA病毒中的一大類群，其宿主范圍廣[7]。聚肌胞苷酸（Poly I：C）是人工合成的雙鏈RNA[8]。Poly I：C 能夠模擬病毒感染后所形成的雙鏈RNA，刺激宿主產(chǎn)生免疫反應(yīng)[9]。因此被廣泛應(yīng)用于病毒引起的宿主免疫應(yīng)答的研究中。外周血淋巴細(xì)胞具有易于分離，對(duì)動(dòng)物傷害小等特點(diǎn)常被用作基因表達(dá)分析的樣本來(lái)闡述免疫問(wèn)題，例如病毒-宿主交互、免疫性狀的生物標(biāo)記[10-11]等。

轉(zhuǎn)錄組測(cè)序（RNA-seq）是指通過(guò)新一代高通量測(cè)序技術(shù)，對(duì)樣品中的所有RNA及整個(gè)轉(zhuǎn)錄組進(jìn)行整體測(cè)序，計(jì)算出不同mRNA的表達(dá)量[12]。Illumina測(cè)序?qū)儆谛乱淮鷾y(cè)序技術(shù)適用于研究基因組不明確的物種的基因表達(dá)[13]，具有測(cè)序通量大、準(zhǔn)確性高、成本低等優(yōu)勢(shì)[14]。RNA-seq技術(shù)的大規(guī)模應(yīng)用，使得非模式生物在基因組層面基因表達(dá)變化的研究有了快速的發(fā)展。目前轉(zhuǎn)錄組測(cè)序技術(shù)已經(jīng)在野生動(dòng)物保護(hù)與利用等相關(guān)研究中被廣泛使用。如通過(guò)轉(zhuǎn)錄組測(cè)序技術(shù)揭示了大熊貓的微衛(wèi)星位點(diǎn)多態(tài)性，篩選出與狀態(tài)特有生理相關(guān)的潛在分子標(biāo)記[15]；利用轉(zhuǎn)錄組測(cè)序技術(shù)分析了嗜水氣單胞菌感染后齊口裂腹魚(yú)免疫相關(guān)基因的表達(dá)情況[16]。目前利用轉(zhuǎn)錄組測(cè)序技術(shù)分析豚鹿的抗病毒免疫機(jī)制的研究尚未見(jiàn)報(bào)道。

本研究利用轉(zhuǎn)錄組測(cè)序技術(shù)對(duì)Poly I：C刺激后的豚鹿外周血淋巴細(xì)胞進(jìn)行測(cè)序分析，發(fā)現(xiàn)Poly I：C能夠引起基因的差異表達(dá)。對(duì)差異基因進(jìn)行富集分析，進(jìn)一步篩選出免疫相關(guān)基因，并使用qPCR方法對(duì)轉(zhuǎn)錄組數(shù)據(jù)的準(zhǔn)確性進(jìn)行了驗(yàn)證，為后續(xù)豚鹿抗病毒免疫機(jī)制的研究奠定了基礎(chǔ)。

1 材料和方法

1.1 材料與試劑

Ficoll Plus細(xì)胞分離液，北京索萊寶科技有限公司；PolyI：C，上海西格瑪奧德里奇有限公司；Trizol RNA提取試劑，美國(guó)Thermo公司；primescriptTMRT reagent kit反轉(zhuǎn)錄試劑，大連寶生物有限公司；SYBRgreen II熒光染料，美國(guó)Bio-rad公司；試驗(yàn)所用引物由北京擎科有限公司合成。

1.2 儀器與設(shè)備

Nanodrop-2000超微量核酸蛋白測(cè)定儀，美國(guó)Thermo Scientific公司；CFX96型熒光定量PCR儀，美國(guó)Bio-Rad公司；ChemiDocTM XRS+型凝膠成像系統(tǒng)，美國(guó)Bio-Rad公司。

1.3 試驗(yàn)方法

1.3.1 樣品采集和處理

6只成都動(dòng)物園內(nèi)圈養(yǎng)的雄性豚鹿，分別無(wú)菌收集抗凝血10 mL。按照Ficoll Plus細(xì)胞分離液說(shuō)明書(shū)分離淋巴細(xì)胞，分離后的細(xì)胞用L-15完全培養(yǎng)基（含10% Gibco血清、150 U/mL青霉素及150 μg/mL鏈霉素）重懸，計(jì)數(shù)后按照每孔3×105個(gè)細(xì)胞的密度平均鋪在6孔板中，分為PolyI：C組和PBS對(duì)照組，每組三個(gè)重復(fù)。無(wú)菌適應(yīng)性培養(yǎng) 6 h。Poly（I：C）組加入 Poly（I：C）的 PBS 溶液，終濃度為 10 μg/mL；對(duì)照組加入等體積無(wú)菌PBS溶液。置于37℃的CO2培養(yǎng)箱中培養(yǎng)6 h后收集細(xì)胞并平均分成兩份，其中一份樣品加入1 ml Trizol吹打混勻后置于-80℃儲(chǔ)存?zhèn)溆?，此樣品用于轉(zhuǎn)錄組測(cè)序。另一份樣品按照Trizol方法[17]提取總RNA，使用NanoDrop2000和瓊脂糖凝膠電泳法檢測(cè)RNA的純度和質(zhì)量后于-80℃保存用于qPCR驗(yàn)證。

1.3.2 測(cè)序文庫(kù)的構(gòu)建和Illumina測(cè)序

總RNA質(zhì)量達(dá)到建庫(kù)要求。轉(zhuǎn)錄組文庫(kù)構(gòu)建和測(cè)序工作由北京諾禾致源科技有限公司完成。測(cè)序平臺(tái)為Hiseq，測(cè)序策略為PE150。

1.3.3 測(cè)序數(shù)據(jù)處理

原始測(cè)序序列（raw reads）按照Cellerino的方法[18]進(jìn)行過(guò)濾，進(jìn)而得到clean reads。采用Trinity法[19]對(duì)過(guò)濾后的數(shù)據(jù)進(jìn)行拼接。取每個(gè)基因中最長(zhǎng)的轉(zhuǎn)錄本作為該基因的代表[20]進(jìn)行后續(xù)的分析。

1.3.4 基因的功能注釋

基因功能注釋數(shù)據(jù)庫(kù)包括：Nr、Nt、Pfam、Swiss-Prot、KEGG/KO、COG/KOG 和 GO[21]。

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

利用RSEM[22]計(jì)算出每個(gè)unigene的表達(dá)水平，并以FPKM值[23]表示。隨后進(jìn)行差異分析,篩選閾值為 log2（fold change）|≥1，P-Value＜0.05，并繪制火山圖以直觀地展現(xiàn)差異表達(dá)基因（DEGs）。

1.3.6 免疫相關(guān)差異表達(dá)基因的篩選

結(jié)合注釋信息篩選免疫相關(guān)的差異表達(dá)基因[24]。利用HemI軟件[25]繪制免疫相關(guān)DEGs熱圖。

1.3.7 qPCR驗(yàn)證

將1.3.1所述用于qPCR驗(yàn)證的RNA樣品按照TAKARA公司反轉(zhuǎn)錄試劑盒方法合成cDNA，利用SYBRgreen II熒光染料進(jìn)行實(shí)時(shí)熒光定量PCR實(shí)驗(yàn)，反應(yīng)體系見(jiàn)表1。qPCR條件為95℃，3 min；95℃，5 sec、引物退火溫度，30 sec、39個(gè)循環(huán)；95 ℃，10 sec；65℃～95℃，每0.5℃恒溫5 s檢測(cè)溶解曲線。本試驗(yàn)所用引物見(jiàn)表2。

表1 qPCR反應(yīng)體系Table 1 Reaction system of qPCR

表2 用于免疫相關(guān)DEGs qPCR驗(yàn)證的引物Table 2 Primes used for qPCR validation of immune-related DEGs

2 結(jié)果與分析

2.1 測(cè)序數(shù)據(jù)統(tǒng)計(jì)與評(píng)估

對(duì)照及Poly I：C處理組的測(cè)序始測(cè)序序列（Raw reads）都在4 600萬(wàn)條以上。經(jīng)過(guò)濾后，各組得到的Clean reads的量占Rawreads的比例都達(dá)到了96%以上。各樣品Q20均在97%以上，Q30百分比均超過(guò)93%，CG含量在53.08%以下（表3）。

表3 測(cè)序數(shù)據(jù)輸出質(zhì)量情況Table 3 The quality of reads

2.2 轉(zhuǎn)錄組拼接及基因功能注釋

經(jīng)過(guò)混合拼接后，轉(zhuǎn)錄本（transcript）及單基因（unigene）長(zhǎng)度分布情況如圖1所示。共獲得了252 783條轉(zhuǎn)錄本和121 204個(gè)unigenes，有77 790條轉(zhuǎn)錄本超過(guò)2 000 bp，30 170條unigene長(zhǎng)度超過(guò)1 000 bp，轉(zhuǎn)錄本平均長(zhǎng)度為1 959 bp，N50=4 131，單基因平均長(zhǎng)度為 1 121 bp，N50=2 000（表 4）。

表4 拼接長(zhǎng)度分布情況Table 4 The distribution of splicing length bp

圖1 轉(zhuǎn)錄本及unigene拼接長(zhǎng)度分布情況Figure 1 Distribution of transcript and unigene splice length

Unigenes在7個(gè)數(shù)據(jù)庫(kù)中成功注釋，45 352個(gè)unigene在Nr中成功注釋，在Nt數(shù)據(jù)庫(kù)中成功注釋了117 078個(gè)unigene，有8 296個(gè)unigene成功注釋在swiss-prot數(shù)據(jù)庫(kù)中；有23 171和7 710個(gè)unigene分別成功注釋在GO和KEGG數(shù)據(jù)庫(kù)中（圖2）。

圖2 Unigene在各個(gè)數(shù)去庫(kù)中的注釋情況Figure 2 Unigene annotations in database

共計(jì)2 488個(gè)unigene在所有數(shù)據(jù)庫(kù)中均成功被注釋，117 722個(gè)unigene至少在一種數(shù)據(jù)中被成功注釋，但有3 482個(gè)unigene并未成功注釋（表5）?；蚬δ茏⑨屖菍?duì)比多個(gè)物種的基因信息完成的（圖 3）。

圖3 注釋所用參考物種分類情況Figure 3 The speciation and proportion of annotation

表5 注釋情況統(tǒng)計(jì)Table 5 The statistics of success rate of annotation

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

差異表達(dá)基因分析共挖掘出402個(gè)差異表達(dá)基因（252個(gè)顯著上調(diào)和150個(gè)顯著下調(diào)），其余基因未出現(xiàn)顯著變化（圖 4）。

圖4 差異表達(dá)基因的火山圖Figure 4 Volcano plot of differential expressed genes

2.4 差異基因GO及KEEG富集分析

GO富集分析將所有差異表達(dá)基因富集到生物學(xué)過(guò)程（biological process，BP），細(xì)胞組分（cellular component，CC）和分子功能（molecular function，MF）三大類別中，在生物學(xué)過(guò)程中，免疫反應(yīng)（immune response）和免疫過(guò)程（immune system progress）等出現(xiàn)較高比例的富集（圖5）。KEGG富集分析將差異表達(dá)基因富集到各種信號(hào)轉(zhuǎn)導(dǎo)通路中，其中Tolllike receptors信號(hào)通路和RIG-1-like receptors信號(hào)通路等出現(xiàn)大量富集（圖6）。

圖5 差異表達(dá)基因GO富集分析Figure 5 Gene Ontology enrichment analysis of DEGs

圖6 DEGs的KEGG富集分析Figure 6 KEGG enrichment analysis for DEGs

2.5 免疫相關(guān)差異基因分析

基于功能注釋和富集結(jié)果，共篩選到28個(gè)免疫相關(guān)基因，利用HEMI軟件將兩組的免疫相關(guān)基因繪制熱圖，以直觀顯示二者變化差異（圖7）。

圖7 免疫相關(guān)DEGs熱圖Figure 7 Heat map of immune-related DEGs

2.6 qPCR驗(yàn)證

隨機(jī)選取9個(gè)免疫相關(guān)基因進(jìn)行實(shí)時(shí)熒光定量分析，以驗(yàn)證轉(zhuǎn)錄組測(cè)序數(shù)據(jù)的準(zhǔn)確性。結(jié)果表明，Poly I：C刺激后所選的9個(gè)基因均出現(xiàn)不同程度的應(yīng)答，與轉(zhuǎn)錄組數(shù)據(jù)相符（圖8）。

圖8 免疫相關(guān)DEGs在轉(zhuǎn)錄組數(shù)據(jù)和qPCR數(shù)據(jù)中的表達(dá)水平Figure 8 Comparison of the fold change expression of DEGs by RNA-seq and qPCR

3 討論與結(jié)論

豚鹿是中國(guó)國(guó)家Ⅰ級(jí)重點(diǎn)保護(hù)的野生動(dòng)物，其基因組信息尚未公布[26]。本研究對(duì)Poly I：C刺激后的豚鹿外周血淋巴細(xì)胞進(jìn)行了無(wú)參轉(zhuǎn)錄組測(cè)序。轉(zhuǎn)錄組測(cè)序數(shù)據(jù)質(zhì)量主要從堿基測(cè)序錯(cuò)誤率（Q20和Q30分別表示測(cè)序錯(cuò)誤率為1%和0.1%的堿基比例）和堿基分布情況進(jìn)行評(píng)估[27]。本次測(cè)序共獲得318 259 890 bp的原始數(shù)據(jù)，過(guò)濾后得到的clean reads占原始數(shù)據(jù)的97.80%，Q20和Q30均占較高比例，且CG含量較低。組裝后的轉(zhuǎn)錄本及unigene長(zhǎng)度均大于300 bp，說(shuō)明測(cè)序建庫(kù)及拼接工作質(zhì)量良好，具有較高的可信度。在基因功能注釋方面，共有117 722條unigene被注釋到至少一個(gè)數(shù)據(jù)庫(kù)中，僅有2.88%的unigene未得到注釋?？傮w注釋情況較好，能夠有效完成基因功能注釋。本研究豐富了豚鹿的基因資源，可為豚鹿的抗病毒免疫研究提供高質(zhì)量的基因數(shù)據(jù)參考。

有研究表明，PolyI：C能夠模擬病毒的雙鏈RNA，進(jìn)而引起機(jī)體免疫反應(yīng)[28]。差異表達(dá)基因數(shù)目統(tǒng)計(jì)結(jié)果表明，Poly I：C刺激后可以引起豚鹿外周血淋巴細(xì)胞基因的差異表達(dá)。差異表達(dá)基因經(jīng)GO富集后被分為生物學(xué)過(guò)程，細(xì)胞組分和分子功能三大類中。生物學(xué)過(guò)程中，DEGs富集比例較高的亞類且與免疫密切相關(guān)的“免疫反應(yīng)”和“免疫系統(tǒng)過(guò)程”等。細(xì)胞組分中，“胞外區(qū)域”占據(jù)較高的比例，可能與胞外病毒識(shí)別有關(guān)。在分子功能中，“細(xì)胞因子激活”“受體結(jié)合”等免疫相關(guān)亞類具有較高的富集比例。豚鹿外周血淋巴細(xì)胞在Poly I：C刺激后在抗病毒相關(guān)免疫調(diào)控信號(hào)轉(zhuǎn)導(dǎo)途徑中出現(xiàn)了復(fù)雜的信號(hào)轉(zhuǎn)導(dǎo)反應(yīng)。尤其是在機(jī)體抗病毒免疫過(guò)程中起著重要作用的RLR和TLR信號(hào)通路[29]。由此可見(jiàn)，Poly I：C能夠有效模擬病毒刺激機(jī)體產(chǎn)生免疫反應(yīng)，為后續(xù)抗病毒免疫相關(guān)基因的篩選奠定基礎(chǔ)。進(jìn)一步篩選免疫相關(guān)的差異基因發(fā)現(xiàn)，Toll樣受體（TLRs）信號(hào)通路中的IRF7、IL6等基因在刺激后有較高的表達(dá)；另外，補(bǔ)體系統(tǒng)及趨化因子信號(hào)通路相關(guān)基因如PTX3，C1RL，CXCL10等也出現(xiàn)上調(diào)。

干擾素調(diào)節(jié)因子7（IRF7）在抗病毒免疫反應(yīng)中具有重要的作用，能夠參與I型干擾素的產(chǎn)生并調(diào)節(jié)先天免疫反應(yīng)[30]。PolyI：C是合成的病毒雙鏈RNA類似物，它可以用來(lái)模擬病毒感染并引發(fā)抗病毒免疫反應(yīng)過(guò)程產(chǎn)生I型干擾素和炎性細(xì)胞因子來(lái)發(fā)揮抗病毒免疫的作用[31]。在本試驗(yàn)中，polyI：C刺激后上調(diào)了IRF7的表達(dá)，并誘導(dǎo)了IL-6的表達(dá)。此外，長(zhǎng)鏈正五聚蛋白（PTX3）能夠識(shí)別多種多樣的病原體，通過(guò)結(jié)合C1q調(diào)節(jié)補(bǔ)體活性，并促進(jìn)巨噬細(xì)胞和樹(shù)突狀細(xì)胞對(duì)病原體的識(shí)別[32]。PTX3能夠結(jié)合并抑制人和鼠巨噬細(xì)胞病毒[33]和多種流感病毒[34]。本研究中polyI：C刺激外周淋巴細(xì)胞后上調(diào)了PTX3和補(bǔ)體受體C1RL基因的表達(dá)，推測(cè)補(bǔ)體系統(tǒng)在豚鹿抗病毒免疫反應(yīng)中起到了重要的作用。

在涉及轉(zhuǎn)錄組技術(shù)的研究中，使用qPCR對(duì)轉(zhuǎn)錄組數(shù)據(jù)的準(zhǔn)確性進(jìn)行驗(yàn)證是常用的方法[35]。因此，本研究從篩選出的免疫相關(guān)基因中隨機(jī)挑選了9個(gè)基因進(jìn)行了qPCR驗(yàn)證，結(jié)果表明，二者數(shù)據(jù)呈現(xiàn)出一致的變化規(guī)律。說(shuō)明轉(zhuǎn)錄組數(shù)據(jù)具有較高的可信度。

本研究利用Poly I：C對(duì)豚鹿外周血淋巴細(xì)胞進(jìn)行刺激，獲得二者的基因表達(dá)譜；并深度挖掘免疫相關(guān)基因，以及其表達(dá)差異。揭示豚鹿在病毒防御中的免疫應(yīng)答機(jī)制，有利于人工繁育豚鹿過(guò)程中定制合理的飼養(yǎng)方案和疾病預(yù)防措施，并且為疫苗佐劑的開(kāi)發(fā)及應(yīng)用提供了理論基礎(chǔ)。