华南地区牛结节性皮肤病病毒全基因组测序及分析

周秀蓉，黄敏霞，徐志宏，温肖会，罗胜军，贾春玲，高 原，翟 颀*，吕殿红*

(1.广东省农业科学院 动物卫生研究所 广东省畜禽疫病防治研究重点实验室/农业农村部兽用药物与诊断技术广东科学观测实验站/岭南现代农业科学与技术广东省实验室肇庆分中心，广东 广州 510640；2.仲恺农业工程学院 动物科技学院，广东 广州 510225;3.珠海市启信达生物科技有限公司，广东 珠海 519000)

牛结节性皮肤病(lumpy skin disease,LSD)是一种对养牛业具有重大经济意义的动物传染病，其病原体是牛结节性皮肤病病毒(lumpy skin disease virus,LSDV)，属于痘病毒科(Poxviridae)，脊索痘病毒亚科(Chordopoxvirinae,ChPV)[1]。与山羊痘病毒(goatpox virus,GTPV)和绵羊痘病毒(sheeppox virus,SPPV)一同归属山羊痘病毒属(Capripoxvirus，CaPV)。最近LSD频繁暴发于东亚和南亚国家，包括孟加拉国、中国、印度、尼泊尔[2]、不丹、越南、缅甸、斯里兰卡、泰国[3]、马来西亚、柬埔寨和老挝等[4]。

2019年我国与印度、孟加拉国几乎同一时期暴发LSD。其中孟加拉国最先于2019年7月暴发LSD，印度首次发生LSD是在2019年8月东部地区[5-6]。我国首例LSD疫情暴发于2019年8月邻近哈萨克斯坦的新疆伊犁州，张敏敏等[7]从患病动物身上采集皮肤组织样本分离出LSDV/China/Xinjiang/2019，并得到其全基因组序列，发现与2017年俄罗斯LSDV/Russia/Saratov/2017高度同源，推测我国LSD疫情为邻国传入，但是该毒株具有在LSDV/Russia/Saratov/2017中未观察到的新进化特征，我国LSD暴发的LSDV分离株的真正来源仍需进一步调查。此后间隔了将近9个月，于2020年6月传播到华南地区，我国多个省份开始密集发生了多起LSD疫情，至今已波及了十几个省份[4,8]。期间，越南也在2020 年 11 月 1 日向 OIE 首次报告LSD疫情[9]。MA等[10]获得了2020年我国广东的1株LSDV分离株(China/GD01/2020)全基因组序列，确定其为1株田间强毒与疫苗的重组毒株，并根据GPCR和PRO30基因遗传进化树推测可能与新疆LSD疫情有共同来源。China/GD01/2020是目前鉴定的第3株LSDV重组毒株，另外2株重组毒株是俄罗斯的LSDV/Russia/Saratov/2017和LSDV/Russia/Udmurtiya/2019[11-12]。但我国2020年发生的其他几起LSD疫情来源尚未知，与2019年新疆LSD疫情来源是否一致仍然不确定，同时LSDV基因组有高度保守性和遗传稳定性，必须对整个基因组进行分析和监测，以进一步发现LSDV的遗传变化特征，掌控国内LSD疫情的传入源头和流行趋势，更加精准地防控LSDV。

本研究自2020年6月我国广东省潮州市和2020年10月广西壮族自治区百色市公布LSDV疫情的牛群中分别采集到疑似LSDV感染的牛皮肤组织样品，经qPCR检测和透射电镜鉴定后，使用高通量测序技术对2份样品进行全基因组测序、序列比较、遗传进化分析和重组事件分析研究，以期为国内疫情溯源提供线索，为国内LSDV的防控策略制定提供一定参考。

1 材料与方法

