流行性感冒(简称流感)是由流感病毒引起的急性呼吸道传染病，历史上曾发生多次世界范围的大流行。根据抗原性不同，人类流感病毒分成甲、乙、丙三型。甲型流感病毒基因易变且易引起流行，乙型流感也常引起局部暴发[1]。近年乙型流感受到全球关注，有研究报道，年龄≥50岁人群乙型流感感染者中，临床病死率可达到2.5% (95%CI，0.7%~7.6%)[2]。

乙型流感病毒基因组具有单负链8个节段，其中血凝素(hemagglutinin，HA)基因位于第4节段，神经氨酸酶(neuraminidase, NA)基因位于第6节段[1]。HA在病毒包膜上以三聚体形式存在，识别并结合宿主细胞表面的受体结合位点(receptor binding sites，RBS)，具有特异性[1]；同时刺激宿主产生保护性中和抗体，是设计流感疫苗的常用靶点[1]。NA由4条糖蛋白链组成，单体头部是NA活性中心，具有型别特异性；此外NA持有神经氨酸酶活性(涉及感染细胞内成熟病毒释放)和抗原性(抗原特异性)。乙型流感8个节段中，以HA和NA基因的生物学功能最为重要，是流感疫情流行的分子基础。

国内山东济宁地区2018-2020年乙型流感Bv毒株流行，属于Bv.1a株系[3]；2021年下半年，广东地区流感疫情频繁，较去年同期每月疫情上升50%~85%。结合本地区的监测工作，初步鉴定为乙型流感活动[4]。为深入了解广东地区2020-2021年乙型流感HA和NA基因的变异和进化特征，本研究开展毒株HA/NA基因分析，以期为乙型流感防治提供分子水平依据。

广东各地区流感监测网络采集流感咽拭并分离流感病毒；时空分布方法筛选流感毒株，其中2020年和2021年分别为8株和26株，合计34株。从GenBank下载全球人类乙型流感毒株的HA/NA基因序列共23对，包括世界卫生组织(World Health Organization，WHO)疫苗株4株[5]，株系代表株5株(2008-2021年)；此外亚洲7株(含中国3株)，欧洲2株和北美5株。

根据2015-2020年全球乙型流感Bv毒株HA和NA基因核苷酸序列，采用Primer Premier 5.0软件(PREMIER Biosoft International，Palo Alto，CA，USA)设计和兼并，合成HA和NA基因扩增引物各3对。

采用QIAGEN公司QIAamp Viral RNA mini Kit(Lot 169018932)提取RNA，置-20 ℃保存；采用Invitrogen的SuperScriptⓇ Ⅲ One-Step RT-PCR System with PlatinumⓇ Taq (Lot 1996333)进行RT-PCR扩增。扩增产物鉴定，Size Marker为100~2 500 bp。

PCR产物纯化采用广州天一辉远基因科技有限公司磁珠法琼脂糖凝胶PCR回收试剂盒，测序用ABI BigDye Terminator v3.1 Cycle Sequencing Kit，在ABI 3730 DNA Analyzer上测序。根据顺义和反义引物的测序结果，采用Lasergene v8.1.3软件(DNASTAR，Madison，WI，USA)拼接核苷酸(Nucleotide，Nt)，验证检测序列的正确性。

采用MEGA 7.026软件的Clustal W进行核苷酸同源性分析，用邻接法(Neighbor-joining)构建HA/NA基因进化树。通过BioEdit v7.2.6(Ibis Biosciences，Carlsbad，CA，USA)计算HA/NA上每个氨基酸(Amino acid，AA)位点的熵值(Entropy value)，并进行评价。通过NetNGlyc 1.0 Server对氨基酸序列进行糖基化位点分析。通过Datamonkey 2.0，分别采用SLAC (Single likelihood ancestor counting)、FEL (Fixed-effects likelihood)和FUBAR (Fast，Unconstrained Bayesian App Roximation)方法，对HA/NA基因进行进化分析。

规整基因序列资料，采用SPSS 23.0对核苷酸和氨基酸变异信息进行独立样本t检验。以P < 0.05为差异有统计学意义。

采用时空分布法筛选2020年(8株)和2021年(26株)毒株(2020年4-12月监测无Bv分离株)，包括潮州5株、东莞1株、广州2株、河源1株、惠州3株、江门2株、茂名2株、梅州1株、清远2株、汕头1株、汕尾2株、韶关5株、深圳3株、湛江1株、中山1株和珠海2株。

