柯萨奇病毒A 组6 型中国株的分子进化和时空传播分析

贾永涛,唐佳婕,王 惠,董长征*

（1.宁波大学 医学部,浙江 宁波 315211;2.浙江省病理生理学技术研究重点实验室,浙江 宁波 315211）

人肠道病毒在物种分类上属于小RNA 病毒科(Picornaviridae)肠道病毒属(Enterovirus),分为A、B、C 和D 四个种型.柯萨奇病毒A 组6 型(Coxsackievirus A6,CV-A6)属于肠道病毒A 种,是手足口病的主要病原体之一.与同属于肠道病毒A种的肠道病毒A 种71 型(Enterovirus A71,EV-A71)和CV-A16 所引起的手足口病不同的是,CV-A6 感染常表现为非典型症状,包括全身性大疱疹、脱甲和睾丸炎等[1-4].CV-A6 和EV-A71 一样,都可引起无菌性脑膜炎和脑干脑炎等中枢神经系统疾病[5-6],造成严重的疾病负担.CV-A6 于1949 年在美国首次被分离,2008 年在芬兰引发大规模疫情[7],此后在欧洲[8-9]、美洲[10-11]和亚洲[12-13]等全球各地引发了广泛的疫情,并逐渐取代EV-A71 成为手足口病流行的优势病原体.自2013 年以来,我国CV-A6感染的病例数呈明显的上升趋势[14],CV-A6 很快成为了上海、深圳和北京等地区手足口病的优势病原体[6,15-16].近几年,广西[17]、福建[18]、江苏[19]、四川[20]和山东[21]等全国大多数省区也都出现了以CV-A6 感染为主的手足口病疫情,CV-A6 已经成为我国手足口病流行的优势病原体.虽然手足口病疫情已经在全国高发十多年,但是更多的研究关注时空聚集性,即哪些地区是手足口病的高风险地区,以CV-A6 为代表的肠道病毒的起源和传播路线尚不清楚.

CV-A6 是一种具有单股正链RNA 的无包膜病毒,病毒的基因组被衣壳包裹.病毒衣壳由四种结构蛋白(VP1、VP2、VP3 和VP4)组成[22],表面分布着能与抗体和受体结合的抗原表位.病毒与受体结合后引起构象改变,VP1 N-端从衣壳内部暴露到衣壳表面,引发衣壳裂解,然后病毒基因组从衣壳中释放到宿主细胞,完成感染宿主细胞过程.EVA71 疫苗已于2016 年由我国研发成功并上市[23],为防控EV-A71 引起的手足口病尤其是重症手足口病做出了巨大贡献,而包括CV-A6 在内的其他肠道病毒均无上市疫苗,也没有相应的抗病毒药物.因此,持续对CV-A6 等肠道病毒的变异进行监测和传播路线进行追踪对防控大规模手足口病疫情非常重要.

本实验室前期发展了肠道病毒抗原表位的生物信息学预测算法,并成功地应用到肠道病毒A种(EV-A71、CV-A16[24]、CV-A6[25]和CV-A10[26])、B 种(CV-B1 等[27]和E18[28])和D 种(EV-D68[29])等代表性肠道病毒的表位预测和研究中,发现了肠道病毒抗原表位分布的基本规律.本文以前期确定的CV-A6 抗原表位[25]为基础,分析其在CV-A6 分子进化中的作用,为病毒基因组变异监测提供参考依据.VP1 是抗原表位分布最多、突变频率最高、公共数据库中序列最丰富的结构蛋白[24],而全基因组序列则较少且病毒株采样地区分布严重不均衡.因此,本文基于CV-A6 中国株VP1 序列,利用Nextstrain平台的augur 病原体生物信息分析包,构建了CV-A6 的时间尺度分子进化树,确定了进化分支(Clade)和分支突变,分析了分支突变与抗原表位之间的联系,并进一步分析了CV-A6 中国株各分支在中国七大地理区域的时空分布和优势分支,绘制了主要分支的时空传播路线.

1 材料与方法

1.1 CV-A6 中国株VP1 序列收集

