杜瑞曉，王茜，劉明琛，高帆，白玉，吳星，毛群穎，梁爭論，卞蓮蓮

中國食品藥品檢定研究院國家衛(wèi)生健康委員會生物技術產(chǎn)品檢定方法及其標準化重點實驗室，北京 102629

手足口病（hand，foot and mouth disease，HFMD）多發(fā)于5歲以下兒童，通常引起以手、足、口咽黏膜和臀部丘疹樣病變?yōu)樘卣鞯妮p癥癥狀，部分病例會引發(fā)與神經(jīng)系統(tǒng)疾病相關的重癥，如無菌性腦膜炎、腦炎和脊髓灰質(zhì)炎樣麻痹［1-5］。HFMD在亞太地區(qū)持續(xù)廣泛流行，新加坡、韓國、馬來西亞、日本、越南、中國大陸和中國臺灣均發(fā)生過大規(guī)模暴發(fā)［4，6-13］，引發(fā)全球關注。腸道病毒71型（enterovirus type 71，EV71）屬于小核糖核酸病毒科腸道病毒屬，是引起HFMD重癥和死亡病例的主要病原體之一。自1981年上海發(fā)現(xiàn)首例HFMD后，國內(nèi)大部分省份陸續(xù)報道了HFMD病例［14］。2007年后，在中國大陸中部和東部陸續(xù)發(fā)生EV71的大范圍流行，如2008年，安徽阜陽暴發(fā)HFMD疫情，造成了大量兒童感染及多個重癥病例死亡，因此，2008年5月，將HFMD列為我國法定報告的丙類傳染病，國內(nèi)企業(yè)啟動單價EV71全病毒滅活疫苗的研發(fā)并得到了迅速發(fā)展［15-18］。2015—2016年，北京科興生物制品有限公司、中國醫(yī)學科學院醫(yī)學生物學研究所、武漢生物制品研究所有限責任公司生產(chǎn)的EV71疫苗相繼獲得國內(nèi)上市許可。3家企業(yè)的EV71疫苗在Ⅲ和Ⅳ期臨床試驗中均表現(xiàn)出良好的保護效果、一致性和有效性［16］。

目前，已報道的EV71流行株主要為B和C基因型，進一步分為B0～B5和C1～C5亞型［19］。C4亞型為2000年后在國內(nèi)流行的絕對優(yōu)勢基因亞型［11］，通常包括C4a和C4b亞型分支，已上市的EV71疫苗均以C4a亞型毒株為疫苗株，已有數(shù)百萬嬰幼兒接種了EV71疫苗，有效防控了國內(nèi)EV71相關的HFMD流行，但目前對EV71疫苗上市后C4亞型進化規(guī)律鮮有報道。因此，本研究對EV71疫苗上市前后，中國大陸EV71流行株分子流行病學規(guī)律進行分析，旨在監(jiān)測EV71的遺傳進化特征，以期為制定HFMD預防和控制策略提供實驗依據(jù)。

1 材料與方法

1.1 數(shù)據(jù)篩選 從國家生物技術信息中心（NCBI，http：／／www.ncbi.nlm.nih.gov／）的GenBank數(shù)據(jù)庫中下載所有已知采樣日期和地理位置（國內(nèi)各省）的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基因區(qū)段。應用RDP v4.36軟件篩選多序列比對的重組病毒序列［21］。以已報道的不同基因型EV71毒株為參考序列，應用iqtree 1.6.2軟件構建最大似然（maximum likelihood，ML）系統(tǒng)發(fā)育樹，核苷酸取代模型設置為SYM+R5，設bootstrap值為1 000以評估ML系統(tǒng)發(fā)育樹拓撲結構的可靠性，并應用FigTree v1.4.4軟件進化ML系統(tǒng)發(fā)育樹的可視化，以獲得國內(nèi)分離的EV71 C4亞型VP1基因序列作為最終數(shù)據(jù)集。將選定的C4亞型毒株數(shù)據(jù)集和參考序列原始毒株BrCr（GenBank登錄號為：U22521）合并，應用iqtree 1.6.2和MrBayes 3.2.6軟件分別構建ML系統(tǒng)發(fā)育樹和貝葉斯系統(tǒng)發(fā)育樹，進一步對國內(nèi)的C4亞型毒株進行分類，并應用Mega 7.0.20軟件計算C4a和C4b亞型毒株VP1基因的進化距離。