广东Bv毒株HA和NA基因片段，分别为1 875 bp和1 556 bp；编码区(ORF，Open reading frame)分别为1 749 bp和1 401 bp。基因比对发现：①疫苗株B/Colorado/06/2017的HA基因ORF丢失Nt529-534序列(aaaaac)；此后2017年以后的毒株[包括2019和2021年疫苗株(2株)、广东毒株(34株)和全球毒株(14株)]均丢失Nt529-537序列(aaaaacgac)；②以疫苗株B/Brisbane/60/2008的NA基因为基准，其ORF共有1 401 bp，在Nt223后插入tcc，出现在B/Bangladesh/8007/2021和B/Singapore/WUH7588/2021毒株。

以4株疫苗株基因序列作为基准，计算2020年和2021年广东毒株的核苷酸同源性，结果如表 1所示。同源性分析揭示：①4株疫苗株中，广东2021年流行毒株HA基因同源性与疫苗株B/Austria/1359417/2021同源性最高[(99.19±0.26)%]，而NA基因同源性却与2019年疫苗株B/Washington/02/2019同源性却最高[(99.23±0.11)%]；②以4株疫苗株HA核苷酸序列为基准，2021与2020年毒株HA基因均有统计学差异(P < 0.01)，提示广东2021年毒株HA基因发生极其明显变异；③以疫苗株B/Washington/02/2019的NA核苷酸序列为基准，2021年毒株与2020年毒株NA基因同源性有统计学差异(P < 0.01)，提示广东2020年毒株NA基因有明显变异；④疫苗株B/Austria/1359417/2021的HA基因与2021年毒株HA基因同源性为(99.19±0.26)%；而与2021年下半年毒株(17株)HA基因同源性高达(99.28±0.25)%，提示广东2021年下半年流行毒株HA基因与当前WHO推荐疫苗株基因依然高度契合。

广东2020-2021年的34株毒株HA基因进化特征，如图 1A所示，结果显示：①广东2020年毒株(8株)与广东2021年毒株(下半年分离，17株)具有不同基因特征：2020年毒株呈Bv.1a.3株系特征，2021年毒株主要呈Bv.1a.3a.2株系特征。；②广东2020年毒株与2021年毒株(上半年，9株)具有疫苗株B/Washington/02/2019基因特征；国外2020年毒株包括欧洲(圣彼得堡)、美洲(加州和宾州)、亚洲(新加坡、爱知)和2021年美洲毒株(德州和新墨州)；③广东2021年毒株(下半年分离)呈现疫苗株B/Austria/1359417/2021基因特征，具体表现国外是B/Paris/9867/2020最早出现(2020-11-25)，随后是国内广东B/Guangdong/073/2021 (2021-01-02)，北京B/Beijing/Chaoyang/12099/2021 (2021-08-31)；另包括国外(密歇根、孟加拉和新加坡)2021年1月后的分离毒株。

广东毒株NA基因进化如图 1B所示，主要特征包括：①广东2020年毒株(8株)与2021年(26株)具有不同基因特征，广东2020年毒株NA基因与同期圣彼得堡、加州、新加坡、宾州、爱知、中国香港的毒株同源性类似；②广东2021年上半年毒株(1~6月，9株)与下半年部分毒株(17株)也具有一定差异；③广东2021年26株NA基因相似度更接近疫苗株B/Washington/02/2019的基因(99.23±0.11)%。

HA基因ORF的1 749 bp编码582个氨基酸(AA)；以疫苗株B/Washington/02/2019的HA基因为基准，50株中，31个氨基酸位点发生置换，其中广东34株毒株有26个置换，全球14株(2020-2021年)有16个置换，疫苗株B/Austria/1359417/2021的HA基因有8个置换。以疫苗株B/Austria/1359417/2021的HA基因为基准，广东和全球毒株置换差异一致。