从 ViPR 病毒数据库[30](https://www.viprbrc.org/brc)和NCBI Nucleotide 数据库[31](https://www.ncbi.nlm.nih.gov/nuccore)下载所有长度大于900 bp的CV-A6 中国株VP1 序列,截止时间为2023 年6月30 日,并收集每条序列的信息,包括GenBank序列号、病毒株名、病毒株分离日期和病毒株分离省市等.剔除实验室适应性毒株、克隆株和信息不全序列,最终获得3 288条CV-A6中国株VP1序列,病毒株分离日期介于1992 年12 月和2022 年11 月之间,来自31 个省、自治区、直辖市和特别行政区.为了清晰地表示病毒株的时空分布和避免采样偏倚,以七大地理区域(华东、华南、西南、华北、华中、东北和西北)进行空间分类.分子进化树构建以CV-A6 原株(Gdula 株[32],GenBank 序列号:AF081297)的VP1 序列作为参考株.

1.2 基于Nextstrain 平台的CV-A6 中国株分子进化和时空传播分析

Nextstrain[33]是一款免费、开源的生物信息分析平台,可针对病原体基因或基因组数据,通过强大的生物信息分析工具和数据可视化功能,监测和分析病原体的分子进化和时空传播路线.本研究基于本地化部署的Nextstrain 平台,利用augur病原体生物信息分析包对CV-A6 中国株进行分子进化和时空传播分析.利用Snakemake流程管理工具[34],逐一执行augur 分析包内的分析命令.align命令调用内嵌的MAFFT 软件[35]对CV-A6 中国株VP1 序列进行多重序列比对,tree 命令调用IQ-tree工具[36]使用最大似然法构建分子进化树(Gdula 作为参考株),refine 命令调用TreeTime 工具[37]构建基于时间尺度的分子进化树并优化进化树和确定时空传播路线,clade 命令标注进化分支和确定分支突变,表位突变标注依据之前的研究结果[25],frequencies 命令构建进化分支和突变的频率分布图.最后利用Nextstrain 平台的auspice 工具对分析结果进行基于JavaScript 的网页互动可视化展示.

2 结果

2.1 CV-A6 中国株的分子进化分析

基于3 288 条CV-A6 中国株VP1 序列构建的时间尺度分子进化树如图1 所示.CV-A6 进化树分为A、B、C 和D 四个进化分支.A 分支为美国分离的原株Gdula,作为参考株;中国株包含B、C 和D 三个分支.B 分支有6 株,其中1 株于1992 年从山东分离,其余5 株于2004—2007 年从广东分离(表1).C 分支仅有1 株,于1996 年从山东分离.D分支最多,共有3 281 株;绝大多数D 分支病毒株属于D2 和D3 两个子分支,而D3 又可细分为D3a和D3b.大约在1944 年,CV-A6 进化形成A和B～D两个分支,B～D 分支是B、C 和D 分支的共同祖先(图1).B～D 分支先后在1968 和1979 年分歧成B、C 和D 三个分支.D2 和D3 子分支于1997 年形成,而D3 进一步在2004 年分歧成D3a 和D3b,其中D3a 是近十来年我国 CV-A6 流行的优势株(3 040/3 288),占病毒株总数的92.5%.随着突变的累积,新的进化分支不断形成,并替换旧的分支成为优势分支.

表1 CV-A6 中国株的时空分布

图1 基于CV-A6 中国株VP1 序列的时间尺度分子进化树

每个进化分支上都存在分支突变.对比本实验室之前的表位预测结果,可以发现CV-A6 大多数分支突变都位于表位上.B～D 分支突变包括S5T、I8V、N10S、T14A、S32T、V242I、S279T和A283T,其中S5T、I8V、N10S、T14A 和S32T位于N-端“呼吸表位”区域[25,27],V242I 位于HI环表位,S279T 和A283T 位于C-端表位.B 分支突变为GH 环表位的K210N.C 分支突变包括位于“呼吸表位”的P3S、T5A、S10G、A14T、T18A、A29T 和V71I,BC 环表位的S97N,C-端表位的T281S 和T283S.D 分支突变包括BC 环表位的Q98L,EF 环表位的G160S 和HI 环表位的I242V.D2(A54S 和V242I)、D3(T5A)、D3a(A5T 和A30V)和D3b(N27S)都位于“呼吸表位”或HI表位(V242I).可以发现,每个进化分支都有位于表位的独特的分支突变.

2.2 CV-A6 中国株的时空分布

3 288 株CV-A6 中国株中,华东地区分离的病毒株数最多,有1134 株,占34.5%;华南和西南次之,各有743 和684 株,分别占22.6%和20.8%.华东、华南和西南共分离了约80%病毒株,是全国主要的分离地区.华北、东北和华中地区各有380、161 和130 株,分别占11.6%、4.9%和3.9%.西北地区最少,只有56 株(1.7%).病毒株分离年份介于1992—2022 年这31 年间,其中1992—2001 年这10 年间只有2 株,2002—2011 年这10 年间有296株,2012—2022 年这11 年间有2 990 株,大多数病毒株都是近十年分离的.华东地区不仅病毒株数最多,进化分支类型也最多,每一种分支都在华东地区获得分离;华南地区B 和D3b 分支分离数最多;其他地区主要分离D2、D3a 和D3b 这几种常见分支,尤其是D3a 这种优势分支.

通过频率分布图,可以简便地观察CV-A6 各分支的时空分布、优势分支以及在全国七大地理区域的分布(图2(a)).2005 年前,B 和C 分支是主要流行分支,病毒株从华东和华南分离.2005 年后,D分支取代了B 和C 分支成为流行分支.D 分支除了主要子分支D2 和D3 外,还有少量其他分支(仍用D 分支命名)病毒株,主要从华东分离.D2、D3a 和D3b 是我国主要分离的CV-A6 分支,D2 出现最早,短暂地成为全国多个地区的优势分支后就被D3a替代,直到2015 年后未再分离到.D3a 约在2008年先后出现在华南、华东和华中地区,很快就成为这几个地区的优势分支,然后成为了全国的优势分支直至目前.其间出现了D3b,在全国低水平地流行.

图2 CV-A6 中国株各分支及VP1 第242 个氨基酸残基突变的时空分布

同样可以采用频率分布图监测变异的时空分布(图2(b)).以VP1 第242 个氨基酸残基为例,B 和C 分支的残基为I;D 分支发生了突变I242V,但D2分支又发生了回复突变V242I,因此D2 分支与B和C 分支一样残基都为I,其他分支(包括D、D3a和D3b)的优势残基均为V(D3a 其中一簇发生回复突变V242I).从全国范围看,1990—2005 年期间流行的是残基I,2005 年后优势残基从I 逐渐变成V.从地区看,残基V 最先在华东成为优势残基,然后是华中和华南,最后是全国其他地区,这是因为D3a 分支最先在这几个地区成为优势分支,而华东在D3a 分支流行之前,还流行了一段时间D 分支.华南虽然比华中更早出现残基V,但因为此时D2分支是华南的优势株,所以残基I 仍是优势残基,直至D2 分支被D3a 代替,残基V 才在华南成为优势残基.

2.3 CV-A6 中国株的时空传播路线

由于早期的一些进化分支(如B、C 和D 分支)病毒株数都非常少,因此本研究重点关注D3a、D2和D3b 这三种主要流行分支的时空传播路线(图3和图4).图3 展示了D3a 分支从华南和华东输出到全国的时空传播路线.2008 年,华南开始出现D3a病毒株,2010 年华东也开始出现D3a 病毒株,这一时期传播路线主要表现为华南输出到华东.2013 年,华南和华东形成相互输出的模式,华南往华东方向输出的病毒株数更多;华东开始往西南、华北、华中和西北输出病毒株.2016 年,华南和华东相互输出、华东输出全国的传播模式得到完全建立.2022 年,西南成为第三大病毒株分离地区,全国建立相互传播网络.

图3 CV-A6 中国株D3a 分支的时空传播路线

图4 CV-A6 中国株D2 分支和D3b 分支的传播路线

虽然病毒株数较少,但D2 分支同样显示了华南和华东相互输出、华东输出全国尤其是华北和华中的传播模式(图4(a)).D3b 分支由于病毒株集中在华南、华东和华北三地且病毒株数较少,全国性的传播模式还不显著,主要表现为华南向华东输出、华南和华东向全国各地区少量输出的传播特征(图4(b)).频率分布图上各地区CV-A6 分支的出现、流行和成为优势分支的过程也间接地反映了这种时空传播路线(图2).

3 讨论

自2008 年安徽阜阳大规模暴发以EV-A71 为主的手足口病疫情以来,手足口病发病率和死亡率长期位居丙类法定传染病之首,对儿童造成严重的疾病负担.但随着感染造成的群体免疫和EV-A71 疫苗于2016 年上市,手足口病的病原谱发生大幅转变: 病原体从EV-A71 和CV-A16 为主转变成CV-A6 和CV-A16 为主,重症病例的主要病原体也从EV-A71 转变成CV-A6[23,38].因此,持续对CV-A6 和CV-A16 等肠道病毒的变异进行监测和传播路线进行追踪,加强疫苗和抗病毒药物的研发成为防控手足口病等疫情的关键措施.由于全基因组测序尚未在病毒监测中得到常规应用,ViPR和NCBI 等公共数据库中CV-A6 中国株基因组序列较少,其地区分布严重不均衡,会严重干扰传播路线分析,因此本研究采用序列比较丰富、分布比较均衡的VP1 进行变异监测和传播路线分析.

CV-A6 中国株分子进化树主要包括B、C、D三个分支以及D2、D3a 和D3b 三个子分支,进化树的拓扑结构与Song 等[14]对全球病毒株分析的结果一致,只是中国未分离到A 分支和D1 等少数子分支.对分支突变的分析发现,绝大多数突变都位于抗原表位或者VP1 N-端“呼吸表位”上,这提示CV-A6 通过表位突变形成了新的进化分支,改变了病毒抗原性,从而突破了宿主的免疫屏障,促进了病毒分支的流行.监测表位上的变异及其流行趋势,对于传染病监测和预警具有较大价值.有意思的是,CV-A6“呼吸表位”集中了最多的分支突变,这与对肠道病毒B 种的研究结果[27-28]相一致,提示肠道病毒通过“呼吸表位”的突变影响病毒的脱衣壳和感染宿主细胞过程,同时逃避抗体的中和作用,增强病毒的适应性.最典型的是D3a 迅速压倒D2 和D3b 成为优势株,与“呼吸表位”突变的作用密不可分(图2).另外,VP1 有几处氨基酸残基发生多次突变,可能这些突变对病毒抗原性影响较大,值得进一步深入研究.例如,VP1第5个残基,B～D 分支发生S5T 突变,C 和D3 分支发生T5A突变,D3a 发生A5T 突变.此类位点还包括第10、14、242 和283 等残基.

CV-A6 中国株第一株是1992 年采样的.2008年安徽阜阳大规模手足口病疫情暴发后,全国各地区才开始普遍采样病毒株.直到近5 年CV-A6才成为肠道病毒的常规监测血清型.因此,CV-A6中国株分离时间集中在近十年.D2、D3a 和D3b 分支是全国各地区主要流行的进化分支,但D2 分支只在全国部分地区短时间内成为优势分支,D3b 分支则一直只占少数,只有D3a 分支长期成为优势分支,这说明D3a 分支在进化上具有独特优势,是否是VP1 表位上的残基差异造成了病毒适应性的差异值得实验研究验证.另外,除了VP1,衣壳上的其他结构蛋白如VP2 和VP3 上的表位变异或者基因组上的其他重要突变也可能是潜在原因.对病毒变异适应性影响的研究有助于对病毒变异的风险进行评估,从而帮助预警传染病疫情.

由于全球对CV-A6 早期的进化分支(如A、B和C 等分支)缺乏监测,许多分支的病毒株数都非常少,病毒株采样国家也严重不均衡(中国采样数最多,超过总数的80.0%),难以确定各个分支在全球范围内的起源和传播关系,因此本研究重点关注CV-A6 尤其是病毒株数非常充足的D3a 分支(占病毒株总数的92.5%)在国内七大地理区域之间的时空传播关系,帮助各地区监测CV-A6 及其他肠道病毒的流行情况.传播路线的确定基于亲缘地理图(phylogeographic map),基本原理是根据时间尺度的分子进化树上分支之间的亲缘关系推断不同地区病毒株之间的亲缘关系,它在新冠病毒、流感病毒和西尼罗河病毒的进化和传播研究中已经得到应用[39-41].虽然采样偏倚可能会对亲缘关系和传播路线的确定造成干扰,但整体上不会对时空间传播路线分析造成严重影响,尤其是对D3a分支.图3 很清晰地反映了华南和华东是D3a 分支在中国的起源地.2008—2010 年间,华南和华东开始出现D3a 病毒株并开始相互传播,全国其他地区直到2012 年才出现病毒株.这是因为温暖湿润、日照数较少的气候有利于肠道病毒繁殖,促进手足口病疫情的发生,因此我国东部和南部地区一直是手足口病高发地区,而西部和北部等干燥少雨地区相对来说疫情较轻[42].特别值得注意的是,只有华东是全国各个地区病毒株的主要输入地,而华南主要输出到华东和西南,具体原因尚不清楚,需要后续深入研究.可能一方面华东与全国人流往来比较密切,另一方面华东的病毒株具有更高的适应性且不断地进化,因此在人群传播中具有优势,更容易成为优势株.传统上西南并不是手足口病疫情的高发地区,但近十年,西南的手足口病报告发病率和病死率都迅速上升[43],病毒株分离数也排在全国第三位,提示要加强对西南疫情的监测和防控,提高医疗救治水平.华南和华东相互输出、华东输出全国这种传播模式不仅在D3a中表现得非常典型,在D2 和D3b 中也得到一定程度的反映.这种传播模式是否反映了肠道病毒甚至其他一些病原体的传播模式值得深入研究.EV-A71 和CV-A6 是不同血清型的肠道病毒,研究证实EV-A71 疫苗和CV-A6 不会出现交叉免疫[23],不会对CV-A6 的时空分布和传播路线造成直接影响.但EV-A71 和CV-A6 等血清型在肠道病毒生态中竞争成为优势血清型,而EV-A71 疫苗的广泛接种会抑制EV-A71 的流行,因此会对CV-A6 的时空分布和传播路线造成间接影响,今后的研究将探索这种间接影响.

综上,本研究发现CV-A6 通过抗原表位上的突变产生了新的进化分支,促进病毒的流行;华南和华东是CV-A6 中国株的主要起源地,华南和华东相互输出、华东输出全国是CV-A6 中国株的主要传播模式,这为CV-A6 引起的手足口病等疫情的监测、预警和防控提供了重要信息.