2019-2021年鴨坦布蘇病毒廣西流行毒株遺傳多樣性分析

熊陳勇，尹彥文，施開創，李 軍，鄭 敏，韋顯凱，馮淑萍，龍 鳳，屈素潔，陸文俊，周洪槿，黃海蓮，謝守玉，黎宗強

(1.廣西大學動物科學技術學院，南寧 530005；2.廣西動物疫病預防控制中心，南寧 530001；3.隆安縣城廂鎮水產畜牧獸醫站，南寧 532799)

坦布蘇病毒(Tembusu virus，TMUV)于1955年首次在馬來西亞庫蚊體內被發現[1]，感染宿主譜廣泛，鴨、鵝、雞、麻雀、蚊蟲等均能被感染[2]。2010年，中國福建、浙江、江蘇等地的蛋種鴨場突發新疫情，以蛋鴨產蛋驟降、卵巢出血性壞死、甚至死亡為主要特征[3-5]，后經證實致病原為鴨坦布蘇病毒(Duck Tembusu virus，DTMUV)[6]，發病率高達90%，病死率為5%～10%[7]，目前中國多地均有發病的報道[8-9]，給中國水禽養殖業造成了巨大的經濟損失，嚴重威脅養禽業的健康發展。

DTMUV是屬于黃病毒科黃病毒屬恩塔亞病毒群的單股正鏈RNA病毒，基因組全長約11 kb，包括5′-及3′-非編碼區(UTR)和1個開放閱讀框(ORF)。ORF編碼3種結構蛋白，分別為核衣殼蛋白C、膜蛋白prM和囊膜糖蛋白E；7種非結構蛋白，分別為NS1、NS2A、NS2B、NS3、NS4A、NS4B和NS5[10-11]。E蛋白是病毒粒子表面主要結構蛋白，也是中和抗體的靶向蛋白，參與調節病毒受體結合、宿主范圍及毒力關鍵因子等重要過程[12]。NS5蛋白是黃病毒中最大的非結構蛋白，且不同病毒間最為保守，相關研究表明，NS5蛋白能抑制干擾素信號轉導作用，進而阻礙天然免疫反應[13]。因此，E和NS5基因是DTMUV分子流行病學研究、病原學診斷及研發疫苗和藥物的理想靶標。

廣西地處候鳥遷徙過冬的沿海區域，毗鄰東南亞國家的邊境口岸，同時也是中國的養鴨大省。為進一步了解近年來DTMUV廣西毒株的流行情況及遺傳進化特征，本研究對2019-2021年采集自廣西各地的病死鴨組織進行DTMUV檢測，并根據毒株檢測年份及來源地選擇部分陽性樣品開展全基因組測序和序列分析，以期了解并掌握DTMUV廣西毒株的遺傳變異新特點，為有效防控DTMUV感染提供參考數據。

1 材料與方法

1.1 主要試劑

核酸提取試劑盒購自蘇州天隆生物科技有限公司；PrimeScriptⅡ 1st Strand cDNA Synthesis Kit、2×TaqPCR MasterMix試劑盒、MiniBEST Agarose Gel DNA Extraction Kit Ver.4.0、pMD18-T載體、大腸桿菌DH5α感受態細胞均購自寶生物工程(大連)有限公司。

1.2 病料檢測

病料來自2019-2021年廣西各地養鴨場送檢的臨床病死鴨，采集肝臟、脾臟、腎臟、肺臟、心臟、氣管和腦等組織樣品，置于含有適量PBS溶液的滅菌離心管中，經組織磨碎儀研磨至糜狀，離心取上清提取總RNA，采用PrimeScriptⅡ 1st Strand cDNA Synthesis Kit反轉錄為cDNA，應用廣西動物疫病預防控制中心獸醫實驗室已建立的DTMUV實時熒光定量RT-PCR方法進行檢測[14]。

1.3 DTMUV全序列擴增與測序

