家蠶感染二分濃核病毒（鎮(zhèn)江株）的數字基因表達譜分析

高 坤，尚夢珂，錢荷英，覃光星，郭錫杰

（江蘇科技大學生物技術學院/中國農業(yè)科學院蠶業(yè)研究所，江蘇鎮(zhèn)江 212018）

家蠶感染二分濃核病毒（鎮(zhèn)江株）的數字基因表達譜分析

高 坤，尚夢珂，錢荷英，覃光星，郭錫杰

（江蘇科技大學生物技術學院/中國農業(yè)科學院蠶業(yè)研究所，江蘇鎮(zhèn)江 212018）

【目的】篩選家蠶與二分濃核病毒（Bombyx mori bidensovirus Zhenjiang strain，BmBDV-ZJ）感染相關的差異表達基因，鑒定家蠶與病毒感染有關的調控基因，為進一步闡明家蠶抗BmBDV-ZJ的分子機制提供理論依據。【方法】采用Illumina高通量測序技術，構建家蠶品種JS口服感染BmBDV-ZJ的數字基因表達譜，為排除個體間差異，以10頭蠶作為一個樣本用于DGE檢測。樣本中基因的差異表達檢測通過嚴格的運算法則進行，對差異檢驗的P值（P value）作多重假設檢驗校正，通過控制FDR（false discovery rate）來決定P值的域值。本研究中，差異表達基因定義為FDR≤0.001且差異倍數在2倍及以上（|log2ratio|≥1）的基因。采用基因本體論（GO）分類體系確定所有差異表達基因可能的功能。用GO計算P值和bonferoni校正。選用校正P值≤0.05作為基因組顯著富集的閾值。WEGO軟件用來視化、比較和繪制GO注釋結果。利用KEGG數據庫進行通路富集分析,進一步確定顯著富集代謝途徑或信號傳導途徑，Q值≤0.05的通路指定為DGEs中的顯著富集通路。通過qRT-PCR方法對部分差異表達基因進行驗證。【結果】感染組和對照組分別得到4 850 663和4 875 307個原始標簽,去除低質量標簽后，分別得到4 757 934和4 788 406個清潔標簽，對應的標簽種類數量分別為62 436和63 680種。兩個文庫間的清潔標簽和清潔標簽種類的數量在不同拷貝區(qū)間分布類似，感染組和對照組樣本的測序量分別為3.5 M和3.7 M，測序深度符合試驗的要求，兩樣本的DGE數據是可信的。將這兩個DGE數據庫的所有清潔標簽與家蠶參考基因庫進行比對，在對照組與感染組中，分別有36.39%和45.30%的清潔標簽可以比對到基因。另有50.02%和43.34%的清潔標簽可以比對到家蠶參考基因組，剩余的未知標簽分別占清潔標簽總數的13.59%和12.35%。共發(fā)現了447個差異表達基因，其中306個上調表達，141個下調表達。分別有218、147和179個差異表達基因涉及GO 3個本體中的分子功能、細胞組分和生物過程。利用KEGG公共數據庫進行Pathway顯著性富集分析，注釋到的基因總數為8 473個。447個差異表達基因經鑒定后，其中的330個基因被歸類到151個KEGG路徑中。差異表達基因顯著富集的Pathway（Q值≤0.05）有19個，其中最顯著富集的是細胞質中DNA識別通路。挑選了24個差異表達基因進行qRT-PCR驗證，其中20個基因的差異表達趨勢與DGE的結果一致。其中在DNA識別通路中共檢測到9個差異表達基因，BGIBMGA009408-TA、BGIBMGA004913-TA、BGIBMGA011753-TA均為編碼RNA聚合酶III的基因，表達量均上調，是對照組的4.3、2.3、3.4倍。【結論】構建了3齡家蠶JS感染BmBDV-ZJ后28 h感染組及對照組幼蟲的數字基因表達譜，Pathway顯著性富集分析和qRT-PCR驗證顯示，家蠶感染BmBDV-ZJ后可能通過啟動胞質內DNA識別通路來感應入侵病毒的異源DNA成分并迅速啟動天然免疫抵御BmBDV-ZJ病毒感染,為研究BmBDV-ZJ侵染家蠶和家蠶抵御BmBDV感染的分子機制打下了基礎。

家蠶；數字基因表達譜；二分濃核病毒；qRT-PCR

0 引言

