楊凡力，王意銀，鄭文成，何彪，江廷磊，李瑩瑩，夏樂樂，馮燁，范泉水，涂長春

1 吉林大學畜牧獸醫(yī)學院，吉林 長春 1300622 軍事醫(yī)學科學院軍事獸醫(yī)研究所，吉林 長春 1301223 成都軍區(qū)疾病預防控制中心，昆明 云南 6501184 湖南省動物衛(wèi)生監(jiān)督所，長沙 湖南 4100075 東北師范大學吉林省動物資源保護與利用重點實驗室，吉林 長春 130024

中國部分地區(qū)蝙蝠攜帶病毒的宏基因組學分析

楊凡力1,2，王意銀3，鄭文成4，何彪2，江廷磊5，李瑩瑩2，夏樂樂2，馮燁2，范泉水3，涂長春2

1 吉林大學畜牧獸醫(yī)學院，吉林 長春 130062
2 軍事醫(yī)學科學院軍事獸醫(yī)研究所，吉林 長春 130122
3 成都軍區(qū)疾病預防控制中心，昆明 云南 650118
4 湖南省動物衛(wèi)生監(jiān)督所，長沙 湖南 410007
5 東北師范大學吉林省動物資源保護與利用重點實驗室，吉林 長春 130024

楊凡力, 王意銀, 鄭文成, 等. 中國部分地區(qū)蝙蝠攜帶病毒的宏基因組學分析. 生物工程學報, 2013, 29(5): 586?600.

Yang FL, Wang YY, Zheng WC, et al. Metagenomic analysis of bat virome in several regions of China. Chin J Biotech, 2013,29(5): 586?600.

蝙蝠攜帶有60多種病毒，其中許多對人有高度致病性。為了解中國蝙蝠攜帶病毒的自然本底、蝙蝠病毒的多樣性和挖掘潛在的病毒病原，通過基于 Solexa高通量測序的病毒宏基因組學技術對從吉林、云南、湖南采集的蝙蝠組織進行病毒組學研究，獲得了 11 644 232條讀長 (Reads)，并拼接出 44 872條重疊序列(Contig)。通過核酸序列注釋發(fā)現(xiàn)，其中8.2% (4 002/44 872) 的重疊序列與病毒相關，能進一步注釋到36個病毒科，包括19種脊椎動物病毒、6種植物病毒、4種昆蟲病毒和4種噬菌體。通過對重疊序列的遺傳進化分析、多序列比對顯示，被注釋為細小病毒、腺聯(lián)病毒、博卡病毒、腺病毒、小雙節(jié)RNA病毒等的重疊序列與已知病毒相似，部分序列卻又呈現(xiàn)出明顯的序列差異。通過對腺病毒和博卡病毒進一步的 PCR擴增證實了此研究方法可靠。旨在了解我國蝙蝠攜帶病毒組的構成，對建立高效的野生動物源人獸共患病的監(jiān)測方法提供參考。

病毒宏基因組學，蝙蝠，病毒組，高通量測序

蝙蝠屬于翼手目 (Chiroptera)，是僅次于嚙齒類的第二大類哺乳動物，占全球哺乳動物種類的20%，且分布廣泛，除南北極外全球各個地區(qū)均有蝙蝠的存在[1]。近年來在蝙蝠體內發(fā)現(xiàn)的病毒逐漸增多，至 2007年為止，至少從蝙蝠體內檢測或分離到 60多種病毒，其中包括一些致病性病毒，如埃博拉病毒 (Ebola virus)、亨德拉尼帕病毒 (Henipaviruses)、馬爾堡病毒 (Marburg virus)、SARS冠狀病毒 (SARS coronavirus)、狂犬病病毒 (Rabies virus) 以及狂犬病相關病毒(Lyssavirus) 等[2]。最近，Bokeloh和 Shimoni蝙蝠病毒、圓環(huán)病毒 (Circovirus)、博卡病毒(Bocavirus)、呼腸孤病毒 (Retrovirus)、星狀病毒(Astrovirus) 和松灣病毒 (Cedar virus) 作為新病毒在蝙蝠體內發(fā)現(xiàn)[3-8]。2002年SARS冠狀病毒疫情爆發(fā)，造成全球8 000多人感染，至少800人死亡[9]。研究表明我國的菊頭蝠攜帶有 SRAS樣冠狀病毒 (SARS-like coronavirus)，很可能是SARS樣冠狀病毒的自然宿主[10]。近期研究者又在加納和歐洲的蝙蝠體內發(fā)現(xiàn)了與 2012年在中東地區(qū)引起2人死亡的人β冠狀病毒2cEMC/2012(Human betacoronavirus 2cEMC/2012) 接近的 β冠狀病毒[11]。因此，蝙蝠作為巨大的病毒貯存庫越來越受到關注。我國自上世紀 80年代開展蝙蝠病毒研究以來，至少從蝙蝠體內檢測或分離到16種病毒，包括基孔肯雅病毒 (Chikungunya virus)、羅斯河病毒 (Ross River virus)、乙型腦炎病毒 (Japanese encephalitis virus，JEV)、呼腸孤病毒和腺病毒 (Adenovirus)[6,12-15]。

