相关疫苗上市后中国大陆肠道病毒71型C4亚型的分子流行特征分析

杜瑞晓，王茜，刘明琛，高帆，白玉，吴星，毛群颖，梁争论，卞莲莲

中国食品药品检定研究院国家卫生健康委员会生物技术产品检定方法及其标准化重点实验室，北京 102629

手足口病（hand，foot and mouth disease，HFMD）多发于5岁以下儿童，通常引起以手、足、口咽黏膜和臀部丘疹样病变为特征的轻症症状，部分病例会引发与神经系统疾病相关的重症，如无菌性脑膜炎、脑炎和脊髓灰质炎样麻痹［1-5］。HFMD在亚太地区持续广泛流行，新加坡、韩国、马来西亚、日本、越南、中国大陆和中国台湾均发生过大规模暴发［4，6-13］，引发全球关注。肠道病毒71型（enterovirus type 71，EV71）属于小核糖核酸病毒科肠道病毒属，是引起HFMD重症和死亡病例的主要病原体之一。自1981年上海发现首例HFMD后，国内大部分省份陆续报道了HFMD病例［14］。2007年后，在中国大陆中部和东部陆续发生EV71的大范围流行，如2008年，安徽阜阳暴发HFMD疫情，造成了大量儿童感染及多个重症病例死亡，因此，2008年5月，将HFMD列为我国法定报告的丙类传染病，国内企业启动单价EV71全病毒灭活疫苗的研发并得到了迅速发展［15-18］。2015—2016年，北京科兴生物制品有限公司、中国医学科学院医学生物学研究所、武汉生物制品研究所有限责任公司生产的EV71疫苗相继获得国内上市许可。3家企业的EV71疫苗在Ⅲ和Ⅳ期临床试验中均表现出良好的保护效果、一致性和有效性［16］。

目前，已报道的EV71流行株主要为B和C基因型，进一步分为B0～B5和C1～C5亚型［19］。C4亚型为2000年后在国内流行的绝对优势基因亚型［11］，通常包括C4a和C4b亚型分支，已上市的EV71疫苗均以C4a亚型毒株为疫苗株，已有数百万婴幼儿接种了EV71疫苗，有效防控了国内EV71相关的HFMD流行，但目前对EV71疫苗上市后C4亚型进化规律鲜有报道。因此，本研究对EV71疫苗上市前后，中国大陆EV71流行株分子流行病学规律进行分析，旨在监测EV71的遗传进化特征，以期为制定HFMD预防和控制策略提供实验依据。

1 材料与方法

1.1 数据筛选 从国家生物技术信息中心（NCBI，http：／／www.ncbi.nlm.nih.gov／）的GenBank数据库中下载所有已知采样日期和地理位置（国内各省）的EV71毒株VP1基因序列（截止2021年1月1日）。应用在线工具CD-HIT（http：／／weizhong-lab.ucsd.edu／cdhit_suite／cgi-bin／index.cgi）删除冗余序列，相似性阈值设置为98%。应用MAFFT v7.222软件比对序列［20］，BioEdit v7.2.5软件进行编辑，截取VP1基因区段。应用RDP v4.36软件筛选多序列比对的重组病毒序列［21］。以已报道的不同基因型EV71毒株为参考序列，应用iqtree 1.6.2软件构建最大似然（maximum likelihood，ML）系统发育树，核苷酸取代模型设置为SYM+R5，设bootstrap值为1 000以评估ML系统发育树拓扑结构的可靠性，并应用FigTree v1.4.4软件进化ML系统发育树的可视化，以获得国内分离的EV71 C4亚型VP1基因序列作为最终数据集。将选定的C4亚型毒株数据集和参考序列原始毒株BrCr（GenBank登录号为：U22521）合并，应用iqtree 1.6.2和MrBayes 3.2.6软件分别构建ML系统发育树和贝叶斯系统发育树，进一步对国内的C4亚型毒株进行分类，并应用Mega 7.0.20软件计算C4a和C4b亚型毒株VP1基因的进化距离。