【研究意義】家蠶二分濃核病毒鎮(zhèn)江株（Bombyx mori bidensovirus Zhenjiang strain，BmBDV-ZJ）是感染家蠶的一種重要病毒，主要通過食下感染，引起家蠶的濃核病。患病蠶呈現空頭，下痢，吐液等癥狀，是嚴重影響蠶桑生產的病毒病之一，每年給養(yǎng)蠶業(yè)造成較大的經濟損失。目前，對于該病毒感染家蠶的致病機理和家蠶抵御該病毒感染的應答機制均不清楚。篩選鑒定家蠶與該病毒感染相關的差異表達基因，進一步研究病毒對家蠶致病的分子機理，對家蠶抗病毒育種和家蠶濃核病的有效防治具有重要意義。【前人研究進展】中國早在1959年就證實了生產上所發(fā)生的家蠶空頭性軟化病是由病毒引起，隨后各國科學家分離得到了該病毒的不同株系，分別命名為伊那株（DNV-1）[1]、佐久株（DNV-2）[2-3]、中國（鎮(zhèn)江）株（DNV-3）[4]。之前認為該病毒均屬于細小病毒科（Parvoriridae）濃核病毒屬（Densovirus，DNV）的蠶濃核病毒（B. mori densonucleosis virus）；后來發(fā)現DNV-2和DNV-3基因組包含兩個節(jié)段的單鏈DNA且不采用滾環(huán)復制（細小病毒科常用的復制方式）[5]，因此將DNV-2和DNV-3重新劃分為一個新設定的二分DNA病毒科二分濃核病毒屬（Bidensovirus，BDV）[6]。本研究所用的病毒株是1981年由中國農業(yè)科學院蠶業(yè)研究所（鎮(zhèn)江）錢元駿等分離得到的家蠶二分濃核病毒（鎮(zhèn)江）株（BmBDV-ZJ），與伊那株（DNV-1）在蠶品種感受性、病毒的物理化學性狀、血清學特性等方面明顯不同[7]。BDV-ZJ主要在家蠶幼蟲中腸柱狀上皮細胞的細胞核中復制，引起幼蟲的濃核病，不形成多角體[8]。其病毒基因組包括兩段不同的DNA單鏈（6 543和6 022 bp），一個病毒粒子只包含其中一條單鏈DNA，因此是兩種不同病毒粒子的混合物[9]。成熟的病毒粒子都會隨感染細胞的破裂而釋放，進而感染鄰近細胞。大部分家蠶品種對BmBDV-ZJ是易感的，錢元駿等[10]通過對380多個家蠶品種的抗BmBDV-ZJ性能比較研究和不同抗性品種進行雜交試驗證實家蠶對BmBDV-ZJ抗性為隱性遺傳，同時也受微效多基因的影響。有的品種即使接種高濃度的病毒也完全不發(fā)病。日本科學家[11]研究發(fā)現了一個與DNV-2易感性有關的家蠶基因nsd-2，抗性品種中該基因開放閱讀框中大約6 000個堿基缺失，導致其編碼的具有12個跨膜域的氨基酸轉運膜蛋白的缺失，且該蛋白僅在家蠶的中腸中表達；而通過轉基因技術修復該基因的缺失后可以使抗性品種對BmBDV-2病毒易感，該研究結果表明病毒識別的位點可能位于nsd-2基因缺失的膜蛋白部分，然而病毒是如何與這種膜蛋白相互作用仍不清楚。裘智勇等[12-13]研究表明，家蠶對BmBDV-ZJ的抵抗性可能與家蠶中腸細胞表面的某些特殊蛋白因子有關。BAO等[14]通過抑制消減雜交技術，研究了兩個不同抗性家蠶品種JS（對BmBDV-ZJ感受性）和NIL（對BmBDV-ZJ非感受性）感染BmBDV-ZJ后基因的相對表達變化，發(fā)現了在抗性品種中有11個基因明顯上調，推測這些基因與NIL的抗BmBDV-ZJ感染復制有關。【本研究切入點】病毒生活史的每一步都受病毒和宿主間分子相互作用的調控，通過家蠶與BmBDV-ZJ病毒感染相關的基因差異表達分析，鑒定可能與病毒復制有關的宿主細胞內調控因子，有助于闡明家蠶抗BmBDV-ZJ的抗病毒機制。【擬解決的關鍵問題】通過對易感家蠶品種JS的差異基因表達譜和差異基因可能參與的信號通路分析，篩選和尋找更多可能與病毒感染相關的差異表達基因。

1 材料與方法

試驗于2013年9月至2015年9月在中國農業(yè)科學院蠶業(yè)研究所病理研究室完成。

1.1 供試材料與試劑

1.1.1 家蠶及病毒株 家蠶品種JS由中國農業(yè)科學院國家蠶種質資源保存中心提供；家蠶二分濃核病毒（BmBDV-ZJ）由中國農業(yè)科學院蠶業(yè)研究所家蠶病理研究室繁殖保存；BmBDV-ZJ抗血清由中國農業(yè)科學院蠶業(yè)研究所家蠶病理研究室制備保存。

1.1.2 主要試劑 PrimeScriptTMRT Reagent Kit、

SYBR Premix Ex TaqTMKit購自TaKaRa公司；總RNA提取試劑盒、Trizol、M-MLV反轉錄酶購自Invitrogen公司；DNA marker DL 2000購自TaKaRa公司；其他均為國產或進口分析純試劑。

1.2 試驗方法

1.2.1 家蠶幼蟲添毒感染 取二分濃核病毒感染蠶的中腸干粉0.4 g，加入8 mL 2% Na2CO3溶液研磨勻漿，靜止4—5 min，加5 mL蒸餾水至50 mL，3 000 r/min離心20 min，取上清加3倍體積0.2 mol·L-1醋酸，處理30 min，用Na2CO3調pH值至6—7，用蒸餾水稀釋至1/1 000備用，過濾除菌備用。