1.2 系統(tǒng)發(fā)生地理學分析 利用C4亞型毒株數(shù)據(jù)集研究EV71在中國大陸的系統(tǒng)發(fā)生地理學史。根據(jù)華中、華東、華北、東北、西北、華南和西南7個區(qū)域?qū)?shù)據(jù)進行標記。應用Tracer v1.6軟件模型比較功能選擇最合適的模型，確定BEAST v1.10.5pre中貝葉斯天際線模型（Bayesian skyline）及未校正的對數(shù)正態(tài)松弛分子鐘模型，該貝葉斯天際線模型同時可用于評估相對遺傳多樣性隨時間的變化［22］。使用貝葉斯馬爾可夫鏈蒙特卡洛（Bayesian Markov chain monte carlo，MCMC）框架時，設置鏈長2億，采樣頻率10 000，并去除10%的老化樣本。設GTR+G（由JModelTest確定為最佳擬合模型）為核苷酸替換模型，對序列進行地理信息標記，并用非對稱取代模型設置貝葉斯隨機搜索變量選擇（Bayesian stochastic search variable selection，BSSVS）分析方法。剩余樣本應用Tracer v1.7.1軟件檢查所有估計的參數(shù)是否有足夠的ESS值（＞200）以達到收斂。應用TreeAnnotator軟件生成最大分支可信度（maximum clade credibility，MCC）樹，并使用Figtree v1.4.4軟件進行圖形化。另外，應用SpreaD3 v0.9.6軟件將系統(tǒng)地理學數(shù)據(jù)可視化，并計算不同地區(qū)之間的遷移路徑、后驗概率及貝葉斯因子（bayes factor，BF），一般認為BF≥10為顯著遷移路徑。

1.3 群體動態(tài)分析 為確定EV71 C4亞型在國內(nèi)的群體動態(tài)分布，分析EV71 C4亞型VP1基因遺傳多樣性隨時間的變化，對包含確切采樣日期（年份）的數(shù)據(jù)進行分析并繪制貝葉斯天際線圖（Bayesian skyline plot，BSP）。BSP圖譜可描繪有效種群規(guī)模（effective population size，Ne）隨時間的變化，揭示相對遺傳多樣性的變化規(guī)律。

1.4 選擇壓力位點分析 應用EasyCodeML軟件中位點特異性模型分析不同時期（2007年前、2008—2010年、2011—2014年和2015—2018年）EV71 C4亞型VP1基因的選擇壓力［23］。通過繪制BSP重建EV71的種群歷史，采用M7-M8模型分析EV71 C4亞型VP1基因遺傳多樣性隨時間的動態(tài)變化。

2 結 果

2.1 數(shù)據(jù)篩選及數(shù)據(jù)集建立 截止2021年1月1日，GenBank數(shù)據(jù)庫中共上傳了7 995個包含地理和時間信息的EV71序列，其中4 125個序列來自中國大陸，經(jīng)篩選最終獲取包含293個序列的中國大陸EV71 C4亞型VP1基因數(shù)據(jù)集。ML系統(tǒng)發(fā)育樹和貝葉斯系統(tǒng)發(fā)育樹的C4亞型數(shù)據(jù)集包含286個C4a亞型序列和7個C4b亞型序列，其中C4b亞型毒株序列僅在1998—2003年有報道，2003年后主要流行毒株均為C4a亞型。見圖1。應用Mega 7.0.20軟件分析C4a和C4b亞型毒株進化距離結果表明，EV71 C4a和C4b亞型毒株VP1序列的組內(nèi)平均距離分別為0.003和0.004，組間平均距離為0.006。