我國蝙蝠有7科30屬120種，分布廣泛，多集中在華中、華南、西南等人口密集且氣候濕潤溫熱的地方[16]。因此應全面掌握我國蝙蝠攜帶病毒的自然本底以預防新發(fā)和再發(fā)人獸共患病。傳統(tǒng)的病毒病原生態(tài)學研究方法，如細胞分離、核酸檢測、血清學試驗等都有一定的局限性。近年來基于高通量測序的病毒宏基因組學方法已經(jīng)成為一個高效的工具用于對自然環(huán)境[17-18]、人[19]和動物[20-21]的病毒組學分析。目前采用病毒宏基因組學研究蝙蝠病毒在全球已有 4則報道[22-25]，他們對蝙蝠糞便的病毒組研究在一定程度上反映了蝙蝠攜帶病毒的多樣性，提供了有價值的數(shù)據(jù)，但目前對蝙蝠病毒的研究并不充分。為了建立蝙蝠病毒宏基因組學研究方法、豐富我國蝙蝠攜帶病毒的數(shù)據(jù)庫，本研究自 2008年至2010年從我國三地采集到4種共241只蝙蝠樣品，通過病毒宏基因組學對這些蝙蝠的腸道和肺臟組織進行病毒組學分析，發(fā)現(xiàn)了細小病毒(Parvovirus)、博卡病毒、小雙節(jié) RNA病毒(Picobirnavirus)、腺病毒等病毒的序列，描繪了上述蝙蝠組織樣本攜帶病毒的情況，并通過對博卡病毒和腺病毒的 PCR驗證，進一步證實了這些病毒的存在，為蝙蝠病毒的病原生態(tài)學研究提供一定的參考。

1 材料與方法

1.1 樣品采集

2008年至2010年共采集了4種241只蝙蝠。2008年 12月在云南省靠近緬甸邊境的拉波寨采集到 180只亞洲長翼蝠Miniopterus fuliginosus；2010年9月在吉林省長春市一座公寓頂層采集到31只薩氏伏翼Hypsugo savii和4只東方蝙蝠Vespertilio sinensis；2010年12月在湖南省慈利縣采集到 26只大蹄蝠Hipposideros armiger。經(jīng)形態(tài)學鑒定蝙蝠種類后立即解剖，將腸道 (含內容物)、肺臟分別取出裝入2 mL凍存管，保存于?80 ℃。此研究嚴格按照軍事醫(yī)學科學院軍事獸醫(yī)研究所動物福利委員會相關規(guī)定執(zhí)行。

1.2 試劑及材料

反 轉 錄 酶 (SuperScriptⅢ Reverse Transcriptase)、DNA 聚合酶 (AccuprimeTaqDNA Polymerase)、RNA酶抑制劑 (RNAseOUT)購自Invitrogen公司；核酸外切酶 (Exonuclease I)、堿性磷酸酶 (Shrimp alkaline phosphatase)、TRIzol (Total RNA 提取試劑) 購自寶生物 (大連) 公司；QIAquick PCR Purification Kit購自QIAGEN 公司；DNA 酶 (Turbo DNase) 購自Ambion公司；核酸酶 (Benzonase Nuclease) 購自Novagen公司；DNA聚合酶Klenow片段購自NEB公司；0.22 μm針頭濾器和PelliconⅡ膜包購自Millipore公司；SM緩沖液 (50 mmol/L Tris，10 mmol/L MgSO4，0.1 mol/L NaCl，pH 7.5)；DNA提取試劑盒購自 QIAGEN公司；PCR反應液購自天根生物公司。

1.3 樣品制備

在生物二級安全柜內將吉林采集的 35只蝙蝠和湖南采集的 26只蝙蝠的肺臟和腸道樣本每只各取大約0.1 g，云南采集的180只蝙蝠隨機挑選出 40只蝙蝠的肺臟和腸道各取大約 0.1 g混合。不同采集地點混合成1組，共3組 (云南組名稱為 pool-1；湖南組為 pool-2；吉林組為pool-3)。每一個組按1：10 (W/V) 加入無菌SM緩沖液，用 Waring組織研磨儀磨碎；研磨液于4 ℃、10 000 r/min離心30 min；上清液先后用0.45 μm和0.22 μm孔徑的PelliconⅡ膜包濾掉組織碎片、細菌和其他雜質；濾液使用截留分子量為100 kDa的PelliconⅡ濃縮膜包進行濃縮；濃縮液用 Beckman 超速離心機 SW55Ti轉頭45 000 r/min離心2 h；沉淀用SM緩沖液懸浮并使用0.22 μm針頭濾器過濾；為降低游離核酸污染，每 116 μL的濾過液添加 14 U的 DNA酶(Turbo DNase)、25 U 的核酸酶 (Benzonase Nuclease)、20 U的核酸外切酶 (RNase I) 以及10×Turbo DNase緩沖液在37 ℃消化1 h；隨后用TRIzol提取核酸 (包括RNA和DNA)，并溶解于20 μL的無RNA酶的水中。

1.4 反轉錄及隨機PCR

每組分別將8 μL上述病毒核酸加入1 μL二甲基亞砜 (DMSO) 和 1 μL 50 μmol/L 的帶有20 bp錨定序列的6mer隨機引物 (pool-1: 5￠-GC CGGAGCTCTGCAGATATCNNNNNN-3￠；pool-2:5￠-GTATCGCTGGACACTGGACCNNNNNN-3￠；pool-3: 5￠-CGCATTGGTCGGCACTTGGTNNNN NN-3￠)。于 72 ℃作用 5 min，然后立即冰浴 2 min；加入 1 μL DTT、1 μL dNTP mixture (10 mmol/L)、20 U RNAseOUT、4 μL 5×緩沖液、200 U Superscript Ⅲ和 ddH2O，25 ℃反應 10 min，50 ℃反應50 min合成cDNA，85 ℃ 1 0 min 滅活反轉錄酶。加入5 U DNA聚合酶Klenow片段37 ℃反應60 min合成雙鏈cDNA并75 ℃滅活10 min，隨后在此體系中加入 2 U堿性磷酸酶 (Shrimp alkaline phosphatase) 和 2.5 U 核酸外切酶(Exonuclease I) 于37 ℃反應60 min，除去多余的引物和游離的核苷酸，并75 ℃滅活10 min。然后取 10 μL 上述模板，加入 2 μL 10 μmol/L 錨定序列引物 (pool-1: 5￠-GCCGGAGCTCTGCAGAT ATC-3￠；pool-2: 5￠-GTATCGCTGGACACTGGA CC-3￠；pool-3: 5￠-CGCATTGGTCGGCACTTGG T-3￠)、5 μL 10×AccuPrime bufferⅠ、1 μL DNA 聚合酶 (AccuprimeTaqDNA Polymerase) 及ddH2O，共50 μL體系進行不依賴序列的單引物擴增 (Sequence-independent single primer amplification，SISPA)：94 ℃ 30 s，55 ℃ 30 s，68 ℃ 1 min ，35個循環(huán)。擴增產(chǎn)物使用QIAquick PCR Purification Kit試劑盒純化并溶于50 μL的TE緩沖液中。

