彭小川，靳遠(yuǎn)，胡明達(dá)，周江林，周靜，岳俊杰，任洪廣，梁龍

軍事科學(xué)院 軍事醫(yī)學(xué)研究院 生物工程研究所，北京100071

根據(jù)其在禽類中的致病性，可將禽流感病毒分為高致病性禽流感和低致病性禽流感[1]。大多數(shù)禽流感病毒都是低致病性的，而且在多數(shù)鳥類宿主中不出現(xiàn)病癥。但是在抗原漂變的作用下，一些亞型如H5、H7進(jìn)化成為高致病性毒株，并且在其宿主內(nèi)保持很高的致病率[2]。對(duì)于H5N1來說，高致病性毒株居多，但也有一些零星的低致病性毒株出現(xiàn)在韓國(guó)、日本和北美地區(qū)[3-4]。高致病性H5N1病毒于1996年在廣東省的一只家鵝中首次被分離鑒定出來[5]，從那時(shí)起，全球就不斷有疫情在家禽、野鳥和人群中暴發(fā)[6]。全球首例人感染H5N1病毒的事件發(fā)生在1997年的香港，共造成18人感染，6人死亡[7]。其后5年都沒有人被感染的情況發(fā)生，直到2003年，又有2個(gè)病例在中國(guó)大陸被檢測(cè)出來[8]。從2003年起，高致病性H5N1病毒感染人類的事件便不斷地在中國(guó)[9]、東南亞[10]、西亞[11]、非洲[12]出現(xiàn)，至今已造成900多人感染。歷史上最大規(guī)模的1918西班牙流感就是由禽流感病毒造成的，全球共計(jì)有四五千萬人喪生[13]。因此，H5N1病毒突然跨越宿主屏障，并在最近十幾年持續(xù)感染人類，引發(fā)了全世界對(duì)下一次流感大暴發(fā)的恐慌。

雖然H5N1病毒還不具備在人和人之間有效傳播的能力，但是其突變和適應(yīng)其他宿主的巨大潛力，仍然值得關(guān)心和研究。無論是病毒的毒力、逃避免疫的能力還是跨宿主傳播的能力等，都和病毒的適應(yīng)性有關(guān)，而病毒的適應(yīng)性與其進(jìn)化過程中產(chǎn)生的突變密切相關(guān)。大多數(shù)對(duì)于H5N1病毒突變的研究，都是針對(duì)某一個(gè)或某幾個(gè)位點(diǎn)的作用展開的，如Peng等發(fā)現(xiàn)H5N1病毒血凝素（hemagglutinin，HA）上的突變K193T能夠增強(qiáng)人型受體的結(jié)合[14]；Qin等發(fā)現(xiàn)HA上的突變E190G能夠減少H5N1病毒的復(fù)制，還能明顯減弱病毒的毒力[15]；Han等證明HA上K193E和G225E這2個(gè)突變可以協(xié)同減弱H5N1的毒力且不會(huì)增強(qiáng)受體結(jié)合能力[16]。

類似的研究還有很多，雖然專一性強(qiáng)，但是缺乏對(duì)病毒適應(yīng)性進(jìn)化過程中突變規(guī)律的整體把控。我們以H5N1病毒HA蛋白為研究對(duì)象，利用生物信息學(xué)手段進(jìn)行系統(tǒng)發(fā)生學(xué)分析，建立了一種基于株頻率和枝頻率以及標(biāo)記分枝的突變回溯方法，找到了導(dǎo)致H5N1病毒HA蛋白適應(yīng)性進(jìn)化的關(guān)鍵突變。進(jìn)一步分析這些突變的頻率變化及頻次特征，發(fā)現(xiàn)這些突變位點(diǎn)多數(shù)為正選擇位點(diǎn)，有許多落在抗原表位上，很好地印證了它們和病毒適應(yīng)性進(jìn)化的關(guān)聯(lián)性。此外，基于我們建立的方法和評(píng)價(jià)手段，我們還找到了以往研究中沒有發(fā)現(xiàn)的一些關(guān)鍵突變，提示這些突變及所在位點(diǎn)可能與病毒的適應(yīng)性進(jìn)化有關(guān)，可以為實(shí)驗(yàn)提供指導(dǎo)。