以疫苗株B/Washington/02/2019的HA基因为基准，置换位点包括位点H137Q/N(广东2021，16/26；全球2021，1/6)、位点S209D/N(广东2020，7/8和广东2021，26/26；全球2020，8/8和全球2021，6/6)、位点A142T、P159L和K215R(广东2021，15/26；全球2020，1/8和全球2021，4/6)、位点R148G、G196E和R291K(广东2021，26/26；全球2020，2/8和全球2021，4/6)、位点V232M和P253Q(广东2021，10/26；全球2020，1/8和全球2021，4/6)。

疫苗株B/Austria/1359417/2021的HA基因特征在于，位点A142T、R148G、P159L、G196E、S209D、K215R和R291K已经发生置换；即在这些位点上，疫苗株B/Austria/1359417/2021的HA氨基酸序列已体现出广东2021毒株特征。比较广东2021年下半年毒株与疫苗株B/Austria/1359417/2021的HA基因差异，在于K071E(2/19)以及T182I、PL184-185SS、A214T和G519R(1/19)置换。

疫苗株B/Washington/02/2019的NA基因ORF具有1 401 bp，编码466个AA；以此NA基因为基准，AA72后插入亮氨酸(L)见于B/Bangladesh/8007/2021和B/Singapore/WUH7588 /2021毒株。50株中，共29个AA位点变异，其中广东(34株)和全球(14株)分别有14和20变异位点；此外，疫苗株B/Austria/1359417/2021的6个位点变异。

以2019年疫苗株NA基因为基准，变异特点包括：①广东2020年毒株仅E313G变异；②2021年毒株包括D53N(13/26)、N59S(13/26)、G233E/R(14/26)、V303I(13/26)和K343E(15/26)；以及D53N、N59S、G233E、V303I和K343E；③少数基因变异还包括F12L、I45M、P51Q、P119T、V249I、I263V、G379E和E406D。

疫苗株B/Austria/1359417/2021的HA蛋白共有潜在糖基化位点10个，即NVT40-42、NCT74-76、NIT160-162、NKT178-180、NQT245-247、NKS316-318、NGT345-347、NIT530-532、NHT543-545和NVS575-577。此次变异中，HA氨基酸D209N变异，导致增加糖基化位点NET209-211，见于B/Guangdong/428/2020、B/Guangdong/523/2020和B/Guangdong/053/2021；PL184-185SS变异，增加NSS183-185糖基化位点，见于B/Guangdong/116/2021。疫苗株B/Austria/1359417/2021的NA蛋白共有潜在糖基化位点4个，NAS56-58、NRS64-66、NGT144-146和NKT284-286；未发现置换变异。

以疫苗株B/Brisbane/60/2008的NA为基准，已报道NA单耐药位点共有29个[7]和双耐药位点4个；在本研究中未出现上述具体氨基酸位点变异。

根据50株毒株HA氨基酸熵值结果，熵值≥0.600包括位点148、165和291(均为0.641)，142、159和215(均为0.680)，137(0.717)，196(0.731)和209(0.779)；NA氨基酸包括位点344(0.680)和234(0.684)。见表 2。

50株HA和NA基因，进行选择性分析(NA未计入插入L072，见B/Bangladesh/8007/2021和B/Singapore/WUH7588/2021)。SLAC和FEL两方法中кHA＜кNA，而FUBAR方法中кHA＞кNA，提示不同计算方法导致HA基因的к值[转换(α)/颠换(β)]不同。SLAC检验ω=dN/dS比值，如果ω＜1，且P(ω＜1)＜0.1则存在负向选择位点，SLAC检验分别位于HA基因73，138，437和497位点和NA基因179，321和335位点。

50株毒株采用FEL检验，选择性分析显示，HA和NA基因分别筛选出28个和14个负向选择位点；采用FUBAR检验，分别筛选出HA基因15个负向选择位点和1个正向选择位点(AA209)，NA基因10个负向选择位点和1个正向选择位点(AA384)。根据多种方法结果一致性的原则，综合评价，HA和NA基因氨基酸负向选择位点分别是73，138，437和497位点，及179，321和335位点。

自2009年流感H1N1pdm出现以后[8]，广东流感各型别(包括H3N2和B型)共同存在流行季节中，交替成为优势型别。本研究中，①2020-2021年流感优势毒株被鉴定为B型Victoria株系；②2020年4-12月监测未分离到流感各亚型毒株，这可能与流感流行特征以及国内新冠疫情的严格防控措施有关系；③2021年下半年Bv毒株多从局部疫情暴发人群中分离。