1.5 Solexa高通量測序

純化后的PCR產(chǎn)物送華大基因公司 (深圳)進行Solexa高通量測序。將3組樣品混合后使用超聲法打斷到180 bp左右，末端分別加上接頭，通過接頭 PCR將 DNA片段連接到芯片上形成“橋”(Bridge) 并制成文庫，隨后應用4種熒光標記的核苷酸通過橋式PCR (Bridge PCR) 進行邊合成邊測序 (Sequencing By Synthesis，SBS)，獲得原始讀長數(shù)據(jù)，去掉長度<100 bp的讀長和錨定引物序列，通過SOAPdenove軟件拼接成重疊序列 (Contig)，隨后進行數(shù)據(jù)庫比對分析。首先將獲得的重疊序列使用 BLASTn和 BLASTx在 GenBank的非冗余核苷酸序列數(shù)據(jù)庫(Nonredundant database) 中進行比對，將E value ＜10?3的序列作為有意義的序列，然后按照注釋到的物種信息去掉細菌及真核類序列，余下的病毒相關序列通過系統(tǒng)進化分析和多序列比對進行鑒定。

1.6 病毒樣序列分析

根據(jù)BLAST注釋信息，挑選E value較小且同源性較高的病毒樣重疊序列使用 DNAStar編輯，整理后用 MEGA5.0[26]軟件進行系統(tǒng)發(fā)生分析，進化樹選擇基于 Maximum Composite Likelihood算法的鄰接法 (Neighbor-joining method，NJ)，Bootstrap檢驗設定為1 000個重復。使用DNAstar軟件包的MegAlign軟件對部分低同源性的病毒樣序列和 GenBank中參考序列的氨基酸序列進行多序列比對，從而進行人工輔助鑒定。

1.7 PCR檢測驗證博卡病毒和腺病毒

使用博卡病毒和腺病毒序列信息，分別設計了套式 PCR對這兩種病毒進行驗證。通過GenBank下載博卡病毒參考序列，使用 Primer Premier 5.0軟件設計檢測博卡病毒的簡并引物。腺病毒的檢測引物參考已發(fā)表的文獻[15]。使用QIAGEN的DNA提取試劑盒對保存的蝙蝠組織樣品 (腸道和肺臟) 單樣提取 DNA，DNA提取后進行PCR檢測。PCR反應體系為：PCR反應液 (天根) 25 μL、ddH2O 19 μL、上游引物(10 pmol/μL) 2 μL、下游引物 (10 pmol/μL) 2 μL、DNA模板2 μL。PCR溫度條件為：94 ℃預變性2 min；然后94 30 s℃、53 30 s℃、72 40 s℃，外套反應和內套反應各 35個循環(huán)；72 ℃延伸7 min。PCR產(chǎn)物經(jīng)電泳鑒定為陽性后，送吉林庫美基因公司測序。將獲得的序列構建系統(tǒng)進化樹，方法與重疊序列建樹方法一致。

2 結果

2.1 測序DNA樣品的制備

經(jīng)過病毒純化、核酸提取后使用ND1000核酸蛋白紫外定量檢測儀對3組樣品的核酸濃度進行定量。所提取 RNA和 DNA總濃度均高于200 ng/μL。進過反轉錄和不依賴序列的單引物擴增，產(chǎn)物條帶為抹帶，分布于2 000 bp以內，主要集中在100 bp至500 bp的區(qū)域 (圖1)。3組的PCR產(chǎn)物經(jīng)純化回收后經(jīng)ND1000核酸蛋白紫外定量檢測儀測定 260/280吸光值比在 1.8至 2.0之間，樣品濃度＞800 ng/μL。符合Solexa測序的要求。

2.2 Solexa測序結果

圖1 不依賴序列的單引物擴增結果Fig. 1 Sequence-independent single primer amplification results of viral nucleic acids. M: DNA marker.