廣西養鴨主產區集中于桂東、南及東南，根據毒株檢測年份及來源地，選取2019年北海、鐘山，2020年鐘山、上思、欽州，2021年蒙圩、白沙、下灣(均屬于貴港市)各1份共計8份經檢測為陽性的病料cDNA為模板(表1)，根據DTMUV參考株序列(GenBank登錄號：KC990542)，應用設計的特異性引物(表2)進行DTMUV全序列擴增。PCR反應體系50 μL：2×Gflex PCR Buffer (Mg2+,dNTP plus) 25 μL，上、下游引物(10 pmol/μL)各1 μL，TksGflex DNA Polymerase 1 μL，cDNA 5 μL，滅菌雙蒸水補足體系。PCR反應程序：94 ℃預變性1 min；98 ℃變性10 s，50 ℃退火15 s，68 ℃延伸2 min，共30個循環；68 ℃延伸10 min。PCR產物經1.2%瓊脂糖凝膠電泳檢測，目的條帶使用膠回收試劑盒回收純化，連接至pMD18-T載體中，轉化大腸桿菌DH5α感受態細胞，每個重組子挑取3個陽性菌落送測序。對測序所獲片段進行拼接，獲得DTMUV基因組全序列。

表1 DTMUV廣西毒株信息

1.4 序列分析

應用BioEdit軟件對廣西毒株及參考毒株(表3)進行序列比對和相似性分析。使用Mega 7.0軟件計算ORF及E、NS5基因最佳核苷酸替換模型為GTR+G+I、TN93+G、GTR+G+I。應用Bepipred Linear Epitope Prediction 2.0在線軟件(http:∥tools.immuneepitope.org/bcell/)預測E蛋白潛在的B細胞抗原表位；應用NetNGlyc 1.0在線軟件(https:∥services.healthtech.dtu.dk/service.php?NetNGlyc-1.0)預測E蛋白N-糖基化位點。以最佳核苷酸替換模型為基礎，采用最大似然法(maximum likelihood)繪制系統進化樹，Bootstrap值定義為1 000次。應用RDP 5(recombination detection program 5)和SimPlot v 3.5.1軟件對廣西毒株及參考毒株全序列進行重組分析。應用BEAST v 1.10.4軟件對E、NS5基因估算遺傳進化速率及最近共同祖先時間(the most recent common ancestor,tMRCA)。

表2 用于擴增DTMUV全基因序列的引物

續表

表3 DTMUV參考毒株

續表

2 結 果

2.1 DTMUV全序列擴增與測序

以選取的8份DTMUV陽性病料cDNA為模板，利用設計的特異性引物分段擴增DTMUV全序列，結果顯示，擴增片段大小與預期相符(圖1)。經測序、拼接得到8株DTMUV全序列，大小均為10 992 bp。經E蛋白氨基酸序列比對表明，8株均為DTMUV野毒株，分別命名為：GXZS01-2019、GXBH01-2019、GXSS01-2020、GXZS02-2020、GXQZ01-2020、GXMW01-2021、GXXW01-2021及GXBS01-2021，為廣西毒株。

1～8，P1～P8片段PCR擴增產物；M,DL5000 DNA Marker1-8,PCR amplification products of P1-P8 segments，respectively；M,DL5000 DNA Marker圖1 DTMUV全序列PCR擴增結果Fig.1 PCR amplification result of DTMUV genome sequence

2.2 核苷酸與氨基酸序列相似性比對

相似性比對析結果顯示，廣西毒株之間ORF及E、NS5基因核苷酸與氨基酸序列相似性分別為97.9%～99.8%和99.1%～99.9%、97.0%～99.9%和99.0%～100%、97.8%～99.9%和99.3%～99.9%(表4)；與國內外參考毒株的序列相似性分別為86.4%～99.3%和96.0%～99.8%、85.8%～99.5%和94.8%～100%、87.1%～99.2%和97.5%～99.9%(表5)。其中，與水禽源TMUV的相似性分別為89.2%～99.2%和96.0%～99.8%、89.1%～99.5%和94.8%～100%、89.7%～99.2%和97.5%～99.9%；與雞源TMUV的相似性分別為86.7%～96.5%和96.1%～99.2%、86.5%～96.5%和96.2%～98.8%、87.3%～96.8%和97.6%～99.6%；與其他宿主源TMUV的相似性分別為86.4%～96.5%和96.0%～99.4%、85.8%～96.7%和96.2%～99.4%、87.1%～96.8%和97.6%～99.7%。表明DTMUV廣西毒株間相似性較高，與水禽源TMUV相似性高于其他宿主源。