1.2 系统发生地理学分析 利用C4亚型毒株数据集研究EV71在中国大陆的系统发生地理学史。根据华中、华东、华北、东北、西北、华南和西南7个区域对数据进行标记。应用Tracer v1.6软件模型比较功能选择最合适的模型，确定BEAST v1.10.5pre中贝叶斯天际线模型（Bayesian skyline）及未校正的对数正态松弛分子钟模型，该贝叶斯天际线模型同时可用于评估相对遗传多样性随时间的变化［22］。使用贝叶斯马尔可夫链蒙特卡洛（Bayesian Markov chain monte carlo，MCMC）框架时，设置链长2亿，采样频率10 000，并去除10%的老化样本。设GTR+G（由JModelTest确定为最佳拟合模型）为核苷酸替换模型，对序列进行地理信息标记，并用非对称取代模型设置贝叶斯随机搜索变量选择（Bayesian stochastic search variable selection，BSSVS）分析方法。剩余样本应用Tracer v1.7.1软件检查所有估计的参数是否有足够的ESS值（＞200）以达到收敛。应用TreeAnnotator软件生成最大分支可信度（maximum clade credibility，MCC）树，并使用Figtree v1.4.4软件进行图形化。另外，应用SpreaD3 v0.9.6软件将系统地理学数据可视化，并计算不同地区之间的迁移路径、后验概率及贝叶斯因子（bayes factor，BF），一般认为BF≥10为显著迁移路径。

1.3 群体动态分析 为确定EV71 C4亚型在国内的群体动态分布，分析EV71 C4亚型VP1基因遗传多样性随时间的变化，对包含确切采样日期（年份）的数据进行分析并绘制贝叶斯天际线图（Bayesian skyline plot，BSP）。BSP图谱可描绘有效种群规模（effective population size，Ne）随时间的变化，揭示相对遗传多样性的变化规律。

1.4 选择压力位点分析 应用EasyCodeML软件中位点特异性模型分析不同时期（2007年前、2008—2010年、2011—2014年和2015—2018年）EV71 C4亚型VP1基因的选择压力［23］。通过绘制BSP重建EV71的种群历史，采用M7-M8模型分析EV71 C4亚型VP1基因遗传多样性随时间的动态变化。

2 结 果

2.1 数据筛选及数据集建立 截止2021年1月1日，GenBank数据库中共上传了7 995个包含地理和时间信息的EV71序列，其中4 125个序列来自中国大陆，经筛选最终获取包含293个序列的中国大陆EV71 C4亚型VP1基因数据集。ML系统发育树和贝叶斯系统发育树的C4亚型数据集包含286个C4a亚型序列和7个C4b亚型序列，其中C4b亚型毒株序列仅在1998—2003年有报道，2003年后主要流行毒株均为C4a亚型。见图1。应用Mega 7.0.20软件分析C4a和C4b亚型毒株进化距离结果表明，EV71 C4a和C4b亚型毒株VP1序列的组内平均距离分别为0.003和0.004，组间平均距离为0.006。

图1 基于中国大陆EV71 C4亚型VP1基因构建的系统进化树Fig.1 Phylogenetic tree based on VP1 gene of EV71 subtype C4 in Chinese Mainland

2.2 中国大陆C4亚型毒株的系统发育地理学分析 基于松弛分子钟和贝叶斯天际线模型的贝叶斯MCMC分析结果表明，谱系C4b亚型的序列更接近C4亚型毒株的根。时间尺度的贝叶斯系统发育分析估算EV71 C4亚型VP1基因序列的平均进化速率为2.25×10-3site／year，95%最高后验概率密度（highest posterior density，HPD）区间为（1.98～2.52）×10-3site／year。最早共同祖先（the most recent common ancestor，tMRCA）为1 989.66，95%HPD区间为1 985.42～1 993.53。华东地区的根后验概率最高（0.799），其次为华南地区（0.201）。C4a亚型毒株的平均进化速率和tMRCA分别为3.1×10-3site／year［95%HPD区间为（1.3～5.3）×10-3site／year］和1 999.03（95%HPD区间为1 996.95～2 000.70），见图2。基于C4亚型数据集的BSSVS分析结果也表明，C4亚型毒株起源于20世纪90年代的华东地区。在国内流行历史中，C4亚型毒株共有9条主要迁移路径：最早的迁移路径为1995年前后从华东输入至华南地区，5条显著迁移路径（BF≥10）分别为2001年从华东地区输出至华中、华北和西南地区，2004年从华东地区输出至东北，2007年从华东地区输出至西北地区，见表1。

图2 基于中国大陆EV71 C4亚型VP1基因构建的MCC树Fig.2 MCC tree constructed based on VP1 gene of EV71 subtype C4 in Chinese Mainland

表1 不同路径的BF和迁移率Tab.1 BF and migration rate of various paths

2.3 种群增长动态2002年以前，C4亚型病毒的种群规模相对稳定，Ne为10.76～12.18；2002—2007年，C4亚型毒株的种群规模扩张，2007年达峰值，Ne为1 120.26，随后保持较高水平，与2007年以来国内报告大规模HAMD疫情的时间点相符；2015年后种群规模未见显著变化。见图3。

