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全序列，经相似性、系统进化树、基因重组及抗原蛋白氨基酸关键位点分析发现，广西毒株与国内主要流行毒株亲缘性较近，抗原蛋白氨基酸序列发生独特的突变，向高致病性群体进化趋势明显，毒株存在重组现象，为进一步制定有效的防控措施提供了参考依据。