埃可病毒18 型構象表位的生物信息學研究

王 惠 ,賈永濤 ,方雨露 ,董長征*

（1.寧波大學 醫學院,浙江 寧波 315211;2.浙江省病理生理學技術研究重點實驗室,浙江 寧波 315211）

人腸道病毒屬小RNA 病毒科(Picornaviridae),腸道病毒屬(Enterovirus),包括腸道病毒 A 種(EV-A)、B 種(EV-B)、C 種(EV-C)和D 種(EV-D).EV-A 代表性血清型包括腸道病毒A 組71 型(Enterovirus A71,EV-A71)和柯薩奇病毒A 組16 型(Coxsackievirus A16,CV-A16)等;EV-B 代表性血清型包括埃可病毒18 型(Echovirus 18,E18)等;EV-C 代表性血清型包括脊髓灰質炎病毒1 型(Polioviurs 1,PV1)、PV2 和PV3 等;EV-D 代表性血清型為EV-D68.E18 可引起手足口病、病毒性腦膜炎和急性胃腸炎等疾病[1-4],通常感染后癥狀較為輕微,但嚴重時也會危及生命[5].自1955 年在美國首次分離病毒毒株后,E18 感染鮮有報道,直到21 世紀E18 開始在全球流行,并且引起多起較大規模的疫情[6-9].2014—2016 年美國腸道病毒監測報告顯示,E18 在最常引起疫情的腸道病毒中排第4 位[10].Chen 等[11]在2015—2016 年從我國6 個省份(山東、河北、山西、黑龍江、江蘇和云南)的手足口病、病毒性腦炎和病毒性腦膜炎病例中分離出34 株E18,測定其結構蛋白(Viral Protein,VP)VP1 序列,并完成了其中6株的全基因組測序.2015—2020年,E18 在我國廣東[4,12]、河北[13]、山東[14]和云南[15]等省大量檢出,提示我國存在E18 暴發風險.

與其他腸道病毒相似,E18有一個長度約為7.4 kb 的單股正鏈RNA基因組,其病毒衣殼由60 個亞單位構成[16-18].每個非對稱亞單位由VP1～VP4 構成,其中VP4 位于衣殼內表面,VP1～VP3 位于衣殼外表面,后者是病毒構象表位的所在區域.抗原表位的確定對掌握病毒的致病機制[19-20]、監測病毒的變異和進化[21]、研發抗病毒藥物[22-23]和疫苗[24-25]都有重要作用.然而,目前尚無針對E18 表位的研究報道.

本研究首先利用前期發展的生物信息學算法[17-18]對E18 構象表位進行系統性預測,分析表位的分布規律;然后基于Nextstrain 平臺構建E18 的分子進化樹,確定進化分支和決定進化分支的標志性氨基酸突變(分支突變);最后探索分支突變與表位以及病毒受體結合位置的關系,分析表位在E18 分子進化過程中的作用.研究結果可為E18 的監測和預警提供參考.

1 材料與方法

1.1 E18 結構蛋白的序列和結構特征分析

