柯薩奇病毒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 引起的手足口病等疫情的監測、預警和防控提供了重要信息.