H7N9是一种甲型流感病毒，原属于低致病性感冒病毒，仅在禽间发现。然而，在2013年初人类感染H7N9病毒的病例在上海和安徽两地率先发现，首次跨越了种族屏障。H7N9病毒内部的高频率基因突变预警着该病毒进一步引发人群间大流行的风险不容忽视。为更好地适应宿主，病毒的进化速度远远高于其宿主^[1]，了解病毒的进化率是探索病毒进化机制的基础。研究表明感染人类的H7N9是一种新型RNA重组病毒，该病毒外壳是H7N9，内部却是H9N2^[2]。H7N9的8个结构基因面临不同选择压力使其进化速度有所不同，而不同的进化速度将直接影病毒致病性的高低及其对人的亲和力大小。目前对病毒的进化率估算主要基于可变的分子时钟模型和蒙特卡罗马尔科夫链（Monte Carlo Markov Chain，MCMC）模型。可变的分子时钟模型允许病毒在不同时间、不同分支的进化速度具有不同的分布。在贝叶斯理论框架下，设定病毒进化速度的先验分布，采用MCMC模型估算其后验分布^{[3 – 4]}。目前，已有研究通过该方法估算出了多种RNA病毒的进化率^{[5 – 6]}。本研究基于H7N9病毒的全基因组序列，结合H9N2基因片段，估算其8个基因结构区及全基因组的进化速度，为探索H7N9病毒的进化与传播提供理论基础。

本研究选取全球共享禽流感数据倡议组织（The Global Initiative on Sharing Avian Influenza Data，GISAID，http://platform.gisaid.org/epi3/frontend）发布的中国大陆、以人类为宿主的H7N9全基因组序列数据。研究表明新型H7N9病毒除外部结构基因片段血凝素（haemagglutinin，HA）和神经氨酸酶（neuraminidase，NA）来自H7N9外，其余6个内部结构基因片段均来自于H9N2^[2]。为此，本研究还收集了中国大陆禽类宿主H9N2序列数据，剔除不完整的或重复的序列，共保留155条全基因组核苷酸序列。

采用Bio Edit 7.0软件对收集的核苷酸序列进行Clustal W多序列比对，并对H7N9和H9N2序列除NA和HA基因片段外的其他6个基因片段进化合并。比对后，通过MEGA 6.06软件，采用最大似然法构建系统发育树，并运用模型优化指数AICc从软件自带的24种备选碱基替换模型中获得最佳模型。运用bootstrap检验验证系统发育树的可靠性，重复次数为1 000次。

首先利用Path-o-gen软件进行root-to-tip线性回归，初步估算H7N9病毒不同基因结构片段及全基因组进化率，为下一步的贝叶斯估算进化率提供先验信息。使用BEAST 1.8.2软件分析3种分子时钟模型（即lognormal、strict、exponential分子时钟模型）下的进化率。其中，使用的碱基替换模型为上述MEGA分析中获得的最佳模型，链长为100个百万步长。采用Tracer 1.6软件（http://tree.bio.ed.ac.uk）对进化速率计算结果进行分析，运用ACIM值筛选出最佳分子时钟模型；若结果分析中统计量的有效抽样大小（effective sample size，ESS） < 200时，将马尔科夫链长增加到300个百万步长。

进化率数据采用SPSS 20.0软件进行分析和处理，进化率均数比较采用多个独立样本比较的Kruskal-Wallis H检验，合并H9N2病毒数据前后各基因片段是否存在差异采用Mann-Whitney U检验，检验水准为α = 0.05。

据国家疾病与预防控制中心2016年12月发表的H7N9禽流感流行分布数据，至2016年8月31日，共有775例实验室确诊的人感染H7N9禽流感病例，出现4次流行高峰，其中春季为高发季节。在第一波疫情阶段共报道了134例H7N9病例，其中有96 % 的病例集中发生在2013年的4和5月，80 % 的病例集中在江苏、浙江及东南沿海城市。在第二波疫情阶段共有306例H7N9病例和128例死亡，疫情较第一波明显严重，其中发病高峰聚集在2014年1 — 2月，疫情进一步南移至广东。在第三波疫情阶段共有219例H7N9病例和98例死亡，广东和浙江为病例数最多的2个省份，占54 %。虽然相对第二波疫情，第三波疫情病例数有所减少，但死亡率却由42 % 上升至45 %。