試驗用家蠶幼蟲JS在標準溫度和光照條件下飼養(yǎng)至3齡起蠶。計數40頭蠶為感染組（添食BmBDV-ZJ），另40頭為對照組（添食滅菌水）。將制備病毒懸浮液均勻涂到桑葉上喂食感染組家蠶，對照組喂食同樣體積的滅菌水涂過的桑葉，12 h后全部喂飼正常桑葉。考慮到個體間差異可能帶來的影響，每10條蠶收集一管提取RNA，共設置3組生物學重復。最后剩余10頭蠶用來檢測病毒感染發(fā)病情況。

1.2.2 病毒感染的確認 在一塊7 cm×10 cm的瓊脂板（1%瓊脂糖）上打7個孔，中間孔加入抗血清，周圍6個孔加待測樣品，將瓊脂板放入補濕飯盒中，在20—30℃條件下孵育1 d觀察有無沉淀帶出現。

1.2.3 數字基因表達譜（digital gene expression, DGE）測序及分析 由華大基因公司提供技術支持，對28 h時間點的感染組及對照組JS家蠶幼蟲各取10 μg總RNA，經過純化、反轉錄、酶切和PCR線性擴增后，使用Illumina HiSeqTM2000進行測序。差異表達基因的篩選及其Gene Ontology（GO）、Pathway顯著性富集分析參照文獻[15-16]進行。

1.2.4 qRT-PCR分析 通過實時熒光定量PCR（qRTPCR）對DGE篩選得到的部分差異表達基因進行驗證。感染組和對照組各3管樣品，以每管樣品提取的總RNA作為模板，用PrimeScriptTMRT Reagent Kit（TaKaRa）進行反轉錄合成第一鏈cDNA。以家蠶Actin 3作為內參基因，使用SYBR? Premix Ex TaqTM（TaKaRa）試劑盒在ABI PRISM? 7300 檢測系統(tǒng)上進行實時定量PCR，每個模板做3次重復。反應條件：94℃ 10 min變性，94℃ 15 s，60℃ 31 s，40個循環(huán)。數據分析參考文獻[17]進行，引物序列見表1。

2 結果

2.1 病毒感染的確定

通過雙向免疫擴散法診斷病毒的感染。在BmBDV-ZJ感染28 h后，所有接種感染的家蠶幼蟲其中腸研磨液與血清孔之間均出現了白色沉淀帶，而與健康蠶（對照組）中腸研磨液及BmCPV、BmNPV感染蠶研磨液均不起反應，表明感染組家蠶幼蟲全部被BmBDV-ZJ感染。

2.2 DGE數據庫分析

采用Illumina高通量測序技術，感染組和對照組分別得到4 850 663和4 875 307個原始標簽，去除低質量標簽后，分別得到4 757 934和4 788 406個清潔標簽，對應的標簽種類數量分別為62 436和63 680種（表2）。兩個文庫間的清潔標簽和清潔標簽種類的數量在不同拷貝區(qū)間分布類似（圖1），感染組和對照組樣本的測序量分別為3.5 M和3.7 M（圖2），測序深度符合試驗的要求，兩樣本的DGE數據是可信的。

圖1 對照組和感染組的清潔標簽和清潔標簽種類在不同拷貝區(qū)的分布Fig. 1 Distribution of total clean tags and distinct clean tags in each library

表1 定量PCR引物Table 1 Primers for qRT-PCR

2.3 標簽比對及標準化處理

將所有清潔標簽與家蠶參考基因庫進行比對，在對照組與感染組中，分別有36.39%和45.30%的清潔標簽可以比對到基因，占總清潔標簽種類數的30.56%和37.41%。另有50.02%和43.34%的清潔標簽可以比對到家蠶參考基因組，剩余的未知標簽分別占清潔標簽總數的13.59%和12.35%（表2）。

表2 兩個樣本中標簽分布表Table 2 Distributions of tags between the two libraries

圖2 數字基因表達譜的測序飽和趨勢圖Fig. 2 Trends of saturation of DGEs

2.4 差異基因分析及其驗證

篩選了發(fā)現錯誤率（false discovery rate）FDR≤0.001且差異倍數在2倍及以上（|log2ratio|≥1）的基因作為差異表達基因，共發(fā)現差異表達基因447個。在感染組中，上調表達基因306個，下調表達基因141個。為了驗證DGE數據的準確性，挑選了24個差異表達基因進行RT-PCR驗證，其中20個基因表達情況與DGE的結果一致（圖3），說明DGE數據基本可以準確反映樣本的基因表達情況。

2.5 GO和KEGG分析

GO共有3個本體，分別描述基因的分子功能、所處的細胞組分、參與的生物過程。GO結果顯示分別有147、179、218個差異基因分配在這3個本體中，其中在核糖核蛋白復合體、核糖體、核糖體亞基、大分子復合物、細胞內組分、蛋白酶體調控組分、代謝過程、基因表達和結構分子活性中基因富集程度較高，校正后的P值≤0.05（表3）。

為了進一步了解基因的生物學功能，利用KEGG公共數據庫進行Pathway顯著性富集分析，注釋到的基因總數為8 473個，447個差異表達基因經鑒定后，其中的330個基因被歸類到151個KEGG路徑中。定義為在差異表達基因中顯著富集的Pathway（Q值≤0.05）有19個（表4），其中最顯著富集的是細胞質中DNA識別通路，其次是RNA聚合酶、核糖體、嘧啶代謝、蛋白酶體、抗原加工與呈遞、代謝途徑等。