1 材料與方法

1.1 數(shù)據(jù)下載和預(yù)處理

從GISAID（Global Initiative on Sharing All Influenza Data）數(shù)據(jù)庫（www.gisaid.org）下載具有完整基因組的H5N1序列，從中提取出1271條HA序列。去除測(cè)序錯(cuò)誤和有明顯缺失的序列后，為了減少采樣偏差，我們將具有相同觀測(cè)時(shí)間、相同地理信息、相同宿主和基因型的序列只保留1條，最終得到一個(gè)包含1002條HA序列的數(shù)據(jù)集。

1.2 多序列比對(duì)和系統(tǒng)發(fā)生樹構(gòu)建

利用MAFFT v7.402[17]對(duì)HA序列集進(jìn)行多序列比對(duì)，然后用MEGA 7.0.25[18]手動(dòng)整理比對(duì)結(jié)果并且只保留HA序列的編碼區(qū)部分。接著，利用RaxML v8.2.4[19]，選用GTRGAMMA模型進(jìn)行1000次bootstraps，構(gòu)建H5N1病毒HA序列的最大似然樹，并選取毒株A/turkey/England/50-92/1991對(duì)系統(tǒng)發(fā)生樹定根。

1.3 株頻率的定義

在甲型流感病毒H3N2預(yù)測(cè)模型的研究中，Luksza等提出了“株頻率”的概念[20]?？紤]到H5N1是禽流感病毒，我們?cè)诒狙芯恐袑?duì)這個(gè)概念做出了一定的改進(jìn)。

首先，我們將每一個(gè)毒株映射到一個(gè)信息向量=(n,g,t,h,l)上，向量的5個(gè)元素分別代表序列的名字、基因型、觀測(cè)年份、宿主及地理信息。不難發(fā)現(xiàn)，我們預(yù)處理過的數(shù)據(jù)和信息向量是一一映射的。然后，我們將每一個(gè)毒株用一個(gè)指標(biāo)k來標(biāo)記，在年份t中觀測(cè)到的毒株的集合用S(t)來標(biāo)記。對(duì)于每一個(gè)毒株k，我們可以定義一個(gè)多樣性dk表示具有完全相同的基因型gk和觀測(cè)年份tk，但宿主信息hk或地理信息lk不同的HA序列的數(shù)量，即具有完全相同的信息向量的序列只計(jì)算1次。這樣，便可以計(jì)算毒株k的株頻率fk為：

1.4 分枝和枝頻率的定義

根據(jù)圖論知識(shí)，系統(tǒng)發(fā)生樹是二元的，由分枝和結(jié)點(diǎn)組成。結(jié)點(diǎn)分為葉結(jié)點(diǎn)和內(nèi)部結(jié)點(diǎn)。顯然，葉結(jié)點(diǎn)就是序列數(shù)據(jù)集在樹上的映射，而它們的父結(jié)點(diǎn)或者更早的祖先結(jié)點(diǎn)構(gòu)成內(nèi)部結(jié)點(diǎn)。內(nèi)部結(jié)點(diǎn)和它們的子代構(gòu)成的子樹稱為進(jìn)化樹的“分枝（clade）”。對(duì)于每個(gè)cladeν，所有該分枝上屬于年份t的毒株集合表示為Sν(t)，Sν(t)?S(t)，并且對(duì)于進(jìn)化樹。由此，cladeν在年份t的“枝頻率”可以定義為：

進(jìn)一步，我們可以得到一個(gè)關(guān)于cladeν的反映各年份間頻率變化的頻率向量：

其中ti∈Ω，Ω是全部年份的集合。

1.5 標(biāo)記分枝的定義

