ABSTRACT
季节性甲型流感病毒进化迅速，然而针对在中国广泛使用的聚合酶酸性蛋白（PA）抑制剂巴洛沙韦马波西酯（BXM）的耐药性出现及其分子基础仍不明确。为解答此问题，研究从国家监测网络和GISAID数据库中收集了2018年至2025年间来自中国内地流感患者的3,938条PA基因序列。通过筛选上市后出现的、至少在两个样本中出现且频率低于50%的PA蛋白N端结构域（PAN）突变，鉴定出25个单点突变以及6个可能与药物压力相关的连锁突变。随后评估了这些突变对BXM敏感性的影响。分析揭示了已知的与BXM敏感性降低相关突变（如L28P、K34R、E198K）的出现，但其流行率仍然很低（2/3850, 0.05%）。值得注意的是，研究鉴定出一个新型替换突变D27G，该突变导致病毒对BXM的敏感性比野生型病毒降低约12.4倍，并且在人气道类器官中显示出比经典耐药突变I38T更高的复制适应性。分子动力学模拟进一步表明PA-D27G减弱了PA与BXM活性形式巴洛沙韦酸（BXA）之间的相互作用。流行病学分析显示D27G突变在中国内地仍很罕见，在四个分离株中检测到（4/1247, 0.32%），在全球分离株中呈零星流行（<0.1%, 9/53,132）。总之，这些结果证明了中国BXM相关耐药性的早期出现，并确定PA-D27G是一个与耐药性相关且保持病毒适应性的突变，强调了持续基因组和流行病学监测的重要性。
Introduction
流感是一种高度传染性的呼吸道病原体，对全球公共卫生构成重大威胁。据世界卫生组织（WHO）统计，流感每年影响超过10亿人，导致300万至500万例严重感染和29万至65万例死亡。由于对奥司他韦等常规抗病毒药物的耐药性日益增加以及普通人群疫苗接种覆盖率低，严重流感肺炎的管理仍然面临重大挑战。作为RNA病毒，流感具有高突变率，其进化速率通常超过大多数RNA病毒，尽管该速率因基因节段、宿主背景和病毒株系等因素而异。具体而言，表面糖蛋白血凝素（HA）和神经氨酸酶（NA）在持续的抗原漂变驱动下进化更快，而内部基因（如聚合酶酸性蛋白PA）的进化速率较低。此外，不同的病毒谱系表现出不同的进化速率，A/H3N2由于更强的抗原漂变和更频繁的选择性清除，其进化速度快于A/H1N1pdm09和乙型流感。宿主背景进一步调节进化动态，因为循环于人类、猪和禽类物种中的病毒经历不同的免疫和细胞环境，这塑造了适应性突变模式。此外，药物诱导的选择压力可导致特定突变，例如神经氨酸酶中的H275Y突变，该突变赋予病毒对奥司他韦的显著耐药性。
巴洛沙韦马波西酯（BXM）是一种新型强效抗流感病毒药物，特异性靶向甲型流感病毒聚合酶酸性蛋白（PA）的N端内切酶结构域（CEN，约第1至200位氨基酸）。通过抑制“帽状结构抢夺”机制，进而阻断病毒RNA合成，BXM发挥其抗病毒作用。与常规抗病毒药物相比，BXM具有若干优势，包括单次口服、起效迅速、广谱疗效以及较低的耐药性发展风险。然而，由于流感病毒的高突变率，尤其是在抗病毒药物施加的选择压力下，耐药突变可能迅速出现。其中，I38T突变是与显著耐药性相关的标志性突变，能将病毒株的耐药性提高10至100倍。2024年7月，WHO发布了跨流感亚型的已记录BXM耐药突变综合列表，其中大多数是通过实验室条件下的耐药性选择试验鉴定出来的。因此，研究循环流感毒株中BXM耐药性的真实世界流行情况至关重要。
中国长期面临流感病毒的持续循环和季节性爆发。在COVID-19大流行期间，社交距离等公共卫生干预措施显著降低了流感活动。然而，随着2022年底COVID-19限制的放松，流感爆发迅速恢复到疫情前水平。BXM于2021年4月27日正式在中国内地获批使用，并在2022/2023流感季节（后COVID时代）开始广泛临床应用。自获批以来，BXM已连续用于三个流感季节。然而，尚不确定在这三个流感季节期间，中国是否出现了BXM耐药株，或者是否出现了与BXM相关的新型耐药突变。
本研究通过国家传染病医学中心的微生物研究网络和GISAID数据库，收集了2018年4月至2025年1月期间的H3N2和H1N1pdm09流感病毒PA基因序列，共计3,938个样本。在2021年4月27日BXM在中国获批后，鉴定出25个潜在耐药单突变和6个连锁突变。观察到D27G突变导致病毒对BXM的敏感性比野生型毒株降低约12.4倍，并在人气道类器官模型中得到证实，其显示出比I38T毒株更高的复制适应性。利用分子动力学模拟研究了D27G突变对PA蛋白与巴洛沙韦酸（BXA）相互作用的潜在结构影响。此外，评估了携带D27G突变分离株的全球流行情况。这些发现突显了与BXM敏感性降低相关突变的出现，鉴定了新型D27G突变，并强调了持续进行BXM耐药性基因组和流行病学监测的重要性。
Materials and Methods
数据收集
本研究中分析的中国H3N2和H1N1pdm09样本来自国家传染病医学中心的全国微生物研究网络和GISAID公共数据库。该网络获得了复旦大学附属华山医院伦理委员会的批准，并获得了所有参与者的知情同意。
从GISAID公共数据库检索了全球H1N1pdm09聚合酶酸性（PA）基因序列，用于PA D27G突变的流行率和系统发育分析，涵盖中国内地（2018年4月27日至2025年1月31日）、美国（2015年10月24日至2025年1月31日）、加拿大（2017年2月19日至2025年1月31日）、香港（2016年2月22日至2025年1月31日）、欧盟（EU）（2018年1月7日至2025年1月31日）、澳大利亚（2017年2月21日至2025年1月31日）和日本（2015年3月14日至2025年1月31日）的序列。收集了每个国家/地区引入BXM前三年的PA序列，直至2025年1月31日。
全基因组测序
对国家传染病医学中心全国微生物研究网络收集的临床样本进行测序。逆转录和PCR扩增（RT-PCR）：选择单个检测中CT值<30且H1N1pdm09和H3N2呈阳性的样本，每两周进行一次逆转录和PCR扩增。使用Takara PrimeScript™ RT试剂盒（Perfect Real Time）进行逆转录。PCR扩增采用合作研发的甲型流感病毒全基因组扩增检测试剂盒。进行两轮多重PCR反应，然后进行两轮纯化以获得文库。合格文库的总量需超过50 ng，目标片段大小范围为300 bp至400 bp。随后使用Illumina NovaSeq 6000平台进行文库制备和测序。
基因组分析
将临床测序FASTQ文件转换为FASTA格式进行下游分析。使用Trim-galore修剪双端文件以去除低质量区域，碱基质量阈值（-q）设为20。使用Bowtie2将修剪后的文件与参考序列进行比对。使用Samtools将SAM文件转换为BAM格式，然后使用RSeQC进行质量控制评估。使用BCFtools生成VCF和一致性序列FASTA文件。比对率低于80%的样本被排除在进一步分析之外。
应用Nextclade Web/CLI对GISAID和临床FASTA样本进行质量控制、比对和氨基酸翻译。用于比对的参考基因组与SNP calling过程中使用的参考基因组相同。为尽量减少不相关突变的选择，选择最近的参考基因组而非WSN/HK68。在Nextclade QC报告中被标记为“bad”或“mediocre”的序列被排除，同时排除包含重复名称、模糊氨基酸（>0.5%）、覆盖率<95%或核苷酸突变超过Q3 + 3IQR的序列。
使用香农熵衡量每个位置的氨基酸多样性。
系统发育树
使用MAFFT v7.525进行多序列比对，使用--auto参数选择最优策略。使用IQ-TREE v2.4.0构建最大似然树，使用-m TEST选项进行模型选择，并进行1000次超快速自举和近似似然比检验（aLRT）重复。使用TreeTime构建时间系统发育树，快速缩放分支，估计进化速率，并从根到尖回归中排除异常值。随后，使用BEAST v10.5.0构建时间标度系统发育树，使用HKY模型和GTR + G4替换模型。先验进化速率基于文献数据和TreeTime结果，H1N1pdm09设为0.0025，H3N2设为0.002。MCMC链长设为1亿次迭代，进行3-4次独立运行。使用LogCombiner合并log和tree文件，burn-in为10%，并确认有效样本大小（ESS）>200。使用TreeAnnotator生成最大后验树。
为减少计算负荷并确保代表国家/地区的序列数量平衡，使用SeqKit去除具有相同序列的重复样本，但保留携带D27G的毒株。对于全球分析，在比对后对序列进行下采样。质量控制后，每年从每个国家/地区随机选择250条序列。如果样本量<250，则纳入该年所有序列用于系统发育树构建。同样应用BEAST构建全球H1N1pdm09时间标度系统发育树。
细胞和试剂
HEK293T和MDCK细胞在补充有青霉素（100 U/mL）、链霉素（10 mg/mL）和10%胎牛血清的DMEM中培养。细胞在37°C、5% CO2的加湿培养箱中维持，达到汇合度时以1:5的比例传代。BXM前药的游离酸活性形式，即巴洛沙韦酸，购自MedChem Express。
质粒和抗体
用于重建A/Hong Kong/68 (H3N2)微型复制子系统的pcDNA质粒和用于A/WSN/1933 (H1N1)反向遗传学的pHW2000质粒由相关教授惠赠。兔抗甲型流感病毒核蛋白抗体、兔抗甲型流感病毒H1N1 HA抗体、兔抗β-肌动蛋白抗体和山羊抗兔IgG（HRP）抗体购自GeneTex。小鼠抗TUBULIN抗体购自Santa Cruz Biotechnology，而山羊抗兔抗体和山羊抗小鼠抗体购自Abcam。
反向遗传学
使用pHW2000八质粒系统生成重组流感病毒（流感A/WSN/1933及其突变变体）。将HEK293T细胞和MDCK细胞与以下质粒共转染。转染4小时后，更换为补充有0.5%胎牛血清（FBS）的DMEM。在转染后48小时收获病毒，并储存于-80°C直至使用。
空斑试验
用病毒感染MDCK细胞单层。孵育1小时后，用PBS洗涤细胞，随后覆盖含有1%低熔点琼脂糖和TPCK处理的胰蛋白酶（0.5 μg/mL）的DMEM，然后在37°C孵育2-3天。用4%多聚甲醛在室温下固定细胞15分钟。去除琼脂糖覆盖层后，用1%结晶紫甲醇溶液染色10分钟。计数每个孔的空斑数。
EC50耐药性测试
使用CellTiter-Glo®细胞活力测定法测定BXA的半数最大效应浓度（EC50）。将MDCK细胞以每孔1×104个细胞的密度接种到96孔板中。每孔接种50 μL病毒悬液（感染复数MOI=0.01），然后加入50 μL系列稀释的BXA。将板在37°C孵育72小时。按照制造商说明进行细胞活力测定。所有药物浓度至少进行三次独立重复实验。使用GraphPad Prism软件（版本8.0）计算EC50值。
微型复制子试验
在转染前一天将HEK293T细胞接种在6孔板中。以1:1:1:1:1的比例转染pcDNA-NP、pcDNA-PA、pcDNA-PB1、pcDNA-PB2和pHW2000-HAmt质粒。转染4小时后，更换为含有系列稀释BXA的完全DMEM。在转染后24小时收获细胞，并在冰上用RIPA裂解缓冲液裂解30分钟。通过12,000 rpm离心15分钟去除细胞碎片。裂解物进行SDS-PAGE并转移至PVDF膜。使用兔抗HA抗体进行免疫印迹，然后与相应的二抗孵育。使用ChemiScope S7系统显示蛋白条带。
病毒生长动力学
在MDCK细胞中分析流感病毒及其突变株的生长动力学。细胞培养和感染条件与空斑试验一致，MOI为0.001。在感染后2至72小时的时间点收集上清液，并通过空斑试验测定病毒滴度。
病毒产量抑制试验
接种MDCK细胞的6孔板用流感病毒（MOI=0.001）感染，并在补充有0.5% FBS、0.5 μg/mL TPCK处理的胰蛋白酶和1 nM BXA的DMEM中培养。在感染后2至72小时的时间点收集上清液，并通过空斑试验测定病毒滴度。
气道类器官的建立、扩增和分化
本研究遵循《赫尔辛基宣言》进行，并获得复旦大学附属华山医院机构审查委员会的批准。在组织收集前获得了所有参与者的书面知情同意。从因非恶性或局限性肺部疾病接受手术切除的患者身上收集肺组织。用于建立类器官的样本取自病变邻近的视觉正常区域。所有组织在切除后立即处理以优化细胞活力。将新鲜肺组织用无菌手术刀切成小碎片，并转移到含有PBS的离心管中。通过300 g离心5分钟使组织碎片沉淀，随后在37°C用胶原酶消化30分钟。将消化后的组织悬浮液用巴斯德吸管机械分离，并通过100 µm细胞筛过滤。加入2毫升胎牛血清（FBS）终止酶反应，通过300 g离心5分钟收集细胞。将细胞沉淀重悬于红细胞裂解缓冲液中，在室温下孵育5分钟，然后用基础培养基中和，并在4°C、300 g下离心5分钟。用基础培养基再洗涤一次后，将最终细胞沉淀重悬于75% Matrigel中，并在预热的12孔悬浮培养板中铺板为40 µL液滴。凝胶化后，液滴在扩增培养基中于37°C、5% CO2的加湿培养箱中孵育。每隔一天更换一次培养基。每1-2周以1:3至1:10的比例传代未分化的类器官，具体取决于使用机械剪切还是酶消化来分割类器官。为了诱导分化为气道类器官，将肺类器官最初在扩增培养基中培养7天，随后在近端分化培养基中再培养7天。每隔一天更换一次培养基。
人气道类器官的感染
用病毒（MOI约为0.1-1）感染人气道类器官，在37°C下孵育2小时。在感染后2、24和48小时收集上清液，通过逆转录定量聚合酶链反应（RT-qPCR）和空斑试验测定病毒滴度。在感染后10和48小时用4%多聚甲醛固定类器官，用于后续免疫荧光染色和流式细胞术分析。
RNA提取和RT-qPCR
为了评估人气道类器官中的病毒生长动力学和宿主免疫反应，在指定时间点收集类器官和上清液。使用MiniBEST Universal RNA Extraction Kit和MiniBEST Viral RNA/DNA Extraction Kit分别从类器官和上清液中提取RNA。使用PrimeScript™ RT reagent Kit（Perfect Real Time）和uni-12及Oligo-dT引物进行cDNA合成。使用TB Green® Premix Ex Taq™ (Tli RNaseH Plus)和基因特异性引物进行实时定量PCR（qPCR），以确定宿主和病毒基因表达水平，这些水平均归一化至GAPDH的表达水平。
免疫荧光染色
为了鉴定病毒感染的细胞，对感染或模拟处理的类器官使用兔抗NP抗体进行染色。染色前，样品用4%多聚甲醛在室温下固定15分钟，用0.5% Triton X-100透化10分钟，并用封闭缓冲液封闭1小时。为了评估病毒在类器官中引起的细胞病变效应（CPE），将感染和模拟处理的样品用兔抗NP抗体和小鼠抗TUBULIN抗体共染色，然后与荧光标记的二抗（山羊抗兔抗体和山羊抗小鼠抗体）孵育。用DAPI和Phalloidin-647分别复染细胞核和肌动蛋白丝。染色后，将整个类器官用Fluoromount-G™抗衰减封片剂封片，并用Leica STELLARIS系统成像。
流式细胞术
为了通过流式细胞术确定病毒感染率，在感染后10小时（MOI=1）使用10× TrypLE在37°C孵育3-5分钟将类器官解离成单细胞。然后将细胞用4%多聚甲醛（PFA）在室温下固定15分钟，用0.1% Triton X-100在4°C透化5分钟，并用兔抗NP抗体染色，随后与山羊抗兔二抗孵育。使用模拟感染的类器官设定门控策略。在BD Fortessa仪器上进行流式细胞术，并使用FlowJo（v10.9）分析数据。
分子动力学模拟
分子动力学（MD）模拟的初始结构包括甲型流感病毒（WSN/1933）PA亚基CEN结构域的野生型（WT）和D27G突变型，每个都与BXA和两个Mn2+离子复合。这些系统分别称为WT–BXA和D27G–BXA。两种结构均使用SWISS-MODEL服务器建模，以晶体结构（PDB ID: 6FS6）的B链为模板。
使用YASARA Dynamics在Ubuntu 20.04 LTS平台上对每个蛋白质-配体系统进行100 ns的MD模拟。设置涉及优化氢键网络和预测pKa值以微调pH 7.4下的质子化状态。然后将每个复合物在周期性边界条件下的立方体模拟池中溶剂化，蛋白质与池边缘之间的距离为10 Å。系统用0.9% NaCl中和，并额外添加Na+或Cl-离子。通过最速下降法和模拟退火进行能量最小化以消除空间冲突。在NPT系综（298 K, 1 atm）下进行100 ns的生产模拟，蛋白质使用AMBER14力场，配体使用GAFF2和AM1BCC，水使用TIP3P模型。范德华相互作用采用8 Å的截断值，而长程静电相互作用采用粒子网格Ewald方法处理。运动方程用多时间步长积分，键合相互作用为1.25 fs，非键合相互作用为2.50 fs。每200 ps保存一次轨迹帧用于后续分析。
为了评估系统的构象稳定性和动态行为，计算了均方根偏差（RMSD）、均方根涨落（RMSF）、回转半径（Rg）和蛋白质二级结构含量。此外，还使用边界元方法结合AMBER14力场估算了结合能。
统计分析
本研究中的数据处理、统计分析和绘图使用R（v4.4.1）、Python（v3.11）和GraphPad Prism（v8.0）进行。使用ggtree R包可视化系统发育树。所有定量数据均表示为平均值±标准差（SD）。通过非线性回归分析估算剂量反应曲线和EC50值。使用Shapiro–Wilk检验评估数据正态性。对于两组之间的比较，对正态分布数据使用非配对双尾Student t检验评估统计学显著性，对非正态分布数据使用Mann–Whitney U检验。对于多组比较，对正态分布数据进行多重t检验，然后采用Benjamini, Krieger and Yekutieli方法进行校正。表示P值 < 0.05，表示P值 < 0.01，表示P值 < 0.001，**表示P值 < 0.0001。
Results
2018年至2025年中国季节性流感PA基因表现出保守性
本研究基于国家传染病医学中心微生物研究网络和GISAID数据库的数据，分析了2018年4月至2025年1月期间在中国内地收集的共计3,938条H3N2和H1N1pdm09的聚合酶酸性（PA）基因序列。其中，484条序列来自该网络，3,454条序列检索自GISAID数据库。由于质量差，共有21个临床样本和67个GISAID样本被排除在进一步分析之外。排除后，保留2,392条H3N2 PA序列和1,458条H1N1pdm09 PA序列用于下游分析。大多数H3N2样本集中在2022年、2023年和2019年，而H1N1pdm09样本主要收集于2023年和2024年。临床样本主要集中在2024年和2025年，补充了最近几年的数据。由于COVID-19相关的社交距离措施，2021-2022年期间没有样本可用。患者的地理分布和人口统计学特征见补充表1。
先前研究表明，PA基因在流感病毒中高度保守。相应地，研究调查了中国内地近年循环流感病毒PA基因和蛋白质的进化模式，发现它们也表现出相对保守的进化模式。使用BEAST分析进化速率显示，H3N2 PA的平均替换速率为每位点每年2.94×10−3次替换（95% HPD: [2.66×10−3, 3.23×10−3]），而H1N1pdm09 PA的速率为2.92×10−3（95% HPD: [2.68×10−3, 3.16×10−3]）。这两个速率均低于HA基因报道的速率。香农