本研究共收集155条全基因组序列，含13 103个碱基。结合H9N2序列数据后，共有290条含982碱基的基质蛋白（matrixprotein，MP）基因，338条含1 497个碱基的核蛋白（nucleoprotein，NP）基因，298条含838个碱基的非结构蛋白（nonstructrual，NS）基因，357条含2 151个碱基的聚合酶A（polymerase A，PA）基因，357条含2 273个碱基的聚合酶B1（polymerase B1，PB1）基因，374条含2 280个碱基的聚合酶B2（polymerase B2，PB2）基因。

比对24种备选碱基替换模型下所得的AICc值，得出GRT + G + I为最佳碱基替换模型。基于该模型采用最大似然方法构建H7N9人类宿主病毒序列样本的系统发育树，结果显示，H7N9病毒的HA、NA基因片段及全基因组序列的系统发育树具有相似的结构。提示感染人类的H7N9禽流感病毒同源性较高。第二波疫情及之后的散发疫情较第一次疫情有些病毒进化分支不同，同一时间同一地点进化关系较强。

图 1

图 1 H7N9病毒的HA片段系统发育树

图 2

图 2 H7N9病毒的NA片段系统发育树

图 3

注：图中未标注代表第一波H7N9禽流感疫情（2013年2月1日—6月30日）；空心正方形代表第二波H7N9禽流感疫情（2013年7月1日 — 2014年7月31日）；实心正方形第二波疫情之后的散发病例。 图 3 H7N9病毒的全基因组系统发育树

运用root-to-tip线性回归方法初步估算H7N9禽流感病毒的进化率，结果显示，H7N9禽流感病毒不同基因结构片段的进化率为10^–2～10^–5次/位点/年，其中主要集中在10^–3次/位点/年水平，合并H9N2序列数据后估算的进化速度较未合并前小。

表 1

表 1 H7N9病毒线性回归进化率（× 10–3次/位点/年） 合并H9N2 全基因组 HA NA MP NP NS PA PB1 PB2 否 7.97 5.92 4.92 3.15 4.63 2.97 2.60 4.07 4.91 是 0.43 3.47 2.63 1.51 1.49 0.32

表 1 H7N9病毒线性回归进化率（× 10–3次/位点/年）

H7N9在最佳分子时钟模型下的进化率估算结果，其中最佳模型采用AICM指标筛选。在未合并H9N2序列分析中，全基因组序列进化速度为 4.60 × 10^–3次/位点/年；8个基因结构区中PA、HA和NA的进化速度较快，PB1、PB2次之，NS、NP和MP进化率较低。在合并H9N2序列估算的进化率中，PB1、PB2和PA进化率较高。根据进化率均数可看出大多数合并后估算的H7N9 进化速度会相对慢些，但总体来说H7N9 病毒除HA和NA外的 6 个基因结构片段进化率普遍在3 × 10^–3～5 × 10^–3次/位点/年。其中，合并H9N2 序列数据前后的PA基因片段的进化率差距最大，但PA基因片段合并H9N2 数据前的抽样误差（5.36 × 10^–5次/位点/年）大于合并后的抽样误差（2.13 × 10^–5次/位点/年）。Kruskal-Wallis H检验表明无论是否合并H9N2碱基序列，基因结构区之间的进化率差异有统计学意义（P < 0.001）。进一步两两比较，发现结构基因片段两两之间差异有统计学意义（ P < 0.001）。提示H7N9禽流感病毒不同基因结构区的进化速度不同。

表 2