尋找H5N1病毒HA進(jìn)化過程中的突變是我們研究的重點(diǎn)。為此，我們使用PAUP 4.0[21]重建系統(tǒng)發(fā)生樹內(nèi)部結(jié)點(diǎn)的最大似然序列，得到結(jié)點(diǎn)間的連接關(guān)系和衍征信息。然后我們定義“標(biāo)記分枝”如下：對(duì)于每一個(gè)年份t，滿足的最大子樹叫做年份t的標(biāo)記分枝，其中N代表標(biāo)記分枝的最小規(guī)模，可根據(jù)需要細(xì)化的程度靈活選取。

1.6 氨基酸頻率圖表及頻率變換

為了評(píng)估尋找到的突變，我們引入他人研究中關(guān)于氨基酸頻率和頻率變換的概念[22]。對(duì)于年份t，如果第j位的氨基酸是ak的序列數(shù)量是n(t,j,ak)，那么第j位的氨基酸是ak的頻率為f(t,j,ak)=n(t,j,ak)/N(t)。因此，各個(gè)突變位點(diǎn)的氨基酸頻率是關(guān)于年份t的函數(shù)，我們可以將它們儲(chǔ)存在一個(gè)氨基酸頻率圖表中。

對(duì)于在j位點(diǎn)發(fā)生的從ak到am的突變，它們?cè)谀攴輙的氨基酸頻率分別為f(t,j,ak)和f(t,j,am)。我們稱j位點(diǎn)氨基酸ak和am在年份t和t+1之間存在頻率變換，如果它們滿足如下2個(gè)條件：

條件①表明在年份t和t+1，位點(diǎn)j的氨基酸主要由ak和am構(gòu)成；條件②用來檢驗(yàn)位點(diǎn)j發(fā)生突變后氨基酸頻率是否有明顯的改變。

如果j位點(diǎn)氨基酸ak和am在年份t和t+1之間存在頻率變換，而且突變后占主導(dǎo)的氨基酸在樹上固定或者基本固定，即固定率，其中n(ν,j,am)是cladeν上位點(diǎn)j是氨基酸am的序列數(shù)量，N(ν)是cladeν上全部序列的數(shù)量，我們稱這是一個(gè)有效變換，反之稱為無效變換。

2 結(jié)果

2.1 HA系統(tǒng)發(fā)生樹及進(jìn)化過程中的突變

基于定義的株頻率、枝頻率和標(biāo)記分枝，我們建立了一個(gè)基于回溯的尋找H5N1病毒HA蛋白進(jìn)化過程中突變的方法。

株頻率計(jì)算結(jié)果見表1。根據(jù)計(jì)算出的株頻率，我們得到樹上全部分枝的頻率向量，選取N=3，從葉結(jié)點(diǎn)向前回溯，計(jì)算出全部標(biāo)記分枝共194個(gè)。分別算出根結(jié)點(diǎn)到194個(gè)標(biāo)記分枝的進(jìn)化路徑和其上包含的氨基酸突變，經(jīng)過篩選去重，最終得到突變435個(gè)。系統(tǒng)發(fā)生樹及進(jìn)化路徑上關(guān)鍵位置的部分突變?nèi)鐖D1。

2.2 氨基酸頻率圖表

氨基酸突變位點(diǎn)的頻率圖表反映了該位點(diǎn)氨基酸頻率隨年份的變化情況，可以方便直觀地體現(xiàn)氨基酸替換的動(dòng)態(tài)過程。圖2的a、b顯示，171、279位點(diǎn)（H5參考系）都經(jīng)歷了明顯的氨基酸替換過程和頻率起伏，符合正選擇的特征，而已有研究也表明這2個(gè)位點(diǎn)都是正選擇位點(diǎn)[23]。相對(duì)地，圖2的c、d反映出498和98位點(diǎn)在整個(gè)時(shí)間軸上沒有明顯的氨基酸頻率波動(dòng)，但在Ven?kata等的研究中卻將其推斷為正選擇位點(diǎn)[23]。另外，他們認(rèn)為圖2的e、f中61和100位點(diǎn)不是正選擇位點(diǎn)，但是從2個(gè)位點(diǎn)的氨基酸替換情況和頻率波動(dòng)幅度以及2個(gè)位點(diǎn)位于HA抗原表位的事實(shí)[24]，我們可以推斷出它們都是正選擇位點(diǎn)。由此看出，頻率圖表能很好地體現(xiàn)自然選擇的過程，對(duì)于正選擇位點(diǎn)能做出有效判斷。