通過高通量測序，獲得了11 644 232個讀長，并進一步拼接成48 872個重疊序列，平均長度為136.7 bp。與NCBI病毒數(shù)據(jù)庫比對，4 002條重疊序列注釋為病毒，占所有重疊序列的 8.2%，其中2 310條可注釋到36個不同的病毒科，另1 692條重疊序列被定義為疑似病毒序列但未被分類。在這36個科中，44.6% (1 030/2 310) 的序列與 19個脊椎動物病毒科相關，包括腺病毒科 (Adenoviridae)、皰疹病毒科 (Herpesviridae)、乳頭瘤病毒科 (Papillomaviridae)、細小病毒科(Parvoviridae)、圓環(huán)病毒科 (Circoviridae)、冠狀病毒科 (Coronaviridae)、小雙節(jié) RNA 病毒科(Picobirnaviridae)、 小 RNA 病 毒 科(Picornaviridae)等；昆蟲病毒相關序列占到了6.9% (160/2310)，包括囊泡病毒科(Ascoviridae)、桿狀病毒科 (Baculoviridae)、多DNA病毒科 (Polydaviridae)、布尼亞病毒科(Bunyaviridae) 4個科；植物病毒序列占 1.8%(41/2310) 包括花椰菜花葉病毒科(Caulimoviridae) 等 6個科；另外，還有43.0%(993/2310) 的序列與噬菌體相關，包括長尾噬菌體科 (Siphoviridae)、短尾噬菌體科(Podoviridae) 等4個科 (圖2)。在脊椎動物病毒序列中，腺病毒、圓環(huán)病毒、腺聯(lián)病毒(Adeno-associated viruses，AAV)、博卡病毒、細小病毒、小雙節(jié)RNA病毒的序列在GenBank中能找到較好的比對結果，即 E value＜10?20，比對得分＞100，且核苷酸或氨基酸序列與已知序列同源性在 60%至 98%，經(jīng)遺傳進化分析和分子生物學鑒定后可以確認為已知病毒和新發(fā)現(xiàn)的病毒序列；而皰疹病毒、乳頭瘤病毒、冠狀病毒等相關序列的很大一部分無論是核酸序列還是氨基酸序列與現(xiàn)有序列同源性低于 60%，且E value在 10?2～10?5之間，目前可將其定義為新病毒相關序列。

圖2 Solexa高通量測序獲得的重疊序列注釋Fig. 2 Annotation of Solexa contigs. The numbers and percentages of each virus species were labeld on left. The figure on right is the number of the contigs, which are classified into 36 viral families.

2.3 細小病毒 (Parvovirus，PV) 重疊序列鑒定

本研究有 2條重疊序列注釋到細小病毒屬(Parvovirus)。其中一條為291 nt的NS1基因片段。通過與GenBank數(shù)據(jù)庫BLASTx比對，發(fā)現(xiàn)與之同源性最高的為小鼠微小病毒 (Mice minute virus)，局部有70%的氨基酸同源性；另一條重疊序列為114 nt的VP1基因片段，與小鼠微小病毒有74%的氨基酸序列同源性。使用相對較長的 291 nt的重疊序列 (BtPV-0866) 構建的系統(tǒng)進化樹。選取細小病毒屬的代表序列作為參考序列的同時用腺聯(lián)病毒 2型 (AAV-2) 作為外群對照。進化樹顯示BtPV-0866處在細小病毒屬的外圍，形成一個獨立分支 (圖 3A)。根據(jù)核苷酸序列推導的氨基酸序列長度為95個殘基，與小鼠微小病毒該區(qū)域氨基酸同源性為56% (圖3B)。利用氨基酸序列進行多序列比對發(fā)現(xiàn)，BtPV-0866在保守基序處與參考序列完全相同，如N端第6～13位的DQKGKGSK、第31～48位僅有2個殘基的區(qū)別，在第66～80位區(qū)域變異較大，與參考序列僅有一個殘基相同；而作為外群對照的 AAV-2氨基酸序列僅在少數(shù)位點與細小病毒屬基序相同。根據(jù)核苷酸序列和氨基酸序列的鑒定，認為所獲得的重疊序列來自新的細小病毒。

圖3 細小病毒重疊序列鑒定Fig. 3 Characterization of PV-like contigs. (A) Phylogenetical analysis of PV-like contig BtPV-0866 with reference sequences based on about 290 nt NS1 gene. (B) Multiple sequences alignment of deduced amino acid sequences of parvovirus NS1. The black boxed aa sequence is majority sequence.

2.4 腺聯(lián)病毒 (Adeno-associated viruses，AAV) 重疊序列鑒定

之前的研究已在棕果蝠、菊頭蝠、鼠耳蝠等蝙蝠糞便內檢測到AAV。本研究中共11條重疊序列被注釋到依賴病毒屬 (Dependovirus)，其中7條 (長度在108～318 nt之間) 與在云南省的大足鼠耳蝠糞便中發(fā)現(xiàn)的蝙蝠腺聯(lián)病毒 YNM 株(BtAAV-YNM)[27]在6個不同的區(qū)域有78%～95%的核苷酸同源性 (圖4A)。如168 nt的重疊序列BtAAV-0854與BtAAV-YNM的Rep基因具有高達95%的核苷酸同源性，而192 nt的BtAAV-2413與BtAAV-YNM的Cap基因具有高達90%的核苷酸同源性。其余4條重疊序列與GenBank數(shù)據(jù)庫中的AAV序列有46%～60%的核苷酸同源性，說明還存在其他與已知的 AAV遺傳距離較大的病毒。對裁剪后 150 nt的Cap基因序列BtAAV-2143和BtAAV-3759的遺傳進化分析顯示BtAAV-2143與BtAdV-YNM處于同一分支，遺傳距離較近；而BtAAV-3759獨立成為一個分支，與AAV代表株遺傳距離較大，屬于依賴病毒屬一個新的種 (圖4B)。綜合對11條疑似AAV重疊序列的分析，雖無法確定這 11條重疊序列來自多少株病毒，但可以確定至少有 1株與BtAAV-YNM同源性較高的AAV和1株與已知的病毒株變異較大的AAV存在這批蝙蝠樣本中。

圖4 腺聯(lián)病毒重疊序列鑒定Fig. 4 Characterization of AAV-like contigs. (A) Sequence comparison of AAV-like contigs with BtAAV-YNM. The percentages above bars are similarity of contigs with BtAAV-YNM. (B) Phylogenetical analysis of AAV-like contigs with reference sequences based on 150 nt Cap gene.