表 2 H7N9分子时钟最佳模型估算进化率（× 10– 3次/位点/年） 基因片段 合并H9N2 $\bar x \pm s$ 95 % CI 最优模型 Pa值 全基因组 4.60 ± 0.38 3.94～5.40 严格时钟 HA 7.43 ± 1.49 4.72～10.4 指数时钟 NA 5.92 ± 0.74 4.52～7.40 严格时钟 MP 否 3.02 ± 0.47 2.07～3.89 严格时钟 < 0.001 是 3.19 ± 0.34 2.57～3.87 指数时钟 NP 否 4.34 ± 0.50 3.38～5.31 严格时钟 NS 否 4.45 ± 0.83 2.91～6.09 严格时钟 < 0.001 是 3.47 ± 0.35 2.81～4.15 指数时钟 PA 否 8.03 ± 1.69 4.85～1.13 指数时钟 < 0.001 是 4.02 ± 0.37 3.36～4.78 指数时钟 PB1 否 5.02 ± 0.96 3.36～7.02 指数时钟 < 0.001 是 4.79 ± 0.40 4.03～5.55 指数时钟 PB2 否 5.60 ± 1.07 3.72～7.85 指数时钟 < 0.001 是 4.25 ± 0.37 3.57～5.01 指数时钟 　　注：a Mann-Whitney U检验。

表 2 H7N9分子时钟最佳模型估算进化率（× 10– 3次/位点/年）

尽管已有研究通过人类及禽类宿主的H7N9禽流感碱基序列探索H7N9病毒的进化率^{[7 – 8]}，但本研究考虑到新型重组的H7N9病毒不是由经典的H7N9病毒进化而来^[9]，故首先独立探索人类宿主H7N9病毒进化率，继而联合禽类H9N2病毒的6个基因结构片段进一步估算H7N9病毒进化率。通过收集来自不同地点、不同时间的人类宿主的H7N9禽流感病毒，构建全基因组及各基因结构片段的系统发育树，表明感染人类的H7N9禽流感毒具有高度同源性。此外，与Lam等^[10]研究发现HA基因片段在第二波疫情中具有更多的进化分支结论类似，本研究发现第二波疫情及之后的散发疫情的H7N9禽流感病毒相对于第一波疫情病毒有所进化，表明禽流感病毒为适应宿主不断地进化^[11]，以逐渐适应宿主的免疫压力^[12]。

经由分子时钟模型估算的H7N9病毒进化率结果显示，H7N9禽流感病毒的进化速度较快，不同基因片段的进化速度不同，这与其所承担的选择压力不同有关。全基因进化率为4.60 × 10^–3次/位点/年，HA和NA基因的进化率分别为7.43 × 10^–3和5.92 × 10^–3次/位点/年，分别比H5N1禽流感病毒HA和甲型N9亚型流感病毒进化快^{[13 – 14]}。表明HA和NA基因片段特别是HA可能具有更高的选择压力。HA与NA作为H7N9禽流感病毒表面主要的2种糖蛋白，其抗原变异性大^[2]，可能是导致其进化速度快于其他基因结构片段的主要原因。同时，有研究表明HA在氨基酸位点的进化使其与人体的亲合力增加^{[15 – 16]}，NA分子水平的进化使H7N9禽流感病毒产生神经酰胺酶抑制剂耐药性的出现^[17]；这2种糖蛋白的快速进化很可能是H7N9禽流感流行的主要原因之一。除HA和NA基因片段外，H7N9病毒的其余6个基因结构片段进化率普遍在3 × 10^–3～5 × 10^{– 3}次/位点/年，PB1、PB2、PA基因片段进化速度相对较快，这与各基因片段发挥的不同结构功能以及其所受的选择压力有关。如有研究表明PB1与H7N9的高致病性有关^{[18 – 19]}，PB1基因片段的进化促进了H7N9病毒从禽类宿主低致病性到人类宿主的高致病性的转变。此外，由于目前人类宿主的H7N9病毒碱基序列较少、年度跨越较短，故合并H9N2病毒序列估算的进化率更可能反映H7N9病毒的进化速度。

国内禽流感病毒从最初H5N1甲型禽流感的发生，到H9N2禽流感病毒的暴发，再到如今新型H7N9禽流感病毒的流行，均暗示着对活禽市场与人群缺乏有效流感病毒系统监测手段。H7N9禽流感病毒是否会在未来进一步在人类中发生大规模流行，仍需提高警惕。