表1 各年份毒株數(shù)量及株頻率分布表

2.3 頻率變換及有效變換

氨基酸頻率圖表清楚地展示了頻率的跳躍式變化，也就是突變位點(diǎn)氨基酸明顯的取代更迭過程，我們將這種年份間的主要氨基酸替換過程稱為“頻率變換”。如果突變后的氨基酸固定或基本固定了至少1年，我們就稱其為“有效變換”，否則叫做“無效變換”。通過年份間最大頻率變化和突變后固定情況，共計(jì)算出有效變換79個(gè)，無效變換97個(gè)。79個(gè)有效變換發(fā)生在36個(gè)位點(diǎn)上，其中20個(gè)位點(diǎn)是已知抗原表位[24]。因此，我們認(rèn)為免疫選擇是突變演替過程的重要推動(dòng)力。

在計(jì)算氨基酸突變后的固定情況時(shí)，我們得到一張全面的氨基酸頻率時(shí)間圖（表2），從中發(fā)現(xiàn)了一個(gè)有趣的現(xiàn)象：在大多數(shù)年份，都有多個(gè)突變位點(diǎn)的氨基酸同時(shí)發(fā)生了進(jìn)化過程中的第一次固定，包含的氨基酸突變共有92個(gè)，分別出現(xiàn)在45個(gè)位點(diǎn)上，其中包含22個(gè)正選擇位點(diǎn)[23]。突變固定隨時(shí)間的連續(xù)性以及正選擇位點(diǎn)的富集證明病毒在進(jìn)化過程中受到持續(xù)不斷的正選擇作用。

2.4 高頻次氨基酸突變

在計(jì)算出的全部氨基酸突變中，我們發(fā)現(xiàn)一個(gè)有趣的現(xiàn)象，有很多位點(diǎn)的氨基酸在進(jìn)化過程中發(fā)生了超過5次突變。經(jīng)統(tǒng)計(jì)，共有29個(gè)位點(diǎn)氨基酸突變次數(shù)超過5次，有14個(gè)位點(diǎn)在HA的抗原表位上。其中文獻(xiàn)報(bào)道過的正選擇位點(diǎn)有16個(gè)[23]，利用本文的氨基酸頻率圖表及有效變換推測(cè)出的正選擇位點(diǎn)有20個(gè)?？梢园l(fā)現(xiàn)，使用我們的方法不但能找到文獻(xiàn)記載的全部正選擇位點(diǎn)，而且抗原表位與我們找到的正選擇位點(diǎn)重合度更高（詳見表3）。隨后我們選取部分高頻次突變標(biāo)記在系統(tǒng)發(fā)生樹上（圖1），發(fā)現(xiàn)高頻次突變很多都發(fā)生在病毒進(jìn)化過程中的關(guān)鍵路徑和關(guān)鍵位點(diǎn)上。此外，將這些高頻次突變放在HA的三維結(jié)構(gòu)中進(jìn)行分析，我們還得到如下幾個(gè)有趣的結(jié)果：①154（N，H，Q，L）、156（K，R，N，T等）、172（A，T）3個(gè)位點(diǎn)近年持續(xù)發(fā)生著突變，其中154、156位點(diǎn)和HA分子對(duì)宿主受體的糖特異性有關(guān)[25]，A156T這一突變有助于病毒適應(yīng)陸生家禽宿主，可以提高病毒在這些宿主中的毒力[26]；②氨基酸替換S145L在2002年首次出現(xiàn)，145位點(diǎn)的絲氨酸是a-2，3唾液酸的受體結(jié)合位點(diǎn)，已有的關(guān)于H5N1病毒HA結(jié)構(gòu)的研究表明，這個(gè)位點(diǎn)絲氨酸到亮氨酸的突變有助于a-2，6唾液酸的結(jié)合[27]；③突變K205R發(fā)生在受體結(jié)合區(qū)，是一個(gè)免疫逃逸突變[28]，后被證實(shí)具有很強(qiáng)的抗原性[29]；④D110N在2003年出現(xiàn)在鳥類宿主中，可以增加HA在人類受體中的親和力[30]。