2.5 小雙節(jié)RNA病毒 (Picobirnavirus，PBV)重疊序列鑒定

本研究獲得的3條PBV相關的重疊序列，1條位于第一節(jié)段與兔小雙節(jié) RNA病毒 (Rabbit picobirnavirus) 具有50%的核苷酸同源性。另外兩條位于第二節(jié)段，與人小雙節(jié) RNA病毒(Human PBV GPBV11) 分別有61%和63%的核苷酸同源性。通過重疊序列BtPBV-1798與其他來源的PBV第二節(jié)段 (編碼RdRp) 303 bp局部序列的進化分析顯示出與現(xiàn)有的 PBV中的Human PBV GPBV11相對距離較近，處于一個進化分支上 (圖5A)。將核苷酸序列翻譯為97個殘基，與Human PBV GPBV11同源性為55%。多序列比對顯示該殘基序列與參考序列具有共同的保守殘基，如位于序列 N端的 QEGGP、DDVKQR和序列中部的AVD、FYQP、TDFTK、HFN等保守基序。而在氨基酸序列變異較大的區(qū)域，如37至49位和第84至96位殘基，重疊序列以及其他參考序列相互之間均有較大的變異(圖 5B)。

圖5 小雙節(jié)RNA病毒重疊序列鑒定Fig. 5 Characterization of PBV-like contig. (A) Phylegenetical analysis of PBV-like contigs with reference sequences based on 303 nt RdRp gene. (B) Multiple sequence alignment of deduced amino acid sequences of PBV. The black boxed aa sequence is majority sequence. The motif QEGGP, DDVKQR, FYQP, TDFTK and HFN are conservative.

2.6 博卡病毒 (bocavirus，BoCV) 重疊序列鑒定及PCR驗證

本研究發(fā)現(xiàn)了一條重疊序列注釋到細小病毒科的博卡病毒屬 (Bocavirus)。BLAST比對發(fā)現(xiàn)這條序列位于博卡病毒NP1基因，與 CnMV核苷酸同源性最高，為62%。通過PCR擴增博卡病毒588 nt的VP1部分片段，在云南的亞洲長翼蝠腸道樣品中檢測出 2份陽性樣品 (圖 6A)。經(jīng)測序兩樣品之間同源性100%，為同一毒株，暫命名為 Bat bocavirus YNMi-1 (GenBank Accession No. KC172378)。該序列與 CnMV有55%的氨基酸同源性。PCR擴增序列與參考序列進行系統(tǒng)進化分析顯示所獲得的 Bat bocavirus YNMi-1序列與CnMV、海獅博卡病毒 (Sea lion BoCV)、貓博卡病毒 (FBoCV) 處于同一個進化分支，但在這個進化分支上 Bat bocavirus YNMi-1和另外3株病毒有相對較大進化距離。YNMi-1與此前在海南鼠耳蝠體內發(fā)現(xiàn)的1株蝙蝠博卡病毒 (Myotis BoCV-1) 處在不同的進化分支且僅有56%的核苷酸同源性 (圖6B)，說明兩株蝙蝠博卡病毒無直接的進化關系即蝙蝠博卡病毒具有基因多樣性。利用翻譯的氨基酸序列與參考序列進行多序列比對顯示，本研究所擴增的 DNA片段序列中部在標尺的第 46～68位(HVR1)、84～133 位 (HVR2)、179～194位 (HVR3)存在3個高度變異區(qū)，序列之間無論是氨基酸的種類還是長度都有很大差異。YNMi-1在 HVR1區(qū)域相比Sea lion BoCV和FBoCV缺失了6個氨基酸殘基；在HVR2和HVR3區(qū)域，YNMi-1相對于主要序列分別額外插入了8個和5個氨基酸殘基 (圖6C)。YNMi-1在中部的保守基序如標尺標記的69～82位與主序列僅有1個殘基的差別，在137～147位殘基完全相同。

2.7 腺病毒 (Adenovirus，AdV) 重疊序列鑒定及PCR驗證

本研究有 69條重疊序列注釋為哺乳動物腺病毒屬 (Mastadenovirus)。BLAST比對發(fā)現(xiàn)這些序列都屬于哺乳動物腺病毒的hexon、polymerase、penton、pX、pVI等基因。其中3條重疊序列與 GenBank中西班牙的一株腺病毒Has070613-2pol基因或hexon基因有分別高達98%，96%和91%的核苷酸同源性。體現(xiàn)了同一蝙蝠腺病毒毒株分布的廣泛性。其余重疊序列與已知的腺病毒序列同源性均低于83%。

通過套式PCR擴增腺病毒270 ntpol基因檢測，從湖南的樣品中檢測出一份腺病毒陽性樣品(圖 7A)，暫命名為 bat adenovirus HNHi-1(GenBank Accession No. KC171020)。經(jīng)比對與犬腺病毒 (CAdV)和豬腺病毒 A (PAdV-A) 具有71%的核苷酸同源性，與蝙蝠腺病毒 A(BtAdV-A) 和蝙蝠腺病毒 B (BtAdV-B) 具有70%的核苷酸同源性。系統(tǒng)進化分析顯示BtAdV-HNHi-1屬于哺乳動物腺病毒與 CAdV、BtAdV-A、BtAdV-B、PAdV-A處在同一個大的分支，相對與犬腺病毒距離最近 (圖 7B)。由于擴增片段較短，新發(fā)現(xiàn)的蝙蝠腺病毒與其他蝙蝠腺病毒的進化關系還有待進一步研究。