圖3 候選基因在對照和感染BmBDV-ZJ的家蠶幼蟲中的差異表達Fig. 3 Differential expression levels of candidate genes in BmBDV-ZJ-infected and control B. mori larvae

3 討論

為了研究BmBDV-ZJ感染家蠶的可能分子機制和宿主的應答反應，基于本研究構建JS家蠶感染組和對照組的兩個DGE文庫分析，篩選得到447個差異表達基因，其中感染組有306個上調表達基因和141個下調表達基因，該差異基因的上調和下調對比結果與之前家蠶感染質型多角體病毒的DGE結果類似，都是上調表達的基因數目明顯高于下調表達的基因數目[15-16,18]。相關研究表明，家蠶在感染大部分病原微生物時的差異表達基因都有上調基因明顯多于下調基因的趨勢，LIU等[19]通過基因芯片研究發(fā)現家蠶感染革蘭氏陽性菌（Serratia marcescens）后有172個基因上調表達，61個基因下調表達；感染革蘭氏陰性菌（Staphylococcus aureus）后125個基因上調表達，104個基因下調表達；感染真菌（Beauveria bassiana）后有133個基因上調表達，24個基因下調表達；而被PBS處理后有201個基因上調表達和40個基因下調表達。可見無論是病毒感染、細菌感染還是真菌感染都會引起家蠶更多基因的上調表達來應對病原的入侵。

表3 基因本體分類Table 3 Gene Ontology classification

表4 顯著富集的PathwayTable 4 Significantly enriched pathways

利用KEGG公共數據庫進行Pathway顯著性富集分析，其中最顯著富集的是細胞質中DNA識別通路（cytosolic DNA-sensing pathway），該通路可以有效識別入侵病毒的異源DNA成分并迅速啟動天然免疫，以及隨后的特異性免疫應答來對病毒進行清除，是機體抵抗病毒感染的重要機制[20-21]。DNA感受器（DNA sensor）是宿主感受外源入侵DNA和免疫防御的橋梁，可以特異地識別病毒等外源DNA進而激活下游的免疫信號途徑，并通過誘導表達抗病毒蛋白來抑制病毒的復制并向周圍的細胞示警[22-23]。目前已經有超過10種DNA感受器被發(fā)現。例如，Toll樣受體9（TLR9）家族[24]、視黃酸誘導基因蛋白RIG-I相關受體（RLRs）家族[25]、DNA依賴的RNA聚合酶III[26]、干擾素誘導蛋白16（IFI16）[27]；DExD-H框解旋酶超家族[28-30]等。RNA聚合酶III可以識別侵入細胞內的DNA病毒，然后以病毒DNA為模板，將病毒信息轉錄合成一種5′三磷酸化的特殊的雙鏈RNA，進而被RIG-I分子所識別。RIG-I可以激活I型干擾素產生，也可以激活NF-κB等信號通路并誘導抗病毒基因的表達，從而抑制病毒等病原體復制[31-32]。家蠶感染BmBDV-ZJ后在DNA識別通路中共檢測到9個差異表達基因，其中BGIBMGA009408-TA、BGIBMGA004913-TA、BGIBMGA011753-TA均為編碼RNA聚合酶III的基因，表達量均上調，是對照組的4.3、2.3、3.4倍（圖2）。其功能可能作為家蠶細胞質中的DNA感受器，通過識別病毒DNA激活家蠶抗BmBDV-ZJ感染的免疫應答反應。

干擾素-γ誘導的溶酶體巰基還原酶（IFN-γ inducible lysosomal thiol reductase，GILT）在脊椎動物適應性免疫中的MHCII類抗原加工和呈遞過程中起著關鍵作用，可以催化未折疊的天然抗原蛋白二硫鍵的斷裂，進而對其進行酶解加工[33]。雖然適應性免疫和干擾素的產生已明確只在脊椎動物中存在，但是GILT基因家族在脊椎動物和無脊椎動物中普遍存在，如在蝦[34]、蟹[35]、鮑[36]、果蠅[37]和線蟲[38]中均發(fā)現了該基因。無脊椎動物中的GILT基因在感染細菌和病毒后也都出現上調表達趨勢，推測其功能可能不同于脊椎動物中的抗原加工和呈遞，而是參與了無脊椎動物中的某種先天性免疫信號途徑。家蠶感染BmBDV-ZJ病毒后基因BGIBMGA003741-TA和BGIBMGA008360-TA均上調表達，且這2個基因都含有GILT保守結構域，推測為家蠶中的GILT基因。其功能是否類似于接頭蛋白干擾素刺激基因，其上調表達是否可以進一步誘導相關抗病毒蛋白的產生還需要進一步試驗驗證。

BGIBMGA013054-TA作為缺血再灌注誘導蛋白的同源基因，其主要功能是調節(jié)細胞內物質裝配和相關轉運蛋白的活性[39]。病毒感染家蠶后，會利用宿主細胞的蛋白合成系統(tǒng)來進行病毒蛋白的復制，這些異源蛋白的產生，加速了細胞的物質運輸過程，以增加宿主細胞對病毒蛋白的多重耐受性。家蠶感染BmBDV-ZJ后該基因的上調表達可能與細胞內一些蛋白的定位有關。