1.1 样品及处理疑似LSDV感染的牛皮肤组织样品分别采自我国广东省潮州市和广西壮族自治区百色市公布LSD疫情的牛群。取部分疑似LSDV感染的牛皮肤组织样品向其中加入适量无菌生理盐水匀浆混悬，制成约10% 组织匀浆，8 000×g离心2 min 后弃沉淀，取上清液用0.45 μm过滤器过滤，收集滤液用于后续试验。

1.2 主要试剂羊痘病毒、LSDV核酸双重荧光PCR检测试剂盒购自郑州中道生物技术有限公司；病毒基因组DNA/RNA提取试剂盒、2×Taq PCR Master mix(KT201)和50×TAE核酸电泳缓冲液购自天根生化科技(北京)有限公司；DNA Marker购自宝生物工程(大连)有限公司；核酸凝胶产物回收试剂盒购自 Axygen 生物技术有限公司。

1.3 主要仪器冷冻高速离心机购自赛默飞世尔(苏州)仪器有限公司；BIO-RAD T100 梯度PCR仪购自伯乐生命医学产品(上海)有限公司；qTOWER3荧光定量PCR仪购自德国耶拿分析仪器股份有限公司；HT7700透射电子显微镜购自日立高新技术公司。

1.4 样品检测将上述样品滤液进行灭活后提取病毒总DNA，利用羊痘病毒、LSDV核酸双重荧光PCR检测试剂盒进行检测，阳性对照和阴性对照均为试剂盒提供，同时设山羊痘疫苗为对照组。

1.5 病毒电镜鉴定将新鲜制备的组织滤液缓慢滴加到400目铜网上，25℃吸附铜网1 min，转移到1滴Tris-EDTA 缓冲液(pH 7.8) 中20 s，2% 磷钨酸染液(pH 7.2)负染10 s后，用滤纸吸掉多余的磷钨酸，自然干燥后，透射电镜下观察病毒粒子形态，观察后将电镜读片进行高压灭菌处理。

1.6 全基因组测序将1.1保存的样品滤液灭活后，重新提取其中病毒基因组DNA送广州和沁生物科技有限公司进行Illumina高通量测序。