圖6 博卡病毒PCR擴增和序列鑒定Fig. 6 Characterization of bocavirus amplicon. (A) The result of bocavirus PCR amplification. M: DNA marker DL2000; lane 1 to 19 represent samples; lane 20: ddH2O. (B) About 585 nt VP1 gene-based phylogenetical analysis of bocavirus PCR amplicon with reference. (C) Multiple sequences alignment of deduced amino acid sequences of bocavirus. The black boxed aa sequence is majority sequence.

圖7 腺病毒PCR擴增和序列鑒定Fig. 7 Characterization of adenovirus amplicon. (A) The result of adenovirus PCR amplification. (B) Phylegenetical analysis of adenovirus PCR amplicon compared with reference sequences based on about 270 nt hexon gene.

3 討論

隨著人類活動區(qū)域的不斷擴大，蝙蝠的生存空間遭到侵占，創(chuàng)造了更多蝙蝠與人類接觸的機會，使得許多高致病性的病毒由蝙蝠傳播給人類或家畜。發(fā)現(xiàn)蝙蝠攜帶的高致病性病毒和新型病毒、監(jiān)控蝙蝠病毒變異，已成為人們進行蝙蝠病毒流行病學研究的重要內容。病毒宏基因組學技術在這一領域的研究中效率最高、覆蓋的病毒面最廣。目前病毒宏基因組學不僅已用在水體、土壤等無機環(huán)境，也在人、浣熊、火雞、海龜?shù)葎游飻y帶病毒的研究中展現(xiàn)出明顯優(yōu)勢，其獲得的結果系統(tǒng)地描繪了它們的病毒組[28]。本研究通過建立組織病毒宏基因組學方法，對蝙蝠的呼吸系統(tǒng)和消化系統(tǒng) (包括內容物) 兩類組織進行病毒組學研究。通過Solxa測序和生物信息學注釋獲得的病毒相關重疊序列占總數(shù)的8.2%，而之前的病毒宏基因組研究病毒序列占總序列數(shù)的0.1%～58%不等[22-25]。并且注釋為病毒的豐度與前面4份蝙蝠病毒宏基因組研究結果也有差異，首先在本研究中噬菌體占了43.0%，屬于僅次于脊椎動物病毒的第二大類的病毒，在以前美國和我國的一例研究結果顯示噬菌體的序列在多個混樣 (pool) 中數(shù)量居第一或第二[22,24]；其次，脊椎動物病毒所占比例最大，高達44.6%，而之前美國和我國的蝙蝠糞便病毒組研究結果顯示的脊椎動物病毒含量相對較少，甚至低于10%。而本研究中的昆蟲病毒和植物病毒序列數(shù)量較少，而在之前的中國蝙蝠糞便的病毒組研究中昆蟲病毒占第二位置[24]。這些區(qū)別可能是因為我們所處理的腸道和肺臟組織相對較多而糞便含量少，而之前的研究集中在蝙蝠排出的糞便，而昆蟲病毒和植物病毒恰恰只存在于蝙蝠的糞便中，并不在蝙蝠體內細胞中增殖。

本研究嘗試利用加在隨機引物前端的錨定引物序列的差異對湖南、云南、吉林三地的蝙蝠病毒進行分組，但由于在測序時使用超聲法將DNA打斷到 180 bp左右，我們所擴增的DNA產(chǎn)物長度大于500 bp的DNA片段被打斷后，隨機 PCR產(chǎn)物中間部分序列丟失了錨定引物序列信息，無法完成分組，故本次高通量測序結果做數(shù)據(jù)處理時放棄了對蝙蝠病毒進行地域上的分組。但病毒株的地理信息我們仍然可以從后續(xù)的PCR驗證中獲得。在以后的研究中，可以通過凝膠回收300 bp以內的隨機PCR產(chǎn)物進行高通量測序來解決這一問題，從而對不同組樣品的測序結果進行區(qū)分。

皰疹病毒 (Herpesvirues)、乳頭瘤病毒(Papillomavirues) 等相關序列與已知病毒有著低于60%甚至更低的同源性，這些低同源性的序列可能源自一些我們尚不知道的病毒，相信隨著病毒數(shù)據(jù)庫的不斷完善，這些類似病毒序列將被注釋到。還有其他一些重疊序列，如牛病毒性腹瀉病病毒 (BVDV)，所獲得重疊序列雖然與BVDV具有高達90%以上的同源性，但同時與動物染色體上的一段基因序列也有著很高的同源性。這可能反映了病毒基因能整合到某些宿主基因組中，這也是Cui等能從蝙蝠腦組織的轉錄組中發(fā)現(xiàn)反轉錄病毒的原因[29]。