圖1 基于中國大陸EV71 C4亞型VP1基因構建的系統(tǒng)進化樹Fig.1 Phylogenetic tree based on VP1 gene of EV71 subtype C4 in Chinese Mainland

2.2 中國大陸C4亞型毒株的系統(tǒng)發(fā)育地理學分析 基于松弛分子鐘和貝葉斯天際線模型的貝葉斯MCMC分析結果表明，譜系C4b亞型的序列更接近C4亞型毒株的根。時間尺度的貝葉斯系統(tǒng)發(fā)育分析估算EV71 C4亞型VP1基因序列的平均進化速率為2.25×10-3site／year，95%最高后驗概率密度（highest posterior density，HPD）區(qū)間為（1.98～2.52）×10-3site／year。最早共同祖先（the most recent common ancestor，tMRCA）為1 989.66，95%HPD區(qū)間為1 985.42～1 993.53。華東地區(qū)的根后驗概率最高（0.799），其次為華南地區(qū)（0.201）。C4a亞型毒株的平均進化速率和tMRCA分別為3.1×10-3site／year［95%HPD區(qū)間為（1.3～5.3）×10-3site／year］和1 999.03（95%HPD區(qū)間為1 996.95～2 000.70），見圖2?；贑4亞型數(shù)據(jù)集的BSSVS分析結果也表明，C4亞型毒株起源于20世紀90年代的華東地區(qū)。在國內(nèi)流行歷史中，C4亞型毒株共有9條主要遷移路徑：最早的遷移路徑為1995年前后從華東輸入至華南地區(qū)，5條顯著遷移路徑（BF≥10）分別為2001年從華東地區(qū)輸出至華中、華北和西南地區(qū)，2004年從華東地區(qū)輸出至東北，2007年從華東地區(qū)輸出至西北地區(qū)，見表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 種群增長動態(tài)2002年以前，C4亞型病毒的種群規(guī)模相對穩(wěn)定，Ne為10.76～12.18；2002—2007年，C4亞型毒株的種群規(guī)模擴張，2007年達峰值，Ne為1 120.26，隨后保持較高水平，與2007年以來國內(nèi)報告大規(guī)模HAMD疫情的時間點相符；2015年后種群規(guī)模未見顯著變化。見圖3。

2.4 自然選擇及適應性分析M7-M8模型的分析結果表明，5個氨基酸位點（98、145、289、291和297）承受正選擇壓力。比對分析C4b亞型序列，發(fā)現(xiàn)VP1基因第98位點存在2種變體（E和K），145號位點存在2種變體（E和Q）；比對分析C4a亞型序列，其在98位點的變體情況與C4b亞型一致，145位均為E；比對2015—2018年毒株序列篩選到VP1基因第289位點承受正選擇壓力（發(fā)生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已成為全球特別是亞太地區(qū)最重要的公共衛(wèi)生問題之一［24-27］。流行病學調(diào)查數(shù)據(jù)顯示，20多種血清型的腸道病毒可引起HFMD，其中EV71、CV16、CV6和CV10是導致HFMD暴發(fā)和流行的主要病原體［15］，EV71感染在重癥和死亡病例中占比最高［16］，成為HFMD疫苗研發(fā)的首要目標。自2015年12月以來，多家疫苗企業(yè)研發(fā)生產(chǎn)的EV71滅活疫苗在國內(nèi)獲批上市，已逐步在適齡嬰幼兒人群中推廣接種。