圖2 279、171、498、98、61、100突變位點(diǎn)氨基酸頻率圖表

?

以上幾個(gè)重要的氨基酸突變?cè)贖A結(jié)構(gòu)上的標(biāo)記如圖3所示。

由此，我們認(rèn)為高頻次突變與病毒適應(yīng)性進(jìn)化密切相關(guān)，尋找病毒進(jìn)化過程中的高頻次突變對(duì)研究病毒的自然選擇過程具有重要意義。

3 討論

我們利用生物信息學(xué)手段對(duì)H5N1病毒HA蛋白適應(yīng)性進(jìn)化過程進(jìn)行研究，嘗試對(duì)其進(jìn)化過程中產(chǎn)生的突變做整體把控和系統(tǒng)分析，因此我們建立了一個(gè)快速尋找病毒進(jìn)化過程中重要突變的方法和基于頻率的評(píng)價(jià)體系。

尋找突變的方法核心在于尋找標(biāo)記分枝，然后向根節(jié)點(diǎn)回溯，而標(biāo)記分枝的大小決定了尋找突變的細(xì)化程度。因?yàn)槲覀兊姆椒ɑ谛蛄械幕拘畔⒑拖到y(tǒng)發(fā)生樹，與病毒自身特征無關(guān)，所以該突變尋找方法同樣適用于其他有明顯年份標(biāo)記的病毒。

表3 高頻次突變位點(diǎn)、種類及適應(yīng)性評(píng)價(jià)

圖3 部分高頻次突變?cè)贖5N1病毒HA蛋白結(jié)構(gòu)上的標(biāo)記

為了對(duì)這些突變做出系統(tǒng)評(píng)價(jià)，我們計(jì)算了它們?cè)诟髂攴莸陌被犷l率和固定比例，引入氨基酸頻率圖表來直觀地反映某一位點(diǎn)氨基酸突變和頻率波動(dòng)情況，提示該位點(diǎn)受到的選擇壓力。將氨基酸頻率圖表與固定情況結(jié)合來看，我們提出了頻率變換和有效變換的概念。據(jù)此，我們形成了一個(gè)氨基酸突變的評(píng)價(jià)標(biāo)準(zhǔn)：發(fā)生了有效變換的位點(diǎn)即年際間氨基酸頻率變換的最大值大于0.4，且突變后氨基酸在某一年發(fā)生了固定或基本固定，我們就認(rèn)為它是正選擇位點(diǎn)且有很大可能是抗原表位；如果只滿足頻率變換大于0.4，該位點(diǎn)有一定可能是正選擇位點(diǎn)，且該突變?cè)跇渖系恼急仍叫?，正選擇位點(diǎn)可能性越低?；谶@個(gè)評(píng)價(jià)標(biāo)準(zhǔn)，我們對(duì)照已有研究結(jié)果，發(fā)現(xiàn)吻合程度很高，且我們找的位點(diǎn)更加詳實(shí)準(zhǔn)確，這為實(shí)驗(yàn)人員提供了很好的突變研究思路。

此外，我們發(fā)現(xiàn)許多位點(diǎn)發(fā)生了高頻次突變，利用上面的方法進(jìn)行評(píng)估后，發(fā)現(xiàn)大多數(shù)高頻次突變位點(diǎn)都和正選擇及表位相關(guān)，而且已有研究已經(jīng)揭示了部分位點(diǎn)的功能。因此，我們認(rèn)為發(fā)生高頻次突變的位點(diǎn)受到的自然選擇壓力更大，搜尋病毒進(jìn)化過程中的高頻次突變對(duì)研究病毒的適應(yīng)性進(jìn)化具有重要意義。