本研究發(fā)現(xiàn)的病毒包括細小病毒、博卡病毒、小雙節(jié)RNA病毒、腺病毒和腺聯(lián)病毒等。博卡病毒屬于細小病毒亞科，博卡病毒屬，目前國際病毒分類命名委員會 (ICTV) 規(guī)定該屬僅有兩個種：牛細小病毒 (Bovine parvovirus，BPV)和犬微小病毒 (Canine minute virus，CnMV)[30]，近年來博卡病毒的宿主范圍從牛、犬擴大到了人、猩猩、豬、海獅、蝙蝠等動物。人博卡病毒能造成兒童的呼吸道感染[31]、胃腸炎和腹瀉[32]；豬博卡病毒則在呼吸道感染的仔豬體內有較高陽性率[33]，因此該病毒具有公共衛(wèi)生意義。小雙節(jié) RNA病毒科目前只有小雙節(jié) RNA病毒屬(Picobirnavirus)，成員為具有兩個節(jié)段的雙鏈RNA病毒，人的小雙節(jié)RNA病毒常在腹瀉病人的糞便中發(fā)現(xiàn)，在家畜、禽類、自然環(huán)境中也能檢測到 PBV[34]。腺病毒科包含 5個屬，均為無囊膜的雙鏈DNA病毒，感染譜廣泛，能在幾乎所有種類的脊椎動物體內發(fā)現(xiàn)[35]。其中人腺病毒的一些血清型能造成呼吸道炎癥和流行性結膜炎[36-37]。本研究通過建立組織病毒宏基因組學方法，對我國三地蝙蝠組織病毒組進行研究，并對一部分結果進行了系統(tǒng)進化分析、多序列比對以及 PCR檢驗。其中一些序列無論核苷酸序列的同源性、序列的遺傳距離還是推導氨基酸序列的特征都支持這些序列可能是新病毒，如細小病毒、博卡病毒、腺病毒、小雙節(jié)RNA病毒、圓環(huán)病毒等。還有一些序列與已知病毒序列有高于90%的同源性，如腺聯(lián)病毒、辣椒溫和斑點病毒、部分腺病毒以及多種噬菌體等。本研究未發(fā)現(xiàn)該地區(qū)蝙蝠攜帶有埃博拉病毒、尼帕病毒、SARS冠狀病毒等已知對人具有高致病性的病毒，但我們不應放松對新發(fā)和再發(fā)蝙蝠源人獸共患病的監(jiān)測。我國蝙蝠分布廣泛、種類眾多，全面地研究我國不同地區(qū)的蝙蝠病毒組為防控蝙蝠將病毒傳播給人類提供了重要信息。

[1]Teeling EC, Springer MS, Madsen O, et al. A molecular phylogeny for bats illuminates biogeography and the fossil record. Science, 2005,307: 580?584.

[2]Calisher CH, Childs JE, Field HE, et al. Bats:important reservoir hosts of emerging viruses. Clin Microbiol Rev, 2006, 19: 531?545.

[3]Freuling CM., Beer M, Conraths FJ, et al. Novel lyssavirus in Natterer's bat, Germany. Emerg Infect Dis, 2011, 17: 1519?1522.

[4]Kuzmin IV, Mayer AE, Niezgoda M, et al. Shimoni bat virus, a new representative of the Lyssavirus genus. Virus Res, 2010, 149: 197?210.

[5]Ge X, Li J, Peng C, et al. Genetic diversity of novel circular ssDNA viruses in bats in China. J Gen Virol, 2011, 92: 2646?2653.

[6]Du L, Lu Z, Fan Y, et al. Xi River virus, a new bat reovirus isolated in southern China. Arch Virol,2010, 155: 1295?1299.

[7]Chu DK, Poon LL, Guan Y, et al. Novel astroviruses in insectivorous bats. J Virol, 2008, 82:9107?9114.

[8]Marsh GA, de Jong C, Barr JA, et al. Cedar virus: a novel Henipavirus isolated from Australian bats.PLoS Pathog, 2012, 8: e1002836.

[9]Shi Z, Hu Z. A review of studies on animal reservoirs of the SARS coronavirus. Virus Res,2008, 133: 74?87.

[10]Li W, Shi Z, Yu M, et al. Bats are natural reservoirs of SARS-like coronaviruses. Science, 2005, 310:676?679.

[11]Augustina A, Heather J, Baldwin, et al. Human betacoronavirus 2c EMC/2012-related Viruses in Bats, Ghana and Europe. Emerg Infect Dis, 2013,19: 456?459.

[12]Zhang HL, Shi HF, Liu LH, et al. Isolation of chikungunya virus from bat in Yunnan province and serological investigations. Chin J Virol, 1989, 5(1):31?36 (in Chinese).

張海林, 施華芳, 劉麗華, 等. 從云南省蝙蝠中分離基孔肯亞病毒及血清抗體調查. 病毒學報,1989, 5(1): 31?36.

[13]Zhao CS, Jiang LH, Yu XL, et al. Isolation of Ross River Virus and its Antibody Prevalance in Hainan Province. Chin J Vet Sci, 1997, 17(3): 241?243 (in Chinese).

趙春生, 蔣廉華, 余興龍, 等. 從海南省蝙蝠腦中分離出1株羅斯河病毒及其血清抗體調查. 中國獸醫(yī)學報, 1997, 17(3): 421?423.

[14]Zhang HL, Zhang YZ, Huang WL, et al. Isolation of Japanese encephalitis virus from brain tissue of bat in Yunnan province. Virol Sin, 2001, 16(1):74?77 (in Chinese).

張海林, 張云智, 黃文麗, 等. 從云南省蝙蝠腦組織中分離出乙型腦炎病毒. 中國病毒學, 2001,16(1): 74?77.

[15]Li Y, Ge X, Zhang H, et al. Host range, prevalence,and genetic diversity of adenoviruses in bats. J Virol, 2010, 84: 3889?3897.

[16]Wang X, Qi D, Hu J. Recent researches on Chiroptera in China. Sichuan J Zool, 2004, 23(2):153?157 (in Chinese).

王曉琴, 齊敦武, 胡錦矗. 中國翼手目研究進展.四川動物, 2004 , 23(2): 153?157.

[17]Djikeng A, Kuzmickas R, Anderson NG, et al.Metagenomic analysis of RNA viruses in a fresh water lake. PLoS ONE, 2009, 4: e7264.

[18]Schoenfeld T, Patterson M, Richardson PM, et al.Assembly of viral metagenomes from yellowstone hot springs. Appl Environ Microbiol, 2008, 74:4164?4174.

[19]Tang P, Chiu C. Metagenomics for the discovery of novel human viruses. Future Microbiol, 2010, 5:177?189.

[20]Lu J, Domingo JS. Turkey fecal microbial community structure and functional gene diversity revealed by 16S rRNA gene and metagenomic sequences. J Microbiol, 2008, 46: 469?477.