2.4 自然选择及适应性分析M7-M8模型的分析结果表明，5个氨基酸位点（98、145、289、291和297）承受正选择压力。比对分析C4b亚型序列，发现VP1基因第98位点存在2种变体（E和K），145号位点存在2种变体（E和Q）；比对分析C4a亚型序列，其在98位点的变体情况与C4b亚型一致，145位均为E；比对2015—2018年毒株序列筛选到VP1基因第289位点承受正选择压力（发生A289T／V突变），见表2。其他正选择压力位点均位于VP1基因的C-末端。

图3 基于EV71 C4亚型VP1基因的BSP分析Fig.3 BSP analysis based on VP1 gene of EV71 subtype C4

表2 应用M7-M8模型进行正选择位点分析的结果Tab.2 Analysis of positive selection sites with M7-M8 models

3 讨 论

HFMD已成为全球特别是亚太地区最重要的公共卫生问题之一［24-27］。流行病学调查数据显示，20多种血清型的肠道病毒可引起HFMD，其中EV71、CV16、CV6和CV10是导致HFMD暴发和流行的主要病原体［15］，EV71感染在重症和死亡病例中占比最高［16］，成为HFMD疫苗研发的首要目标。自2015年12月以来，多家疫苗企业研发生产的EV71灭活疫苗在国内获批上市，已逐步在适龄婴幼儿人群中推广接种。

本研究采用基于聚类的贝叶斯法分析系统地理发育和种群动力变化规律，重建了EV71 C4亚型在国内的时空演化，结果发现，C4b亚型毒株仅在1998—2003年有报道，2003年后的主要流行菌株均属于C4a亚型。在亚洲其他国家或地区，如马来西亚、日本、越南和中国台湾等，EV71流行株包括C1、C2、C4和B5等亚型，且优势毒株在传播中可发生型别转换［28-29］，如在2009—2012年，中国台湾的主要基因群经历了2次变化，先是从B至C，然后从C至B［30］。提示应持续监测国内EV71流行株的型别变化，防止自然感染或疫苗免疫驱动病毒基因群转变而导致大规模HFMD的暴发。

本研究采用贝叶斯系统发育分析估算国内EV71 C4亚型的起源，可追溯至1990年前后的华东地区。华东地区是国内EV71的起源地和核心输出地，对EV71的流行和演化起到重要作用。EV71 C4亚型主要通过5条显著迁移路径，从华东地区向华中、华北、西南、东北和西北各地扩散，表明EV71在国内存在复杂的空间动力变化，在人员流动频繁的华东地区更加显著。BSP分析表明，国内EV71 C4亚型的种群自2002年以来显著扩大，并于2007年达峰值，随后保持在相对较高的水平。近年，EV71暴发或散发的文献报道相对较少，然而Ne呈较高水平，提示EV71仍有暴发的可能，应扩大疫苗在5岁以下婴幼儿人群的免疫接种。

对EV71 C4a和C4b亚型的VP1基因序列进行选择压力分析，发现5个潜在的正选择压力位点，其中第98位和145位为EV71的重要免疫原性、抗原性和毒力位点［31-33］，在C4b亚型中为145E／Q，在C4a亚型中为145E。对于C4a亚型毒株，145位点不再承受正选择压力，推测其在经历自然选择后已成为优势位点。另一个重要的压力位点是VP1基因第98位，暴露在VP1衣壳的表面，具有潜在的免疫原性［33］。另外，还有3个位点位于VP1的C-末端，特别是疫苗上市后流行株序列筛选到的A289T／V突变位点。SUN等［34］认为，A289T突变与重症HFMD密切相关：当该位点氨基酸为A时，EV71感染引起的神经症状显著增加；当氨基酸为T时，EV71神经毒性较低。A289V突变与HFMD症状的相关性尚未明确。本研究结果提示，2015年后流行株在该位点承受正选择压力，王磊等［35］曾报道，江浙沪地区2009—2014年流行株已出现A289T／V突变。该突变株能否在疫苗广泛应用后成为优势株需进一步关注。

EV71疫苗的研发和上市为国内预防和控制EV71感染导致的HFMD发挥了重要作用［36］，而肠道病毒的遗传多样性和HFMD的多病原体特征对全面防控HFMD提出了新的挑战。国内EV71疫苗接种率仍处于较低水平，为在婴幼儿人群中建立免疫屏障，应考虑将EV71疫苗纳入国家计划免疫。本研究重建了国内C4亚型的系统发生地理学和种群动态历史，筛选出VP1区承受正选择压力的5个氨基酸位点，揭示了疫苗上市前后C4亚型的进化信息，为EV71疫苗的有效应用和HFMD的防控奠定了基础。本研究基于2021年之前国内EV71流行株VP1基因开展分析，由于数据量有限，未能结合不同地区EV71疫苗覆盖率差异较大的特征进行更为细致的研究，今后HFMD的病原监测工作中，将持续关注EV71遗传变化情况，动态评估疫苗的保护效力。