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

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

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

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

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

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

1 材料與方法

1.1 主要試劑

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

1.2 病料檢測

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

1.3 DTMUV全序列擴(kuò)增與測序

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

表1 DTMUV廣西毒株信息

1.4 序列分析

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

表2 用于擴(kuò)增DTMUV全基因序列的引物

續(xù)表

表3 DTMUV參考毒株

續(xù)表

2 結(jié) 果

2.1 DTMUV全序列擴(kuò)增與測序

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

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

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

相似性比對析結(jié)果顯示，廣西毒株之間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)；與國內(nèi)外參考毒株的序列相似性分別為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蛋白氨基酸序列分析

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

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

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

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

2.4 遺傳進(jìn)化分析

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

圖2 基于DTMUV ORF核苷酸序列的系統(tǒng)進(jìn)化樹Fig.2 Phylogenetic analysis based on ORF nucleotide sequences of DTMUV

圖3 基于DTMUV E基因核苷酸序列的系統(tǒng)進(jìn)化樹Fig.3 Phylogenetic analysis based on E gene nucleotide sequences of DTMUV

2.5 DTMUV全序列重組分析

應(yīng)用RDP5和SimPlot軟件對廣西毒株和參考毒株全序列進(jìn)行重組分析，結(jié)果顯示，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基因核苷酸序列的系統(tǒng)進(jìn)化樹Fig.4 Phylogenetic analysis based on NS5 gene nucleotide sequences of DTMUV

2.6 E、NS5基因遺傳進(jìn)化速率與tMRCA估算結(jié)果

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

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基因遺傳進(jìn)化速率和tMRCA估算

3 討 論

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

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

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

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

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

4 結(jié) 論

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