[21]Day JM, Ballard LL, Duke MV, et al. Metagenomic analysis of the turkey gut RNA virus community.Virol J, 2010, 7: 313.

[22]Donaldson EF, Haskew AN, Gates JE, et al.Metagenomic analysis of the viromes of three North American bat species: viral diversity among different bat species that share a common habitat. J Virol, 2010, 84: 13004?13018.

[23]Li L, Victoria JG, Wang C, et al. Bat guano virome:predominance of dietary viruses from insects and plants plus novel mammalian viruses. J Virol, 2010,84: 6955?6965.

[24]Ge X, Li Y, Yang X, et al. Metagenomic analysis of viruses from bat fecal samples reveals many novel viruses in insectivorous bats in China. J Virol,2012, 86: 4620?4630.

[25]Wu Z, Ren X, Yang L, et al. Virome analysis for identification of novel Mammalian viruses in bat species from Chinese provinces. J Virol, 2012, 86:10999?11012.

[26]Tamura K, Peterson D, Peterson N, et al. MEGA5:molecular evolutionary genetics analysis using maximum likelihood, evolutionary distance, and maximum parsimony methods. Mol Biol Evol,2011, 28: 2731?2739.

[27]Li Y, Ge X, Hon CC, et al. Prevalence and genetic diversity of adeno-associated viruses in bats from China. J Gen Virol, 2010, 91: 2601?2609.

[28]He B, Tu CC. The advances and applications of viral metagenomics. Acta Vet et Zoo Sin, 2012,43(12): 1865?1870 (in Chinese).

何彪, 涂長春. 病毒宏基因組學的研究現(xiàn)狀及應用. 畜牧獸醫(yī)學報, 2012, 43(12): 1865?1870.

[29]Cui J, Tachedjian M, Wang L, et al. Discovery of retroviral homologs in bats: implications for the origin of mammalian gammaretroviruses. J Virol,2012: 86: 4288?4293.

[30]International Committee on Taxonomy of Viruses.Master, Species, List 2011 v2. [EB/OL]. [2012-12-15]http://ictvonline.org/virusTaxonomy.asp?bhcp=1

[31]Allander T, Tammi MT, Eriksson M, et al. Cloning of a human parvovirus by molecular screening of respiratory tract samples. Proc Natl Acad Sci USA,2005, 102: 12891?12896.

[32]Campe H, Hartberger C, Sing A. Role of human bocavirus infections in outbreaks of gastroenteritis.J Clin Virol, 2008, 43: 340?342.

[33]Zhai S, Yue C, Wei Z, et al. High prevalence of a novel porcine bocavirus in weanling piglets with respiratory tract symptoms in China. Arch Virol,2010, 155: 1313?1317.

[34]Ganesh B, Banyai K, Martella V, et al.Picobirnavirus infections: viral persistence and zoonotic potential. Rev Med Virol, 2012, 22:245?256.

[35]Benko M, Harrach B. Molecular evolution of adenoviruses. Curr Top Microbiol Immunol, 2003,272: 3?35.

[36]McNeill KM, Ridgely Benton F, Monteith SC, et al. Epidemic spread of adenovirus type 4-associated acute respiratory disease between U.S. Army installations. Emerg Infect Dis, 2000, 6: 415?419.

[37]Jawetz E. The story of shipyard eye. Br Med J,1959, 1: 873?876.

December 19, 2012; Accepted: March 5, 2013

Changchun Tu. Tel/Fax: +86-431-86985862; E-mail: changchun_tu@hotmail.com NSFC云南省聯(lián)合基金 (No. U1036601) 資助。

Metagenomic analysis of bat virome in several Chinese regions

Fanli Yang1,2, Yiyin Wang3, Wencheng Zheng4, Biao He2, Tinglei Jiang5, Yingying Li2,Lele Xia2, Ye Feng2, Quanshui Fan3, and Changchun Tu2

1College of Animal Science and Veterinary Medicine,Jilin University,Changchun130062,Jilin,China
2Institute of Military Veterinary,Academy of Military Medical Sciences,Changchun130122,Jilin,China
3Center for Disease Control and Prevention,Chengdu Military Region of PLA,Kunming650118,Yunnan,China
4The Animal Health Inspection of Hunan Province,Changsha410007,Hu’nan,China
5Jilin Key Laboratory of Animal Resource Conservation and Utilization,Northeast Normal University,Changchun130024,Jilin,China

viral metagenomic, bat, virome, next-generation sequencing

Supported by: NSFC-Yunnan Province Joint Fund (No. U1036601).

時間：2013-04-10 網(wǎng)絡出版地址：http://www.cnki.net/kcms/detail/11.1998.Q.20130410.1005.003.html

Bats are important reservoir animals and more than 60 viruses have been identified in bats with many of them highly pathogenic to human. In order to understand the natural background, genetic diversity of bat viruses in China and discover potential viral pathogens, Solexa sequencing based viral metagenomics focusing on bats tissues was established and to analyze the virome of bats collected from Jilin, Yunnan and Hunan province. By Solexa sequencing, 116 442 324 useful reads were obtained and assembled into 4 872 contigs, of which 8.2% (4 002/4 4872) were annotated to 36 viral families, including 19 vertebrate virus families, 6 plant virus families, 4 insect virus families and 4 phages. Further contigs analyses showed that some adenovirus, bocavirus, picobirnavirus, parvovirus contigs sequences were similar with known viruses. However, part of them shared limited identities to these viruses implying the discovery of new viruses. Moreover,PCR validation of adenovirus and bocavirus confirmed the results obtained by viral metagenomics. This study aimed to understand bat virome in China by viral metagenomics and could be helpful to establish effective surveillance on wildlife-associate zoonoses.

(本文責編 郝麗芳)