2.3 E蛋白氨基酸序列分析

應用Bepipred Linear Epitope Prediction 2.0在線軟件預測廣西毒株的E蛋白B細胞抗原表位，結果共發現16個潛在抗原表位，其中第157位的丙氨酸抗原相關指數最高。應用NetNGlyc 1.0在線軟件預測N-糖基化位點，結果發現1個糖基化位點，為154位氨基酸NYSA序列。E蛋白氨基酸序列比對顯示，與疫苗株FX2010相比，廣西毒株共有11處位點發生變異，第38、122及277位氨基酸處為不同毒株普遍發生變異；第43、150、153、326、403、464及487位氨基酸處為部分DTMUV廣西毒株特有的突變，而GXBH01-2019株在第280位氨基酸處出現缺失(表6)。表明廣西毒株E蛋白既有共性又存在獨特的氨基酸突變，但均未涉及抗原表位與糖基化位點。

表4 廣西毒株之間ORF及E、NS5基因序列相似性比對

表5 廣西毒株與參考毒株ORF及E、NS5基因序列相似性比對

表6 廣西毒株與疫苗株E蛋白氨基酸序列比對結果

2.4 遺傳進化分析

基于ORF及E、NS5基因核苷酸序列對廣西毒株與國內外參考毒株構建系統進化樹，結果顯示，ORF及E、NS5基因均可分為1、2、3及TMUV群(圖2～4)。分布1群的毒株為來自泰國的DK/TH/CU-DTMUV2007株；分布2群的毒株來自中國、泰國，其中2群又劃分為2.1和2.2亞群，廣西毒株(圖2～4中標注▲的毒株)均集中于2.1亞群，廣西的GX2013E株也屬于該群，而廣西的GX2011、GX2015株則屬于2.2亞群，表明DTMUV廣西流行毒株進化趨勢不完全一致；分布3群的為山東SD14株；蚊蟲源及雞源TMUV毒株組成TMUV群。廣西毒株GXZS01-2019、GXZS02-2020、GXBH01-2019與CHN-JL株遺傳關系較近，而其他廣西毒株與AQ-19、H株遺傳關系較近，與相似性比對結果一致。基于E、NS5基因的系統進化樹與ORF相似。

圖2 基于DTMUV ORF核苷酸序列的系統進化樹Fig.2 Phylogenetic analysis based on ORF nucleotide sequences of DTMUV

圖3 基于DTMUV E基因核苷酸序列的系統進化樹Fig.3 Phylogenetic analysis based on E gene nucleotide sequences of DTMUV

2.5 DTMUV全序列重組分析

應用RDP5和SimPlot軟件對廣西毒株和參考毒株全序列進行重組分析，結果顯示，ziYY150901、HB2016、GXBH01-2019及GXZS02-2020株均有重組信號(圖5)，其中，ziYY150901株重組的主要親本為ZJ201506株，相似性為99.7%，次要親本為SDHS株，相似性為99.5%；HB2016株重組的主要親本為ZJ201506株，相似性為99.5%，次要親本為GX2015株，相似性為100%；GXBH01-2019株重組的主要親本為GXZS01-2019株，相似性為99.8%，次要親本為GXXW01-2021株，相似性為99.9%；GXZS02-2020株重組的主要親本為GXZS01-2019株，相似性為99.8%，次要親本為GXQZ01-2020株，相似性為99.8%。

圖4 基于DTMUV NS5基因核苷酸序列的系統進化樹Fig.4 Phylogenetic analysis based on NS5 gene nucleotide sequences of DTMUV