BmBDV-ZJ通常引起的是慢性病，感染的中腸上皮細胞不像感染BmBDV-1的中腸上皮細胞那樣容易脫落，而是通過增加細胞數目，引起上皮組織折疊，在10—20 d間死亡，少數幼蟲可以化蛹。因此本研究分析發(fā)現家蠶感染BmBDV-ZJ后，很多與細胞增殖相關的轉錄因子和蛋白都發(fā)生了明顯上調，如叉頭框轉錄因子G1（BGIBMGA010297-TA，forkfead transcription factor G1）、富含亮氨酸重復序列蛋白47（BGIBMGA013072-TA，leucine-rich repeat-containing protein 47-like protein）、賴氨酸去甲基化和組氨酸脫氫的雙功能酶（BGIBMGA014563-TA，bifunctional lysine-specific demethylase and histidyl-hydroxylase NO66）等基因。初步推測家蠶可以通過對這些基因的上調表達促進細胞增殖來增加細胞的數目，進而替代因感染BmBDV-ZJ而失去正常生理功能或凋亡的中腸柱狀細胞。

BGIBMGA007315-TA是家蠶的泛素特異蛋白酶7（Ubiquitin specific protease 7，USP7），定量結果顯示為下調，與DGE結果相反。該基因編碼的一種泛素化酶，可以水解Mdm2（p53的E3泛素連接酶）對p53蛋白進行去泛素化，保護p53不被S26蛋白酶體降解，進而調節(jié)p53與Mdm2的穩(wěn)定性[40]。相關研究表明，EB病毒的核抗原1蛋白可以通過與USP7結合，破壞p53的穩(wěn)定性，有利于病毒的潛伏感染[41]。BmBDV-ZJ感染家蠶后病毒是否會與家蠶的USP7相互作用，USP7下調表達是否也會影響到家蠶p53蛋白的穩(wěn)定性還需要進一步驗證。

綜上所述，差異表達基因最顯著富集的是細胞質DNA識別通路，該通路涉及的基因BGIBMGA009408-TA、BGIBMGA004913-TA、BGIBMGA011753-TA、BGIBMGA003741-TA、BGIBMGA008360-TA和BGIBMGA013054-TA將進一步進行功能驗證，明確DNA識別通路在家蠶抗BmBDV-ZJ病毒感染中的先天免疫機制。其他大部分差異基因功能未知，因此后期主要工作是對這些未知功能的差異表達基因做進一步鑒定及其參與的抗病毒機制研究。

4 結論

構建了3齡家蠶JS感染BmBDV-ZJ后28 h感染組及對照組幼蟲的數字基因表達譜，Pathway顯著性富集分析和qRT-PCR驗證顯示，家蠶感染BmBDV-ZJ后可能通過啟動胞質內DNA識別通路來感應入侵病毒的異源DNA成分并迅速啟動天然免疫抵御BmBDV-ZJ病毒感染。

[1] BANDO H, CHOI H, ITO Y, KAWASE S. Terminal structure of a densovirus implies a hairpin transfer replication which is similar to the model for AAV. Virology, 1990, 179(1): 57-63.

[2] BANDO H, CHOI H, ITO Y, NAKAGAKI M, KAWASE S. Structural analysis on the single-stranded genomic DNAs of the virus newly isolated from silkworm: the DNA molecules share a common terminal sequence. Archives of Virology, 1992, 124(1/2): 187-193.

[3] BANDO H, HAYAKAWA T, ASANO S, SAHARA K, NAKAGAKI M, IIZUKA T. Analysis of the genetic information of a DNA segment of a new virus from silkworm. Archives of Virology, 1995, 140(6): 1147-1155.

[4] IWASHITA Y, CHUN C Y. The development of a densonucleosis virus isolated from silkworm larvae, Bombyx mori, of China//AKAI H, KING R C, MOROHOSHI S. The Ultrastructure and Functioning of Insect Cell. Tokyo: Society for Insect Cells Japan, 1982: 161-164.

[5] HAYAKAWA T, ASANO S, SAHARA K, IIZUKA T, BANDO H. Detection of replicative intermediate with closed terminus of Bombyx densonucleosis virus. Archives of Virology,1997, 142: 393-399.

[6] ADAMS M J, CARSTENS E B. Ratification vote on taxonomic proposals to the International Committee on Taxonomy of Viruses. Archives of Virology, 2012, 157: 1411-1422.

[7] 錢元駿, 郭錫杰, 胡雪芳, 黃可威, 渡部仁. 我國和日本家蠶DNV的血清學關系. 蠶業(yè)科學, 1985, 11(4): 241-242. QIAN Y J, GUO X J, HU X F, HUANG K W, WATANABE H. The serological relationship between China isolate densonucleosis virus and Japan isolate densonucleosis virus. Acta Sericologica Sinica, 1985, 11(4): 241-242. (in Chinese)

[8] 郭錫杰, 錢元駿, 胡雪芳, 王紅林. 我國家蠶濃核病毒(DNV)寄生部位研究. 蠶業(yè)科學, 1985, 11(2): 93-98. GUO X J, QIAN Y J, HU X F, WANG H L. Studies on locations of Bombyx mori densonucleosis virus (China isolate) invasion. Acta Sericologica Sinica, 1985, 11(2): 93-98. (in Chinese)