分别以4株疫苗株HA基因为基准比较广东2020与2021年毒株之间差异，发现4株疫苗株的基准的2年间HA基因同源性均有统计学差异(P < 0.01)，NA基因在2020年(与2020年前比较)也发生变异，但仅与B/Washington/02/2019基准差异明显，提示广东Bv流行株HA基因在2021年(与2020年比较)发生显著性变异。而且广东Bv流行株HA基因变异特征，主要呈现于2021年下半年的疫情毒株。杭州2014-2019年流行Bv的Bv.1a株系(类似疫苗株Bv/Colorado/06/2017)，HA和NA基因同源性分别为94.67%和97.13% [9]。

此次广东毒株变异(以疫苗株B/Washington/02/ 2019的HA基因为基准)见于A表位[V132I (1/26，202 1)、H137N/Q (16/26，2021)、A142T(15/26，2021)；R148G(26/26，2021)]、B表位[G156R(1/26，2021)、P159L (15/26，2021)；N162K(26/26，2021)]；D表位[S209N(7/8，2020)、S209D(25/26，2021)；T211A(5/8，2020)]。与此比较，杭州Bv.1a株系毒株变异见于A表位(L125I、H131N、I132V)、B表位(N141D/G、G148R)和D表位(I195V)[9]；云南地区2017-2018年流行乙型流感By/Bv毒株，其中Bv毒株HA基因具有V161I突变特征，同时S14N和V195I发生变异[10]。

与2019年和2021年疫苗株HA基因比较，广东2020年毒株呈Bv.1a.3株系特征，2021年毒株主要呈Bv.1a.3a.2株系特征[5]，但变异幅度更大。WHO已于2021年9月更换疫苗组Bv毒株为B/Austria/1359417 /2021，与该疫苗株比较，2021年广东株在表位和RBS区域变异包括H137Q/N、R148G和P253Q。以疫苗株B/Washington/02/2019的NA基因为基准，2021年广东毒株NA基因5个氨基酸(D53N、N59S、G233E/R、V303I和K343E)发生明显变异(50.0%~57.7%)。此变异幅度既超越疫苗株B/Austria/135941 7/2021 (2021-10 -09，仅K343E)，也超越2021年上半年的北京毒株B/Beijing/Chaoyang/12099/2021(2021-08-31，仅D53 N、N59S和G233E)。

10个HA和4个NA基因潜在糖基化位点，经分析均无明显变异趋势。NA基因已知药物敏感位点未发现变异。基因氨基酸位点的熵值反映该位点变异频率[6]。本研究中，50株毒株的HA和NA分别存在9个(148/165/291/142/159/215/137/196/209)和2个(344/234)易突变位点。

由于流感病毒缺乏RNA聚合酶校对功能，HA和NA基因易发生点变异[11]，常使用错义替代与同义替代的比值(dN/dS)评估氨基酸位点所受的自然选择压力。当dN/dS＞1时，表示该位点为正向选择位点；dN/dS=1时，该位点为中性进化不受选择；dN/dS＜1时，则表示该位点为负向选择位点[6]。本研究中，50株Bv经SLAC、FEL和FUBAR 3种方法均筛选出HA和NA基因存在负向选择位点，分别位于73，138，437和497位点，及179，321和335位点。此次筛选出的负向选择位点，提示乙型流感病毒Bv毒株HA和NA基因在广东地区的自主变异进化明显，已超越免疫压力因素的作用。

在人类流感病原体中，甲型流感病毒(H1N1、H3N2)的HA/NA基因表现变异频繁，通常存在基因正向选择位点[12]。此次广东Bv毒株出现多个氨基酸的负向选择位点，在流行病学和基因进化中的意义有待于进一步探讨。然而，结合流感疫情现状，同时关联核苷酸/氨基酸的同源性比较、抗原表位的变异、负向选择位点筛选，可以推断广东乙型流感Bv毒株具有引发流感局部/洲际疫情的流行病学意义。总之，本研究提示，①Bv毒株分子(HA等)变异是导致广东2020-2021年流感局部暴发的原因；②2021年下半年的流感进一步变异，可能引起2022年Bv毒株局部流行；③尽管WHO的2021年9月推荐的疫苗株依然有效，但广东2021年Bv毒株的变异，提示WHO应考虑更新疫苗株的可能性。