本研究采用基于聚類的貝葉斯法分析系統(tǒng)地理發(fā)育和種群動力變化規(guī)律，重建了EV71 C4亞型在國內(nèi)的時空演化，結果發(fā)現(xiàn)，C4b亞型毒株僅在1998—2003年有報道，2003年后的主要流行菌株均屬于C4a亞型。在亞洲其他國家或地區(qū)，如馬來西亞、日本、越南和中國臺灣等，EV71流行株包括C1、C2、C4和B5等亞型，且優(yōu)勢毒株在傳播中可發(fā)生型別轉(zhuǎn)換［28-29］，如在2009—2012年，中國臺灣的主要基因群經(jīng)歷了2次變化，先是從B至C，然后從C至B［30］。提示應持續(xù)監(jiān)測國內(nèi)EV71流行株的型別變化，防止自然感染或疫苗免疫驅(qū)動病毒基因群轉(zhuǎn)變而導致大規(guī)模HFMD的暴發(fā)。

本研究采用貝葉斯系統(tǒng)發(fā)育分析估算國內(nèi)EV71 C4亞型的起源，可追溯至1990年前后的華東地區(qū)。華東地區(qū)是國內(nèi)EV71的起源地和核心輸出地，對EV71的流行和演化起到重要作用。EV71 C4亞型主要通過5條顯著遷移路徑，從華東地區(qū)向華中、華北、西南、東北和西北各地擴散，表明EV71在國內(nèi)存在復雜的空間動力變化，在人員流動頻繁的華東地區(qū)更加顯著。BSP分析表明，國內(nèi)EV71 C4亞型的種群自2002年以來顯著擴大，并于2007年達峰值，隨后保持在相對較高的水平。近年，EV71暴發(fā)或散發(fā)的文獻報道相對較少，然而Ne呈較高水平，提示EV71仍有暴發(fā)的可能，應擴大疫苗在5歲以下嬰幼兒人群的免疫接種。

對EV71 C4a和C4b亞型的VP1基因序列進行選擇壓力分析，發(fā)現(xiàn)5個潛在的正選擇壓力位點，其中第98位和145位為EV71的重要免疫原性、抗原性和毒力位點［31-33］，在C4b亞型中為145E／Q，在C4a亞型中為145E。對于C4a亞型毒株，145位點不再承受正選擇壓力，推測其在經(jīng)歷自然選擇后已成為優(yōu)勢位點。另一個重要的壓力位點是VP1基因第98位，暴露在VP1衣殼的表面，具有潛在的免疫原性［33］。另外，還有3個位點位于VP1的C-末端，特別是疫苗上市后流行株序列篩選到的A289T／V突變位點。SUN等［34］認為，A289T突變與重癥HFMD密切相關：當該位點氨基酸為A時，EV71感染引起的神經(jīng)癥狀顯著增加；當氨基酸為T時，EV71神經(jīng)毒性較低。A289V突變與HFMD癥狀的相關性尚未明確。本研究結果提示，2015年后流行株在該位點承受正選擇壓力，王磊等［35］曾報道，江浙滬地區(qū)2009—2014年流行株已出現(xiàn)A289T／V突變。該突變株能否在疫苗廣泛應用后成為優(yōu)勢株需進一步關注。

EV71疫苗的研發(fā)和上市為國內(nèi)預防和控制EV71感染導致的HFMD發(fā)揮了重要作用［36］，而腸道病毒的遺傳多樣性和HFMD的多病原體特征對全面防控HFMD提出了新的挑戰(zhàn)。國內(nèi)EV71疫苗接種率仍處于較低水平，為在嬰幼兒人群中建立免疫屏障，應考慮將EV71疫苗納入國家計劃免疫。本研究重建了國內(nèi)C4亞型的系統(tǒng)發(fā)生地理學和種群動態(tài)歷史，篩選出VP1區(qū)承受正選擇壓力的5個氨基酸位點，揭示了疫苗上市前后C4亞型的進化信息，為EV71疫苗的有效應用和HFMD的防控奠定了基礎。本研究基于2021年之前國內(nèi)EV71流行株VP1基因開展分析，由于數(shù)據(jù)量有限，未能結合不同地區(qū)EV71疫苗覆蓋率差異較大的特征進行更為細致的研究，今后HFMD的病原監(jiān)測工作中，將持續(xù)關注EV71遺傳變化情況，動態(tài)評估疫苗的保護效力。