[9] WANG Y J, YAO Q, CHEN K P, WANG Y, LU J, HAN X. Characterization of the genome structure of Bombyx mori densovirus (China isolate). Virus Genes, 2007, 35: 103-108.

[10] 錢元駿, 胡雪芳, 孫玉昆, 戴仁鳴. 家蠶濃核病毒的研究. 蠶業(yè)科學, 1986, 12(2): 89-94. QIAN Y J, HU X F, SUN Y K, DAI R M. Studies on Bombyx mori densonucleosis virus. Acta Sericologica Sinica, 1986, 12(2): 89-94. (in Chinese)

[11] ITO K, KIDOKORO K, SEZUTSU H, NOHATA J, YAMAMOTO K, UCHINO K, KALYEBI A, EGUCHI R, HARA W, TAMURA T, KATSUMA S, MITA K, KADONO-OKUDA K. Deletion of a gene encoding an amino acid transporter in the midgut membrane causes resistance to a Bombyx parvo-like virus. Proceedings of the National Academy of Sciences of the United Dtates of America, 2008, 105(21): 7523-7527.

[12] 裘智勇, 李木旺, 沈興家, 郭錫杰. 家蠶對濃核病毒(鎮(zhèn)江株)抵抗性和感受性品種的中腸組織蛋白比較分析. 蠶業(yè)科學, 2008, 34(2): 244-249. QIU Z Y, LI M W, SHEN X J, GUO X J. Comparative analysis of proteins extracted from midgut of silkworm strains susceptible and non-susceptible to Bomby mori densovirus (Zhenjiang strain). Acta Sericologica Sinica, 2008, 34(2): 244-249. (in Chinese)

[13] 裘智勇, 李木旺, 覃光星, 劉挺, 沈興家, 郭錫杰. 家蠶對濃核病毒中國鎮(zhèn)江株抵抗性機制的初步研究. 蠶業(yè)科學, 2007, 33(4): 596-601. QIU Z Y, LI M W, QIN G X, LIU T, SHEN X J, GUO X J. Primary studies on mechanism of silkworm (Bombyx mori) resistance to densovirus China (Zhenjiang) strain. Acta Sericologica Sinica, 2007, 33(4): 596-601. (in Chinese)

[14] BAO Y Y, LI M W, ZHAO Y P, GE J Q, WANG C S, HUANG Y P, ZHANG C X. Differentially expressed genes in resistant and susceptible Bombyx mori strains infected with a densonucleosis virus. Insect Biochemistry and Molecular Biology, 2008, 38(9): 853-861.

[15] GAO K, DENG X Y, QIAN H Y, QIN G X, HOU C X, GUO X J. Cytoplasmic polyhedrosis virus-induced differential gene expression in two silkworm strains of different susceptibility. Gene, 2014, 539: 230-237.

[16] GAO K, DENG X Y, QIAN H Y, QIN G X, GUO X J. Digital gene expression analysis in the midgut of 4008 silkworm strain infected with cytoplasmic polyhedrosis virus. Journal of Invertebrate Pathology, 2014, 115(1): 8-13.

[17] 高坤, 鄧祥元, 裘智勇, 覃光星, 郭錫杰. 家蠶感染質型多角體病毒 (BmCPV)后中腸組織差異蛋白質分析. 中國農業(yè)科學, 2013, 46(13): 2796-2807. GAO K, DENG X Y, QIU Z Y, QIN G X, GUO X J. Comparative analysis of differential proteins from midgut of silkworm induced by cytoplasmic polyhedrosis virus infection. Scientia Agricultura Sinica, 2013, 46(13): 2796-2807. (in Chinese)

[18] GUO R, WANG S M, XUE R Y, CAO G L, HU X L, HUANG M L, ZHANG Y Q, LU Y H, ZHU L Y, CHEN F, LIANG Z, KUANG S L, GONG C L. The gene expression profile of resistant and susceptible Bombyx mori strains reveals cypovirus-associated variations in host gene transcript levels. Applied Microbiology and Biotechnology, 2015, 99: 5175-5187.

[19] LIU F, LING E, WU S. Gene expression profiling during early response to injury and microbial challenges in the silkworm, Bombyx mori. Archives of Insect Biochemistry and Physiology, 2009, 72(1): 16-33.

[20] YANAI H, SAVITSKY D, TAMURA T, TANIGUCHI T. Regulation of the cytosolic DNA-sensing system in innate immunity: a current view. Current Opinin in Immunology, 2009, 21(1): 17-22.

[21] MANSUR D S, SMITH G L, FERGUSON B J. Intracellular sensing of viral DNA by the innate immune system. Microbes and Infection, 2014, 16(12): 1002-1012.

[22] RATHINAM V A, FITZGERALD K A. Innate immune sensing of DNA viruses. Virology, 2011, 411(2): 153-162.