2.6 E、NS5基因遺傳進化速率與tMRCA估算結果

應用BEAST v 1.10.4軟件，采用貝葉斯聚結分析方法，對廣西毒株與國內外參考毒株的E、NS5基因估算遺傳進化速率與tMRCA，結果顯示，E、NS5基因的遺傳進化速率分別為1.31×10-3和1.30×10-3替換/(位點·年)；tMRCA分別為88.90和88.91年(表7)。表明E、NS5基因的進化速率相近，二者進化比較同步。

A～D，ziYY150901、HB2016、GXBH01-2019和GXZS02-2020株；1，主要親本；2，次要親本A-D，ziYY150901,HB2016，GXBH01-2019 and GXZS02-2020 strains，respectively；1，Major parent;2，Minor parent圖5 廣西毒株與參考毒株全序列重組分析Fig.5 Recombination analysis of complete genome of DTMUV between Guangxi and reference strains

表7 DTMUV E、NS5基因遺傳進化速率和tMRCA估算

3 討 論

廣西地處亞熱帶氣候，與東南亞國家接壤，同時也是候鳥遷徙過冬地區。TMUV宿主譜廣泛，除感染鴨、鵝、雞等家禽外，也能感染麻雀、鴿子等鳥類及蚊蟲[15-16]，而鳥類、蚊蟲作為TMUV傳播媒介，存在散播病毒及與沿海地區養殖的水禽交叉感染的風險。另外，有報道稱，養鴨場從業人員的DTMUV血清陽性檢出率超過70%，可能引發公共衛生安全問題[17-18]。本研究采集2019-2021年廣西各地病死鴨的肝臟、脾臟、腎臟、卵巢等組織病料，利用本實驗室已建立的實時熒光定量RT-PCR方法進行檢測顯示，DTMUV與DuCV陽性檢出率較高。推測可能的原因有兩方面：一是由于DuCV等病原體引起機體免疫抑制，從而導致DTMUV疫苗接種效果不佳，達不到預期的保護作用；二是田間DTMUV流行毒株可能發生變異，疫苗株交叉保護力下降。因此，為進一步了解近年來DTMUV廣西毒株遺傳分子特征，本研究根據毒株檢測年限及來源地，從廣西養鴨主產區送檢病料中選取8份陽性提取總RNA，利用設計的特異性引物擴增DTMUV全序列，并進行遺傳多樣性分析。

相似性比對結果顯示，與8株廣西毒株ORF相似性最高的分別為H株(5/8)、CHN-JL株(2/8)及AQ-19株(1/8)，在遺傳距離上也得出同樣的結論，其中，H株為Feng等[19]2019年從3周齡的金定鴨體內分離到的隸屬2.1亞群強毒株。廣西毒株之間ORF核苷酸序列相似性最高達99.8%，與參考毒株的核苷酸序列相似性最高達99.3%，其中，與廣西GX2011、GX2013E及GX2015株的ORF核苷酸序列相似性分別為95.9%～96.2%、96.8%～97.0%及95.6%～95.9%，而E、NS5基因的相似性和ORF相當，表明廣西毒株之間相似性較高。農海連[20]分離的4株廣西DTMUV毒株與參考毒株的核苷酸序列相似性高于96%，氨基酸序列相似性高于98%；李剛[21]分離的山東DTMUV毒株與2010-2012年分離株的ORF氨基酸序列相似性為96.8%～98.6%，與2017年分離株的氨基酸序列相似性為98%以上；劉烈發等[22]從江西省某蛋鴨養殖場分離到1株DTMUV，與參考毒株核苷酸序列相似性在98.5%以上。表明國內毒株之間相似性普遍較高，遺傳變異程度小，與本研究結果相似。

