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

中國(guó)食品藥品檢定研究院國(guó)家衛(wèi)生健康委員會(huì)生物技術(shù)產(chǎn)品檢定方法及其標(biāo)準(zhǔn)化重點(diǎn)實(shí)驗(yàn)室，北京 102629

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

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

1 材料與方法

1.1 數(shù)據(jù)篩選 從國(guó)家生物技術(shù)信息中心（NCBI，http：／／www.ncbi.nlm.nih.gov／）的GenBank數(shù)據(jù)庫(kù)中下載所有已知采樣日期和地理位置（國(guó)內(nèi)各省）的EV71毒株VP1基因序列（截止2021年1月1日）。應(yīng)用在線工具CD-HIT（http：／／weizhong-lab.ucsd.edu／cdhit_suite／cgi-bin／index.cgi）刪除冗余序列，相似性閾值設(shè)置為98%。應(yīng)用MAFFT v7.222軟件比對(duì)序列［20］，BioEdit v7.2.5軟件進(jìn)行編輯，截取VP1基因區(qū)段。應(yīng)用RDP v4.36軟件篩選多序列比對(duì)的重組病毒序列［21］。以已報(bào)道的不同基因型EV71毒株為參考序列，應(yīng)用iqtree 1.6.2軟件構(gòu)建最大似然（maximum likelihood，ML）系統(tǒng)發(fā)育樹，核苷酸取代模型設(shè)置為SYM+R5，設(shè)bootstrap值為1 000以評(píng)估ML系統(tǒng)發(fā)育樹拓?fù)浣Y(jié)構(gòu)的可靠性，并應(yīng)用FigTree v1.4.4軟件進(jìn)化ML系統(tǒng)發(fā)育樹的可視化，以獲得國(guó)內(nèi)分離的EV71 C4亞型VP1基因序列作為最終數(shù)據(jù)集。將選定的C4亞型毒株數(shù)據(jù)集和參考序列原始毒株BrCr（GenBank登錄號(hào)為：U22521）合并，應(yīng)用iqtree 1.6.2和MrBayes 3.2.6軟件分別構(gòu)建ML系統(tǒng)發(fā)育樹和貝葉斯系統(tǒng)發(fā)育樹，進(jìn)一步對(duì)國(guó)內(nèi)的C4亞型毒株進(jìn)行分類，并應(yīng)用Mega 7.0.20軟件計(jì)算C4a和C4b亞型毒株VP1基因的進(jìn)化距離。

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

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

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

2 結(jié) 果

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

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

2.2 中國(guó)大陸C4亞型毒株的系統(tǒng)發(fā)育地理學(xué)分析 基于松弛分子鐘和貝葉斯天際線模型的貝葉斯MCMC分析結(jié)果表明，譜系C4b亞型的序列更接近C4亞型毒株的根。時(shí)間尺度的貝葉斯系統(tǒng)發(fā)育分析估算EV71 C4亞型VP1基因序列的平均進(jìn)化速率為2.25×10-3site／year，95%最高后驗(yàn)概率密度（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ū)的根后驗(yàn)概率最高（0.799），其次為華南地區(qū)（0.201）。C4a亞型毒株的平均進(jìn)化速率和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。基于C4亞型數(shù)據(jù)集的BSSVS分析結(jié)果也表明，C4亞型毒株起源于20世紀(jì)90年代的華東地區(qū)。在國(guó)內(nèi)流行歷史中，C4亞型毒株共有9條主要遷移路徑：最早的遷移路徑為1995年前后從華東輸入至華南地區(qū)，5條顯著遷移路徑（BF≥10）分別為2001年從華東地區(qū)輸出至華中、華北和西南地區(qū)，2004年從華東地區(qū)輸出至東北，2007年從華東地區(qū)輸出至西北地區(qū)，見表1。

圖2 基于中國(guó)大陸EV71 C4亞型VP1基因構(gòu)建的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 種群增長(zhǎng)動(dòng)態(tài)2002年以前，C4亞型病毒的種群規(guī)模相對(duì)穩(wěn)定，Ne為10.76～12.18；2002—2007年，C4亞型毒株的種群規(guī)模擴(kuò)張，2007年達(dá)峰值，Ne為1 120.26，隨后保持較高水平，與2007年以來(lái)國(guó)內(nèi)報(bào)告大規(guī)模HAMD疫情的時(shí)間點(diǎn)相符；2015年后種群規(guī)模未見顯著變化。見圖3。

2.4 自然選擇及適應(yīng)性分析M7-M8模型的分析結(jié)果表明，5個(gè)氨基酸位點(diǎn)（98、145、289、291和297）承受正選擇壓力。比對(duì)分析C4b亞型序列，發(fā)現(xiàn)VP1基因第98位點(diǎn)存在2種變體（E和K），145號(hào)位點(diǎn)存在2種變體（E和Q）；比對(duì)分析C4a亞型序列，其在98位點(diǎn)的變體情況與C4b亞型一致，145位均為E；比對(duì)2015—2018年毒株序列篩選到VP1基因第289位點(diǎn)承受正選擇壓力（發(fā)生A289T／V突變），見表2。其他正選擇壓力位點(diǎn)均位于VP1基因的C-末端。

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

表2 應(yīng)用M7-M8模型進(jìn)行正選擇位點(diǎn)分析的結(jié)果Tab.2 Analysis of positive selection sites with M7-M8 models

3 討 論

HFMD已成為全球特別是亞太地區(qū)最重要的公共衛(wèi)生問(wèn)題之一［24-27］。流行病學(xué)調(diào)查數(shù)據(jù)顯示，20多種血清型的腸道病毒可引起HFMD，其中EV71、CV16、CV6和CV10是導(dǎo)致HFMD暴發(fā)和流行的主要病原體［15］，EV71感染在重癥和死亡病例中占比最高［16］，成為HFMD疫苗研發(fā)的首要目標(biāo)。自2015年12月以來(lái)，多家疫苗企業(yè)研發(fā)生產(chǎn)的EV71滅活疫苗在國(guó)內(nèi)獲批上市，已逐步在適齡嬰幼兒人群中推廣接種。

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

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

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

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