[23] 邢雅玲, 鄭洋, 王凱, 陳曉娟, 陳忠斌. 病原DNA識別及其誘導天然免疫調節(jié)機制研究進展. 生物化學與生物物理進展, 2011, 38(12): 1099-1105. XING Y L, ZHENG Y, WANG K, CHEN X J, CHEN Z B. The cellular recognition of pathogenic DNA and the related regulation of innate immunity. Progress in Biochemistry and Biophysics, 2011, 38(12): 1099-1105. (in Chinese)

[24] EBIHARA N, CHEN L, TOKURA T, USHIO H, IWATSU M, MURAKAMI A. Distinct functions between toll-like receptors 3 and 9 in retinal pigment epithelial cells. Ophthalmic Research, 2007, 39(3): 155-163.

[25] CHOI M K, WANG Z C, BAN T, YANAI H, LU Y, KOSHIBA R, NAKAIMA Y, HANGAI S, SAVITSKY D, NAKASATO M, NEGISHI H, TAKEUCHI O, HONDA K, AKIRA S, TAMURA T, TANIGUCHI T. A selective contribution of the RIG-I-like receptor pathway to type I interferon responses activated by cytosolic DNA. Proceedings of the National Academy of Sciences of the United States of America, 2009, 106(42): 17870-17875.

[26] CHIU Y H, MACMILLAN J B, CHEN Z J. RNA polymerase III detects cytosolic DNA and induces type I interferons through the RIG-I pathway. Cell, 2009, 138(3): 576-591.

[27] VEERANKI S, CHOUBEY D. Interferon-inducible p200-family protein IFI16, an innate immune sensor for cytosolic and nuclear double-stranded DNA: regulation of subcellular localization. Molecular Immunology, 2012, 49(4): 567-571.

[28] ZHANG Z Q, YUAN B, BAO M S, LU N, KIM T, LIU Y J. The helicase DDX41 senses intracellular DNA mediated by the adaptor STING in dendritic cells. Nature Immunology, 2011, 12(10): 959-965.

[29] KIM T, PAZHOOR S, BAO M S, ZHANG Z Q, HANABUCHI S, FACCHINETTI V, BOVER L, PLUMAS J, CHAPEROT L, QIN J, LIU Y J. Aspartate-glutamate-alanine-histidine box motif (DEAH)/ RNA helicase A helicases sense microbial DNA in human plasmacytoid dendritic cells. Proceedings of the National Academy of Sciences of the United States of America, 2010, 107(34): 15181-15186.

[30] ZHANG X, BRANN T W, ZHOU M, YANG J, OGUARIRI R M, LIDIE K B, IMAMICHI H, HUANG D W, LEMPICKI R A, BASELER M W, VEENSTRA T D, YOUNG H A, LANE H C, IMAMICHI T. Ku70 is a novel cytosolic DNA sensor that induces type III rather than type I IFN. Journal of Immunology, 2011, 186(8): 4541-4545.

[31] ABLASSER A, BAUERNFEIND F, HARTMANN G, LATZ E, FITZGERALD K A, HORNUNG V. RIG-I dependent sensing of poly (dA-dT) through the induction of an RNA polymerase III transcribed RNA intermediate. Nature Immunology, 2009, 10(10): 1065-1072.

[32] MELCHJORSEN J, RINTAHAKA J, S?BY S, HORAN K A, POLTAJAINEN A, ?STERGAARD L, PALUDAN S R, MATIKAINEN S. Early innate recognition of herpes simplex virus in human primary macrophages is mediated via the MDA5/MAVS- dependent and MDA5/MAVS/RNA polymerase III-independent pathways. Journal of Virology, 2010, 84(21): 11350-11358.

[33] ARUNACHALAM B, PHAN U T, GEUZE H J, CRESSWELL P. Enzymatic reduction of disulfide bonds in lysosomes: Characterization of a gamma interferon inducible lysosomal thiol reductase (GILT). Proceedings of the National Academy of Sciences of the United States of America, 2000, 97(2): 745-750.

[34] DE ZOYSA M, LEE J. Molecular cloning and expression analysis of interferon-γ inducible lysosomal thiol reductase (GILT)-like cDNA from disk abalone (Haliotis discus discus). Journal of Invertebrate Pathology, 2007, 96(3): 221-229.

[35] HUANG W S, DUAN L P, HUANG B, ZHOU L H, LIANG Y, TU CL, ZHANG F F, NIE P, WANG T. Identification of three IFN-gamma inducible lysosomal thiol reductase (GILT)-like genes in mud crab Scylla paramamosain with distinct gene organizations and patterns of expression. Gene, 2015, 570: 78-88.

[36] KONGTON K, PHONGDARA A, SRITHAWORN M T, WANNA W. Molecular cloning and expression analysis of the interferon-γinducible lysosomal thiol reductase gene from the shrimp Penaeus monodon. Molecular Biology Reports, 2011, 38: 3463-3470.

[37] KONGTON K, MCCALL K, PHONGDARA A. Identification of gamma-interferon-inducible lysosomal thiol reductase (GILT) homologues in the fruit fly Drosophila melanogaster. Developmental and Comparative Immunology, 2014, 44: 389-396.

[38] HASTINGS K T, CRESSWELL P. Disulfide reduction in the endocytic pathway: immunological functions of gamma-interferoninducible lysosomal thiol reductase. Antioxidants and Redox Signaling, 2011, 15(3): 657-668.