1.7 基因组组装、注释和序列比较用在线BLAST工具进行验证并筛选组装序列(scaffolds)，参考毒株20L81_Bang-Thanh/VNM/20(MZ577-076)的全基因组序列，使用Primer Premier 5.0软件设计特异性引物以扩增不同scaffolds间隔区序列(表1)，分别以2株病毒总DNA为模板进行PCR扩增，胶回收后，送至生工生物工程(上海)股份有限公司进行测序，最后利用重叠序列拼接scaffolds，进一步组装分离株全基因组序列。使用GATU软件(https://4virology.net/virology-ca-tools/gatu/)参照China/GD01/2020(MW355944)的全基因组序列对分离株基因组进行编码基因功能注释。使用EditSeq软件、SnapGene 3.1.1软件和MegAlign软件进行分离株全基因组序列之间比较分析。

1.8 病毒全基因组进化分析使用在线BLAST工具与GenBank数据库中的全基因组序列进行比对，下载所有可用的LSDV全基因组序列，以及部分GTPV和SPPV全基因组序列，使用MAFFT在线工具(https://mafft.cbrc.jp/alignment/server/index.html)进行全基因组比对，并利用邻近法构建遗传进化树，设置bootstrap值为1 000。

表1 拼接基因序列的特异性引物

1.9 病毒重组分析参照GenBank数据库中可用的LSDV全基因组序列和我国推荐使用的GTPV-AV41疫苗株，与GD02和GX01构成比对，使用RDP 4.100软件中的RDP、Bootscan、MaxChi、GENECONV、Chimaera和 3Seq方法进行分析，识别基因组内潜在的重组事件。

2 结果

2.1 样品检测结果根据样品核酸检测结果(图1A、B)，判定为LSDV阳性，阴、阳性及对照均正常。

1.待检广东潮州核酸样品；2.待检广西百色核酸样品；3.GTPV疫苗核酸样品；4.阳性对照

2.2 病毒电镜鉴定经透射电镜观察，病毒粒子呈卵圆形或砖形，广东潮州样品病毒粒子平均大小为300 nm×283 nm，广西百色样品病毒粒子平均大小为267 nm×233 nm，2份样品的病毒粒子符合LSDV的细胞内成熟病毒(intracellular mature virus，IMV)(图2A、C)和LSDV细胞外包膜病毒(extracellular enveloped virus，EEV)(图2B、D)的外观和大小。

A.放大3万倍的广东潮州样品病毒粒子；B.放大3万倍的广东潮州样品病毒粒子;C.放大1.5万倍的广西百色样品病毒粒子；D.放大4万倍的广西百色样品病毒粒子

2.3 测序数据组装和注释对高通量测序数据进行比对筛选后，设计引物进行间隔序列的PCR扩增(图3)，测序后利用重叠序列组装拼接scaffolds，获得了2个大小分别为150 160，150 665 bp的全基因组连续序列，并将其分别正式命名为China/GD02/2020株(以下简称为GD02)和China/GX01/2020株(以下简称为GX01)。以China/GD01/2020 (MW355944)为参考基因组，利用GATU软件对其分别进行开放阅读框(open reading frame，ORF)功能注释，GD02分离株有149个ORF，GX01分离株有147个ORF，将注释后的GD02和GX01全基因组序列上传GenBank数据库后分别获得登录号为OM803091和OM803092。

M.DL2000 DNA Marker；1.GD-1引物的PCR扩增产物；2.GD-2引物的PCR扩增产物；3.GD-3引物的PCR扩增产物；4.GD-4引物的PCR扩增产物；5.阴性对照；6.GX-1引物的PCR扩增产物；7.GX-2引物的PCR扩增产物；8.GX-3引物的PCR扩增产物；9.GX-4引物的PCR扩增产物；10.GX-5引物的PCR扩增产物

2.4 基因组序列间比较经BLAST工具比对分析，发现GD02分离株和GX01分离株有99.99%的相似性，其中GX01分离株包含GD02分离株所有的基因序列，存在3处差异，主要在基因组两端(图4)。与参考基因组GD01株序列进行比较分析，GD02株和GX01株在基因组的5′端非编码区有25个核苷酸(nt)缺失，GD02额外缺失15 nt；在基因组3′端非编码区插入46 nt，GX01额外插入39 nt。在基因组中间编码区域有4个基因发生碱基突变和3处片段发生碱基缺失。其中ORF001和ORF103分别发生1个氨基酸替换，ORF071发生同义氨基酸突变，ORF155有2 nt突变，导致1个氨基酸替换，且翻译提前终止。编码Kelch 样蛋白的ORF144发生1 nt 缺失造成移码突变，1个氨基酸被替换，产生截短蛋白；另外在ORF143和ORF144之间的编码基因间隔区有1 nt缺失，ORF154和ORF155之间的编码基因间隔区有7 nt突变。GD02在中间编码区域没有发生碱基插入，GX01在1个未注释ORF(ORF22和ORF24之间)发生插入T，导致组氨酸替换5号氨基酸位点的亮氨酸，以及丝氨酸替换8号位点的缬氨酸，且使得翻译提前终止于该位置，中断了1个大小为219 bp 的ORF读取(表2)。

图4 GD02与GX01的基因组序列差异

表2 GD02分离株与GX01分离株的核苷酸序列改变及对编码序列的影响

2.5 全基因组遗传进化分析在GenBank数据库下载43株可用的LSDV全基因组序列信息，同时下载3株SPPV毒株和4株GTPV毒株，包括国内广泛使用的山羊痘弱毒疫苗GTPV_AV41毒株。由遗传进化分析(图5)发现，CaPV毒株之间同源性很高，超过97%，LSDV毒株与GTPV毒株较SPPV毒株的亲缘性更近，LSDV毒株可以大致分为三大分支，14株疫苗株或疫苗样毒株(vaccine/vaccine-like group) 大分支、19株田间分离株(field group)大分支，以及被夹在两大分支中间的12株疑似重组毒株(suspected recombinant group)，本研究的GD02分离株和GX01分离株也归属这一类(带三角形标记)，其中有3株毒株已确定为重组毒株(带圆形标记)。但是在疫苗株大分支和田间株大分支的中间，这12株疑似重组毒株没有共同形成1个大分支，存在4个小亚支，其中GD02株和GX01株与源于越南的4株LSDV分离株、GD01/2020株、中国台湾毒株Taiwan/2020和中国香港HongKong/2020毒株同源性都很高，它们共同组成同一个小分支。

带三角形红色字体为本研究的2株毒株，带圆形绿色字体的是此前研究中已确定的3株重组毒株

2.6 病毒重组分析对46个山羊痘病毒属的全基因组序列使用RDP 软件分析，得到大多数预测事件都将上述9株疑似重组毒株和3株已确定为重组毒株的基因组序列推测为重组体，共显示了73个不同的假定重组事件，且所有这些事件均具有统计学意义，通过RDP软件包中的RDP、GENECONV、Bootscan、Maxchi、Chimaera和3Seq至少一种方法计算，P值均小于0.05。图6显示了12个重组体的重组区域(浅色区域)，重组模式可以大致分为4种，重组模式Ⅰ由5株中国地区的LSDV毒株与4株越南LSDV毒株组成，重组模式Ⅱ、Ⅲ、Ⅳ分别为俄罗斯2017年LSDV分离株、哈萨克斯坦2018年LSDV分离株和俄罗斯2019年LSDV分离株。重组模式Ⅰ的主要亲本为南非疫苗株(Neethling vaccine LW 1959)，次要亲本主要为摩洛哥从牛获得的田间株(LSD/KSGP 0240/Morocco/2017)，次要亲本也有南非野毒(Neethling Warmbaths LW)、以哈萨克斯坦2016年LSDV野毒分离株改造的疫苗分离株(Neethling-RIBSP vaccine)，以及俄罗斯2015年LSDV野毒(LSDV/Russia/Dagestan/2015)。值得注意的是，它们都有1个假设重组事件是以GTPV-AV41为次要亲本。

图6 GD02和GX01等12个重组体的重组区域(浅色)

3 讨论

本研究使用商品化检测试剂盒鉴定了来自华南地区的2份疑似LSDV感染的皮肤结节组织为LSDV阳性，这2份样品分别来源于广东省潮州市和广西壮族自治区百色市公布LSD疫情的牛群。痘病毒是一类最大的动物病毒，病毒粒子呈砖形或椭圆形，LSDV大小约为290 nm×270 nm[13]。本研究将处理后的病毒样品通过透射电镜观察病毒粒子表面形态，捕捉到LSDV的2种病毒粒子——IMV和EEV，并观察到有明显的囊膜和外膜结构特征。IMV外膜表面具有清晰的纹理结构，是IMV外膜表面分布的若干不规则管状结构蛋白；EEV囊膜下的外膜横截面有纹理结构且有一定厚度[14]。本研究拍摄的LSDV电镜图丰富了LSDV病毒粒子微观形态图像，根据标尺估算图中病毒粒子大小发现GD02株病毒粒子平均大小为300 nm×283 nm，GX01株病毒粒子平均大小为267 nm×233 nm，其中GD02株病毒粒子符合一般的LSDV粒子大小，GX01株病毒粒子略微偏小，可能是病毒处于不同的生长时期，而电镜所拍到病毒粒子过少造成的统计偏差。值得注意的是，2株病毒粒子的大小都超过0.22 μm过滤器的滤过范围，属于大型病毒粒子，这提示在病毒分离培养时需要使用0.45 μm过滤器，而不是一般病毒分离所使用的0.22 μm过滤器。

本研究利用高通量测序获取了2株华南地区LSDV全基因组序列，将其分别命名为China/GD02/2020(OM803091)和China/GX01/2020(OM803092)，这为国内LSDV流行毒株的遗传进化和病毒溯源研究奠定了基础。LSDV基因组序列有高度保守性和遗传稳定性，GD02和GX01核苷酸同源性为99.99%，差异主要存在于基因组两端。与参考毒株China/GD01/2020相比，核苷酸相似性均为99.9%，在编码区存在4个氨基酸替换，造成2个截短ORF的产生。对GD02和GX01全基因组序列进行比对和遗传进化分析，发现我国2020年暴发的多起LSD疫情(广东、广西、中国台湾和中国香港)来源一致，与越南2020年LSDV分离株亲缘性最高，同源性为100%，其次依次是俄罗斯2017年LSDV分离株(Saratov/2017)以及哈萨克斯坦2018年的LSDV分离株(KZ-Kostanay-2018)，与后两者的遗传距离相等。进一步使用RDP软件分析全基因组序列，确定与GD02和GX01处于同一小分支的2020年中国香港、中国台湾以及越南的LSDV分离株都是疫苗毒和田间毒的重组毒株，并且有共同来源的亲本，其具体亲本来源无法确定，可能是多个亲本的重组，其中一个可能的亲本是我国正在使用的GTPV-AV41减毒活疫苗。俄罗斯一直使用SPPV减毒活疫苗来防控LSD，但在分析2017年和2019年LSDV重组毒株的重组事件中未见LSDV与SPPV发生重组[12,15]。目前并没有更多发生LSDV和GTPV重组的证据，这需要之后更多病毒基因组信息的监测分析。此外，哈萨克斯坦2018年的LSDV分离株(KZ-Kostanay-2018)也是重组毒株，哈萨克斯坦使用肯尼亚兽医疫苗生产研究所(KEVEVAPI)生产的lumivax®疫苗[16]，遗传进化分析显示该分离株与俄罗斯2017年和2019年的LSDV重组毒株处于不同分支。我国首起LSD疫情于2019年8月发生在新疆伊犁哈萨克自治州察布查尔锡伯自治县，距哈萨克斯坦边境约20 km，距俄罗斯边境约780 km[8]，此起疫情发生在哈萨克斯坦2017年开始使用LSDV同源疫苗的接种活动之后[17]。根据张敏敏等[7]构建的全基因组遗传进化树，确定LSDV/Xinjiang/2019分离株也是1株LSDV重组毒株，由于LSDV/China/Xinjiang/2020分离株的全基因组序列还没有在GenBank公开，因此并不能准确判断我国2020年LSD疫情是否与2019年新疆LSD疫情来源一致。

鉴于LSDV在我国的大范围流行，以及其他国家使用减毒活疫苗免疫失败的案例[18-19]，有必要重新评估减毒活疫苗的安全性和有效性，同时加快研发更加安全有效的新型疫苗。目前国外已有学者研发出一种安全有效的油性佐剂灭活LSDV疫苗，已在实验室和野外条件下测试了该疫苗对牛的安全性、免疫原性和有效性，但仍需对免疫持续时间等方面进行进一步的探索和完善，尚未商品化生产[20]。为了控制LSDV在我国的进一步传播，迫切需要在全国范围内调查LSDV的真实流行情况，并对更多的LSDV分离株进行全基因组测序，以便分析国内LSDV流行毒株的遗传演化过程，为进一步的LSD疫情溯源和防控策略制定提供依据。

本研究通过对疑似LSD的样品进行qPCR检测和透射电镜鉴定，确定了来自我国华南地区的2份样品均为LSDV阳性。通过LSDV全基因组测序、序列比对、遗传进化分析和重组分析，获得了2株华南地区LSDV的全基因组序列，并确定它们都是重组毒株，进一步确定2020年中国香港、中国台湾和越南的LSDV分离株都属于重组毒株，为国内LSDV防控策略的制定提供了参考。