從RCSB PDB 數據庫[26](https://www.rcsb.org)和NCBI Nucleotide數據庫[27](https://www.ncbi.nlm.nih.gov/nuccore)中分別下載E18 天然成熟顆粒的結構蛋白的三維結構文件(PDB ID: 6HBG[16];病毒毒株名: Metcalf,為簡便以下用PDB ID 指代病毒毒株)和其氨基酸序列(Accession ID: AAL37163).將PDB 文件和氨基酸序列導入在線工具ESPript 3.0[28](http://espript.ibcp.fr)中注釋二級結構信息.使用PyMOL[29]繪制病毒衣殼的表面結構圖.所有軟件和在線工具都采用默認參數.

1.2 E18 構象表位的生物信息學預測

實驗室前期在Borley 算法[30]基礎上發展了人腸道病毒構象表位的生物信息學預測算法,并成功地應用于腸道病毒A 種(EV-A71、CV-A16[17])和D 種(EV-D68[31])的表位預測.預測算法主要包括以下3 步(詳見文獻[17-18]):

(1)在PDB 文件中刪除VP4 并生成復合鏈(圖1),將復合鏈整體看作一個蛋白質,代替單個結構蛋白進行表位預測.

圖1 E18 復合鏈和衣殼表面結構

(2)利用3 個表位預測工具Epitopia[32](http://epitopia.tau.ac.il)、Ellipro[33](http://tools.iedb.org/ellipro)以及DiscoTope[34](http://www.cbs.dtu.dk/services/DiscoTope)分別預測復合鏈的表位.三者的閾值均使用默認參數(0.174、0.3、-10.7).采用投票法,將同時被3 種工具預測為表位的氨基酸殘基作為一致性表位.

(3)提取中心鏈Chain 1(一個由VP1～VP3 組成的亞單位)上的一致性表位,并篩選處于病毒衣殼相對暴露面的殘基(其Cα到衣殼中心的距離超過所有Cα到衣殼中心的平均距離),獲得最終的預測結果.

1.3 基于Nextstrain 的E18 分子進化分析

Nextstrain[35]是免費的在線(https://nextstrain.org)或本地開源平臺,針對病毒基因/基因組數據,通過強大的分子進化分析和數據可視化功能,追蹤病毒的分子進化和時空傳播路線.本研究通過本地化部署Nextstrain平臺,利用Nextstrain平臺的augur 病原體生物信息分析包,分別對E18 的VP1和全基因組序列進行分子進化分析,步驟如下:

(1)從VIPR 病毒數據庫[36](https://www.viprbrc.org/brc)和NCBI Nucleotide 數據庫[27](https://www.ncbi.nlm.nih.gov/nuccore)中下載所有長度＞700 bp的E18 VP1 序列和長度＞7 000 bp 的基因組序列,剔除信息不全序列(無病毒毒株分離國家或時間信息等).采用MAFFT 7[37]在線平臺(https://mafft.cbrc.jp/alignment/server)進行多重序列比對,人工剔除異常序列(例如多重序列比對后通過人工核查發現存在蛋白質翻譯錯誤,序列中存在超過3 個模糊堿基).最終獲得285 條E18 的VP1 序列和63 條基因組序列.

(2)利用augur 病原體生物信息分析包進行生物信息學和分子進化分析.align 命令調用內嵌的MAFFT 軟件[37]進行多重序列比對,tree 命令調用IQ-tree工具[38]構建分子進化樹(最大似然法),refine命令調用TreeTime 工具[39]構建基于時間尺度的分子進化樹并對進化樹進行優化,clade 命令標注進化分支(clade)和確定分支突變.由于表位僅分布在VP1～VP3 上,因此只標注了VP1～VP3 上的分支突變.最后利用Nextstrain 平臺的auspice 工具進行基于JavaScript 技術的網頁互動可視化展示.

(3)將獲得的分支突變與表位、病毒受體結合位置利用RIVEM 工具[40-41]繪制在E18 衣殼表面結構圖上,即足跡圖(roadmap),分析分支突變與表位以及受體結合之間的關系.新生兒 Fc 受體(Neonatal Fc Receptor,FcRn)是埃可病毒共用的脫衣殼受體,埃可病毒借助FcRn 感染宿主細胞并釋放病毒基因組.雖然E18 的結合位置尚未確定,但E6 和E30 與FcRn 的結合位置已經確定[19-20],可以通過序列比對的方式確定E18 的模擬結合位置.

2 結果

2.1 E18 的結構蛋白特征和構象表位預測結果

E18 病毒顆粒6HBG 的結構蛋白VP1、VP2 和VP3 的長度分別為287、260 和239 個氨基酸殘基,由8 條反向平行的β 鏈(βB-βI)構成β 桶(β-barrel),鏈之間的部分為環區(loop),兩端為 N-端(Nterminus)和C-端(C-terminus),如圖1和圖2所標注.E18 的衣殼表面結構如圖1 所示,衣殼表面有峽谷(canyon)、峽谷兩側的“邊緣”(rim)、“平臺”(puff)和“突起”(knob)等結構特征以及五倍軸、三倍軸和二倍軸等三維結構標記.

E18 構象表位預測結果見表1 及圖1 和圖2.E18 共有27 個氨基酸殘基預測為表位,分布在VP1(BC 環、DE 環、HI 環和C-端)、VP2(EF 環和HI 環)和VP3(N-端knob 區域、BC 環和C-端).與EV-A[17]和EV-D[31]相似,E18 的構象表位也聚集成三簇(表1 和圖1): site 1、site 2 和site 3,分別位于峽谷的“北側邊緣”區域、峽谷南側的“平臺”區域、峽谷南側的“突起”區域和三倍軸區域.其中VP1 BC 環和C-端、VP2 EF 環是E18 表位的主要構成區域.

表1 E18 構象表位的預測結果

圖2 E18 結構蛋白的一級結構和二級結構注釋

2.2 E18 的分子進化分析

利用Nextstrain 平臺的augur 病原體生物信息分析包,分別構建了基于E18基因組和VP1序列的時間尺度的分子進化樹,分別簡稱為基因組進化樹(圖3(a))和VP1進化樹(圖3(b)),兩者具有一致的拓撲結構,都分為A、B、C 三個進化分支,C 分支又分為C1 和C2 兩個子分支,原株Metcalf 不屬于任何一個分支.以序列數較多的VP1 進化樹為例,大約1946 年Metcalf 從A、B、C 三個分支的共同祖先A～C 中分歧出來,并于1955 年在美國首次分離.1970 年分子進化樹首先分歧出A 分支;1979 年又分歧出B 和C 分支;最后在1989 年C 分支分歧成C1 和C2 兩個子分支.C 分支尤其是C2 分支是目前流行的病毒分支.A 分支由2 株來自中國的病毒毒株構成.B分支有1株來自埃塞俄比亞,其余都來自印度.C1 分支主要來自法國、德國、俄羅斯、瑞典和澳大利亞.C2 分支全球分布廣泛,其中分離病毒毒株數最多的國家為中國(168 株)、法國(14 株)、澳大利亞(14 株)、日本(13 株)和美國(7 株).

每個進化分支都有若干個決定進化分支的分支突變.以包含VP1～VP3 的基因組進化樹為例(圖3(a)),A 分支的分支突變為VP1 C-端的T271A 和D275E,B 分支的分支突變為VP1 C-端的K271S、A285V 和T286S.C1 的分支突變包括VP1 C-端的D275E 以及R6K、VP2 EF 環的S159P、VP3 的V2I和H182N.C2 的分支突變包括VP1 C-端的G257S、A262V、A285V 以及I42L 和I92V,VP2 的T74S,VP3 knob 的V58I 以及N11T.從圖3(a)可以發現,分子進化樹上多數(61.2%)的分支突變都位于表位處,即分支進化伴隨著表位處的氨基酸突變,尤其是VP1 C-端是突變熱點區域.此外,還有VP2 的EF 環.VP1 進化樹(圖3(b))也表現為高度一致,絕大多數分支突變在兩種進化樹上分布一致,但由于序列數目差異較大(285 比63),個別分支突變產生了差異,例如VP1 進化樹B～C 分支上的突變D129E 成為基因組進化樹C 分支上的突變.

圖3 E18 的分子進化樹

利用足跡圖(圖4)將上述分支突變與表位、受體結合位置都標注在E18 的衣殼表面結構圖上,分析三者之間的關系可以發現,大多數E18衣殼表面的分支突變(黑色實線圈)都位于表位上或附近,而且表位的3 個site 均被突變覆蓋;受體FcRn 模擬結合位置(白色虛線橢圓)僅有2 個突變.這提示表位處的突變產生了新的進化分支,但盡量避開受體結合位置,以免影響病毒與宿主細胞的結合能力.這兩個突變位點對受體結合能力的影響需要進一步實驗研究.

圖4 E18 足跡圖

3 討論

病毒構象表位測定的金標準是冷凍電鏡技術[19-20,23],但冷凍電鏡技術門檻高,現階段難以常規應用.免疫逃避實驗是傳統測定構象表位的常用方法[42-43],但只能通過突變體確定部分表位,同時也費時費力.生物信息學算法具有高通量和相對準確的特點,能夠為實驗性研究提供候選靶標,輔助實驗性研究.傳統的生物信息學表位預測算法具有普適性,但如果直接用來預測腸道病毒的結構蛋白的表位,由于未考慮“結構蛋白嵌在衣殼中”這個生物學結構特征,預測結果必然有大量的假陽性[17,30,44].本研究使用實驗室前期發展的表位預測算法,其最大優勢在于將生物信息學算法和腸道病毒的衣殼結構特征密切結合,大大提高了算法的可靠性[17].在不同腸道病毒種型(EV-A 和EV-D)表位預測的實際應用中,通過對比實驗結果可以確認算法具有較高的準確性[17,31].當然,生物信息學算法不可避免地存在自身的缺陷和限制.例如,算法僅能對已測定三維結構的腸道病毒進行表位預測;算法的精確性受到三維結構精確性的嚴重影響;算法的預測原理基于表位的病毒學假設(衣殼表面突出區域更易與抗體或受體結合),而這種假設需要在實踐中通過實驗性研究不斷驗證和改進,算法也隨之更新.

作為EV-B 代表性血清型之一的E18,它的構象表位與EV-A(EV-A71 和CV-A16[17])和EV-D(EV-D68[31])一樣呈現三簇分布規律.VP1 BC 環和C-端、VP2 EF 環是腸道病毒共有的表位,VP1 GH環是EV-A 特有的表位,VP1 DE 環是EV-D 特有的表位,E18 則未發現特異性表位.即使表位分布高度相似,但表位上的氨基酸突變決定了不同血清型腸道病毒的抗原性差異[17,45].與EV-D68[31]相似,VP1 C-端和VP2 EF 環是E18 的突變熱點區域,每一個E18 進化分支的形成都伴隨著氨基酸突變.目前尚未見能中和E18 的單克隆抗體(單抗)的相關報道,但同屬EV-B 的E30 已報道有兩個單抗4B10和6C5[46].單抗4B10 主要結合E30 峽谷南側的VP1 C-端和VP2 EF 環,單抗6C5 結合位點包括峽谷北側的VP1 BC 環、DE 環、EF 環和HI 環,而這些位點大多數都是E18 的表位區域,提示這些表位區域易與單抗結合.

本研究分別基于E18 的VP1 和基因組序列構建了時間尺度的分子進化樹.VP1 進化樹包含的序列數較多(285 條),能夠更加準確地確定進化關系和推斷分歧時間,但VP1 進化樹僅能獲得VP1 上的分支突變.基因組進化樹雖然序列數較少(63 條),但能獲得VP1～VP3 上的分支突變.基因組進化樹和VP1 進化樹具有高度一致的拓撲結構和VP1 上的分支突變,一方面證明本研究確定的分支突變高度可靠,另一方面基因組進化樹補充了VP1 進化樹缺乏的VP2 和VP3 上的分支突變,尤其是VP2 EF 環和VP3 knob 這兩個重要表位區域.參照Chen 等[11,13]研究結果,將E18 分為A、B、C 三個主要分支,原株Metcalf 未劃入任何分支.相較于A 分支和B 分支,C 分支是E18 的主要流行分支,96.8%分離的病毒毒株都屬于C 分支.C 分支分為C1 分支和C2 分支,其中C1 分支在2013 年后未再流行,而C2 分支病毒毒株廣泛分離于近年來的病毒性腦膜炎和手足口病疫情中,提示C2 分支可能具有較強的傳染力和毒力,但尚未有E18傳染力和毒力的進一步研究報道.分子進化樹(圖3)顯示,E18 進化分支的形成伴隨著表位處的氨基酸突變,尤其是VP1 C-端和VP2 EF 環是突變熱點區域.之前的研究發現[4-5,11]多個VP1 上的多個重要突變位點,包括R6K、N10D、R84N、M104L、Y215F、I216V、V262T/A 和D275E 等.VP1 BC 環是E18重要的表位區域(圖1 和圖2).有研究發現[17,31],VP1 BC 環是EV-A 和EV-D 共有的重要表位.BC環上的R84N 是C 分支的分支突變,A 和B 分支均為R,幾乎所有C 分支均為N(2 株為S),進一步驗證了表位VP1 BC 環對腸道病毒的重要性.K6R 和T262V 是C2 的分支突變,D275E 則是C1 的分支突變(圖3).從Nextstrain平臺構建的VP1 進化樹上可以清晰看到,D10N、Y215F、I216V、S257G 和H287R 都是C2 分支中國株的標志性氨基酸突變,M104L 則是其中一大簇病毒毒株上的標志性突變.E18 的分子進化分析,一方面從另一個角度間接地驗證了表位預測的準確性,另一方面也提示表位處的重要氨基酸突變及相應病毒毒株的時空傳播路線是腸道病毒監測和預警分析的重點.

對于流感病毒和冠狀病毒,由于接種疫苗和廣泛感染產生群體免疫,使得病毒處于正選擇進化壓力下,因此不僅表位突變頻繁,而且受體結合區域也存在廣泛突變,從而產生免疫逃避[47-48].而對于腸道病毒,由于新生兒不斷補充易感群體[49],使得多數腸道病毒并沒有面臨顯著的進化壓力,因此相對突變速率要慢一些(E18 的VP1 約為5.6×10-3替換·(位點·a)-1,與其他腸道病毒相似),在表位區域存在較弱的正選擇壓力,其他區域則偏向中性進化[50].E18 受體結合區域也更加保守,其足跡圖提供了佐證.但隨著E18 的廣泛流行,人群血清中普遍存在中和抗體,病毒面臨越來越強的正選擇壓力,表位也會突變得更加頻繁,需要密切監測其抗原性是否會發生較大改變,出現類似EVA71 那樣的大規模疫情[51].

生物信息學在流感病毒的優勢病毒毒株預測方面取得了巨大成功[52-54],為流感疫苗的研發和準備提供了重要的技術支撐.對腸道病毒表位生物信息學的研究,可為進一步通過實驗鑒定構象表位、病毒的監測和預警以及抗病毒藥物和疫苗的研發提供重要支持.