[39] PROKOPENKO O, MIROCHNITCHENKO O. Ischemia-reperfusioninducible protein modulates cell sensitivity to anticancer drugs by regulating activity of efflux transporter. American Journal of Physiology-Cell Physiology, 2009, 296: C1086-C1097.

[40] SHENG Y, SARIDAKIS V, SARKARI F, DUAN S, WU T, ARROWSMITH C H, FRAPPIER L. Molecular recognition of p53 and MDM2 by USP7/HAUSP. Nature Structral and Molecular Biology, 2006, 13(3): 285-291.

[41] HOLOWATY M N, FRAPPIER L. HAUSP/USP7 as an Epstein-Barr virus target. Biochemical Society Transactions, 2004, 32(5): 731-732.

（責任編輯 岳梅）

Digital Gene Expression Analysis of Silkworm Infected by Bombyx mori Bidensovirus Zhenjiang Strain

GAO Kun, SHANG Meng-ke, QIAN He-ying, QIN Guang-xing, GUO Xi-jie

(College of Biotechnology, Jiangsu University of Science and Technology/Sericultural Research Institute, Chinese Academy of Agricultural Sciences, Zhenjiang 212018, Jiangsu)

【Objective】 The objective of this study is to screen differentially expressed genes in the Bombyx mori larvae infected with BmBDV-ZJ (B. mori bidensovirus Zhenjiang strain) and identify regulatory genes related to the virus infection and the host response so as to provide important clues for better understanding of the mechanism of B. mori resistance against BmBDV-ZJ infection. 【Method】 The differential gene expression profiles in JS B. mori larvae after oral infection with BmBDV-ZJ were constructed using Illumina Genome Analyzer platform. In order to exclude the effects of individual differences, 10 larvae were dissected and pooled as one sample for digital gene expression (DGE) analysis, respectively. The differential expression detection of genes across samples was performed using a rigorous algorithm method. False discovery rate (FDR) was used to determine the P value threshold in multiple tests and analyses. The significance of the gene expression difference was obtained through a FDR≤0.001 and the absolute value of log2ratio≥1. The gene ontology (GO) classification system was used to determine the possible functions of all differentially expressed genes. P value was calculated by GO (http://www.geneontology. org/) and corrected by Bonferoni. A corrected P value≤ 0.05 was selected as a threshold for significant enrichment of the gene sets. WEGO (web gene ontology annotation plot) software was used for visualizing, comparing and plotting GO annotation results. Pathway enrichment analysis was conducted to further identify the significantly enriched metabolic pathways or signal transduction pathways by using the KEGG database. Pathways with a Q value≤0.05 were designated as significantly enriched pathways in DGEs. Then some of the differentially expressed genes were verified by quantitative real-time PCR (qRT-PCR). 【Result】 Totally, 4 850 663 and 4 875 307 raw tags were generated in the control and BmBDV-ZJ infected DGE (digital gene expression) libraries, respectively. There were 4 757 934 and 4 788 406 clean tags corresponding to 62 436 and 63 680 distinct clean tags were filtered from the raw tags. The distribution of the total and distinct tags over the different tag abundance categories showed highly similar patterns in each DGE library. The sequencing depths reached approximately 3.5 and 3.7 million in the two DGE libraries, respectively, which satisfied the requirement for the experiment. So the two DGE libraries were reliable. The tag sequences of the two DGE libraries were mapped to the reference database of B. mori. In the control and BmBDV-ZJ-infected DGE library, 36.39% and 45.30% of the clean tags were mapped to a gene in the reference database, 50.02% and 43.34% of the clean tags could be mapped to genome of B. mori, while 13.59% and 12.35% of the clean tags were unknown tags. A total of 447 differentially expressed genes were detected, of which 306 were upregulated and 141 were downregulated. There were 218, 147, 179 differentially expressed genes have GO categories according to molecular function, cellular component and biological process, respectively. KEGG (http://www.genome.jp/kegg) ontology assignments were used to classify the functional annotations of the identified genes. Among the differentially expressed genes, 330 were mapped to 151 pathways in the KEGG database. Nineteen terms was significantly enriched (Q value≤0.05) and the cytosolic DNA-sensing pathway was significantly enriched. Moreover, 24 differentially expressed genes were verified using qRT-PCR, showing that 20 genes were concordant in the expression with DGE. Among the 9 differentially expressed genes related to cytosolic DNA-sensing pathway, BGIBMGA009408-TA, BGIBMGA004913-TA, BGIBMGA011753-TA, which were the DNA-directed RNA polymerase III genes in B. mori, were all up-regulated in the BmBDV-ZJ infected B. mori with 4.3, 2.3, 3.4-fold change, respectively. 【Conclusion】 The results of this study may serve as a basis for future research not only on the molecular mechanism of BmBDV-ZJ invasion but also on the mechanism of B. mori resistance against BmBDV-ZJ infection.

Bombyx mori; digital gene expression; Bombyx mori bidensovirus; quantitative real-time PCR

2016-04-08；接受日期：2016-06-02

國家自然科學基金（31402141）、江蘇省自然科學基金（BK20140508）

聯系方式：高坤，E-mail：gkunjn2002@126.com。通信作者郭錫杰，Tel：0511-84401328；E-mail：guoxijie@126.com