經在線軟件預測，E蛋白共有16個B細胞抗原表位及1個N-糖基化位點，為E蛋白的關鍵位點。E蛋白是DTMUV重要的抗原結構蛋白，可誘導宿主產生免疫應答，使機體產生病毒中和抗體，在宿主的免疫壓力下，E蛋白將會發生遺傳變異[23]。據報道，E蛋白第156位氨基酸突變會破壞第154位氨基酸處的糖基化位點，進而影響DTMUV的復制和傳播[24]。E蛋白氨基酸序列分析顯示，廣西毒株與參考毒株在第156位氨基酸處均未出現變異。與疫苗株FX2010相比，E蛋白出現了多個位點的突變，如廣西毒株GXZS01-2019、GXZS02-2020及GXBH01-2019與參考毒株CHN-JL在第43位氨基酸處均出現相同突變，而在E基因系統進化樹中，該3個毒株與CHN-JL株遺傳距離也最近；第157位氨基酸處發生突變的毒株為雞源、蚊蟲源TMUV，以及系統進化樹中分群屬于3群的SD14毒株，屬于2群的廣西毒株與其他水禽源TMUV在該位點均未突變，而第157位氨基酸處的抗原指數最高，可能與病毒的宿主范圍或分群有關。此外，有7個氨基酸位點是部分廣西毒株特有的突變。廣西毒株的E蛋白突變位點均未涉及關鍵位點，至于這些位點突變的影響尚需進一步研究。

基于ORF核苷酸序列構建的系統進化樹顯示，廣西毒株分布于2.1亞群，同樣源自廣西的GX2013E株也屬于該群，而廣西的GX2011、GX2015株則分布于2.2亞群。Zhu等[25]研究報道，2.1亞群是中國新出現的一種強毒力型毒株，對雛鴨的致病性高于2.2亞群。表明DTMUV廣西流行毒株呈現不同的進化趨勢，更多的毒株向高致病性的群體進化。另外，廣西毒株與參考毒株的遺傳距離也不盡相同，GXZS01、GXZS02、GXBH01與CHN-JL株遺傳距離較近，而另外5個毒株與H、AQ-19株距離較近。廣西的GX2013E與HD1-2013株遺傳距離較近，而廣西的GX2011、GX2015與GD06株遺傳距離較近。E、NS5基因的系統進化樹分布與ORF相似。表明廣西毒株進化趨勢出現分歧，呈現遺傳多樣性的特征。

DTMUV全序列重組分析顯示，共有4株有重組信號，其中，ziYY150901、HB2016株重組位點相似，位于E-NS1-NS2A區域，GXBH01-2019株位于NS1-NS2A-NA2B區域，GXZS02-2020株位于E-NS1區域。表明重組位點位于E蛋白及其臨近的非結構蛋白。E蛋白作為宿主的重要保護性抗原，對病毒的吸附、融合及受體結合等發揮著至關重要的作用[24,26]。非結構蛋白在病毒核酸復制、蛋白合成加工與病毒粒子組裝等方面起關鍵性作用[27]。因此，這些基因發生重組需引起密切關注。貝葉斯遺傳進化速率分析顯示，1955-2021年E、NS5基因進化速率相近。Ninvilai等[28]估算E基因的進化速率為1.507×10-3替換/(位點·年)；Yu等[29]研究報道，E基因的進化速率為5×10-4替換/(位點·年)。本研究中，E基因的進化速率介于二者之間，可能與參考毒株的種類和數量及分析軟件參數設置等因素有關。E基因最大可信進化分支樹顯示，TMUV約出現在1 934.9年，并且經數次進化形成不同群體。通過軟件計算得出，廣西毒株早在2 016.5年就已存在(未發表)。因此，需持續密切關注DTMUV流行病學調查與分子遺傳特征分析，及時了解掌握DTMUV在廣西地區流行新態勢與新特點，為有效防控DTMU感染制定措施提供數據支持。

4 結 論

本研究擴增8株DTMUV全序列，經相似性、系統進化樹、基因重組及抗原蛋白氨基酸關鍵位點分析發現，廣西毒株與國內主要流行毒株親緣性較近，抗原蛋白氨基酸序列發生獨特的突變，向高致病性群體進化趨勢明顯，毒株存在重組現象，為進一步制定有效的防控措施提供了參考依據。