2018-2021年菲律宾达沃地区狂犬病病毒的新基因组趋势

摘要

由狂犬病病毒(RABV)引起的狂犬病仍然是菲律宾的一个重大公共卫生问题，尽管人们努力控制它。要在2030年前消灭狂犬病，有效的监测策略至关重要。在这项研究中，我们使用来自达沃市和附近省份的基因组序列检测了RABV在达沃地区的进化和系统动力学。我们将RABV ARTIC方案用于牛津纳米孔高通量测序，以在有限的资源下优化工作流程效率。对比2019年6月至2021年6月采集的新病毒样本(n= 38)以2018年6月至2019年5月的基线样本(n= 49)，在系统发育树中观察到新的亚分支，表明与以前未检测到的旧变异体有差异。大多数新病毒属于亚洲SEA4_A1.1.1谱系，但也发现了菲律宾从未报道过的新(SEA4_B1和SEA4_B1.1)和新兴(SEA4_B1.1_E1)谱系。基线研究报告了来自相同地区的RABV分离株的系统地理学聚类。然而，这种模式在当前的生物监测中被打破了，在原始集群之外的区域检测到了变异。此外，我们的发现揭示了达沃市和邻近省份之间的重要传播途径，与基线研究中观察到的主要市内传播形成对比。这些结果强调了持续和及时的基因组监测的必要性，以监测遗传多样性的变化和新菌株的出现，以及追踪传播途径的变化。实施具有成本效益的下一代测序工作流程将有助于将基因组监测整合到狂犬病控制计划中，特别是在资源有限的情况下。不同部门之间的合作可以增强当地实验室和专家在基因组技术和分析方面的能力。

关键词:

犬狂犬病病毒; 全基因组测序; 纳米孔技术; 生物监测; 菲律宾

1.介绍

狂犬病是由狂犬病毒(RABV)引起的危及生命的急性脑脊髓炎狂犬病病毒主要通过家犬的咬伤传播给人类。狂犬病是一种古老的疾病，但仍然是最被忽视的疾病之一，特别是在发展中国家。它每年导致大约60，000人狂犬病死亡；50%的死亡是15岁以下的儿童，主要在亚洲和非洲。世界卫生组织正在支持一项到2030年消除犬类传播的人类狂犬病的运动。

菲律宾是东南亚主要的狂犬病流行国家，每年大约有200到300例由犬传播的狂犬病病例。为了预防和控制人类狂犬病，卫生部于2007年颁布了第9482号共和国法，即2007年《反狂犬病法》。《反狂犬病法》要求地方政府部门实施犬的大规模疫苗接种，控制干预措施，如扣押和绝育，现场控制和处置未登记、流浪犬和未接种疫苗的犬，并开展狂犬病预防和控制信息和教育活动。达沃市是该国最大的城市之一，拥有大约163，000只犬。2006年，达沃市强制实施了犬类大规模疫苗接种，而自2011年以来，已实施了控制干预措施，如对流浪犬进行扣押和绝育。尽管加强了这些控制措施的实施，近年来该市的狂犬病发病率并没有显著下降。

根除狂犬病病毒的一个主要挑战是缺乏一些国家的可靠监测数据，如该疾病流行的菲律宾。监测对于构建知情的传染病管理和控制规划至关重要[8].基因组数据正在成为一种有效的监测工具，通过了解病原体群体、确认疑似病例、区分毒株、识别感染因子以及阐明病毒的进化史和传播模式。由于高通量测序技术在过去二十年中的突破，产生微生物基因组不再被认为是困难或昂贵的。因此，全基因组测序(WGS)已被视为诊断学的明显和不可避免的未来。

2014年，牛津纳米孔技术公司(ONT)推出了一款名为MinION的便携式平价测序设备。MinION已被用于对人类和其他动物临床样本中的多种微生物进行测序，包括细菌、病毒和真菌[13].鉴于其小巧、灵活的物理存在和长期读取能力，MinION已被用于许多重要的微生物学研究，包括检测埃博拉病毒、寨卡病毒和耐万古霉素肠球菌。Brunker等人开发了狂犬病病毒基因组监测的测序工作流程提供的基因组数据揭示了肯尼亚、坦桑尼亚和菲律宾狂犬病传播的重要而独特的见解。

Bacus等人通过Sanger方法使用传统的全基因组测序，生成了2018年6月至2019年5月期间菲律宾达沃市及其邻近地区狂犬病病例的基线基因组数据，突出了犬狂犬病的本地动态。为了进一步监测这些地区狂犬病病毒的基因组进化和系统动力学，利用牛津纳米孔MinION测序平台，从2019年6月至2021年6月进行了基因组监测。参考基线数据和新确定的传播途径，新病毒的遗传多样性和进化位置提供了对该地区狂犬病干预措施的长期影响的更深入的见解，使当地政府能够评估和改进现有战略。

2.材料和方法

2.1.样本收集和处理

2.2.RABV基因组富集

2.3.文库准备

2.4.纳米孔测序

2.5.生物信息学数据处理和分析

3.结果

3.1.样本特征

在这项研究中，总共处理了81份在达沃CVO进行的荧光抗体试验(FAT)中检测出狂犬病毒阳性的报告样本。在81个样本中，30个样本是在2019年6月至12月获得的，35个是在2020年获得的，16个是在2021年1月至6月获得的(图1a)。达沃地区位于菲律宾棉兰老岛南部，由一个主要城市达沃市和五个省组成:达沃-德奥罗省、北达沃省、南达沃省、东达沃省和西达沃省。样本是从达沃市(76.54%)和邻近的南达沃省(4.94%)、北达沃省(14.81%)和达沃-德奥罗省(3.70%)报告的狂犬病病例中收集的(图1b–D)。

在81个FAT阳性脑样品中，66个(81.5%)通过RABV N RT-PCR得到证实，并通过多重PCR进行RABV全基因组富集。Brunker等人建立的RABV MinION测序方案的修改在本研究中采用(表1)。改变RNA提取试剂盒和天然条形码试剂盒以适应可用的资源。引入了其他修改来改进测序工作流程，该流程首先在试运行中进行了测试。对于多重PCR，最初一批纯化扩增子的核酸浓度太低而无法检测。随后采用触地PCR来增加引物特异性，其中反应温度从65℃降低0.1℃，直到达到62.5℃，从而产生浓度范围为1.39至108 ng/μL的纯化扩增子。来自2019年的最老样品(n= 9)从文库制备中排除，因为在PCR纯化后无法检测到核酸浓度。

表1.Brunker等人(2020)对RABV MinION测序工作流程的修改。

在文库制备的试运行期间，最初将约5 ng输入DNA用于末端制备和条形码连接反应(表1)。然后将条形码样本汇集在一起，总计30 ng遗传物质用于衔接子连接。在接头连接之前，合并的扩增子被纯化，导致浓度仅为0.7 ng/μL。在随后的运行中，输入的DNA从5 ng增加到50 ng。为了使纯化过程中的DNA损失最小化，使用短片段缓冲液(SFB)和80%乙醇洗涤珠子。还通过使用热循环器而不是加热块来应用精确的培养温度。通过这些修饰，四个批次中条形码扩增子库的核酸浓度被记录为0.84至3.78 ng/μL(表2)。总共57个样品(86.4%的样品输入)进行了纳米孔测序。

表2.合并的扩增子和文库的总核酸浓度。

a计数中不包括阴性质控品；b样品一式两份进行处理。

在测序过程中，碱基调用模式关闭时，运行时间延长至16小时或更长，而不是碱基调用模式打开时的6小时，以提高读数(表1)。为了最大限度地利用资源，流通池被回收用于两种用途。在使用相同的流动池连续排序的批次中，尽管进行了流动池清洗步骤，仍观察到来自前一轮的污染痕迹。例如，记录了第二批阴性对照的数千个读数，足以产生覆盖率> 70%的共有序列，尽管在文库制备期间合并前核酸浓度< 1 ng/μL。有问题的阴性质控品与上次运行的RABV阳性样品具有相同的天然条形码。为了研究该批次中可能存在的污染，使用bammix对重复使用的流动池中运行的样本的二元比对图(BAM)文件进行了进一步的核苷酸混合物检查。样品(n= 11)与前一批次重叠的条形码被标记为在多个位置具有可能的核苷酸混合物。因此，这些序列被丢弃，不用于下游分析。然而，重要的是要注意，在干净的样品中也观察到了核苷酸混合物，尽管只有几百个读数。考虑到运行时使用了独特的条形码和流动池，干净的序列仍包括在分析中。

3.2.RABV全基因组序列

产生了46种狂犬病病毒的共有序列(80.7%的测序样品)。由于非常低的序列阅读深度，一个样品通过RABV-GLUE未被确认为RABV，而其余45个样品被鉴定为属于亚洲SEA4分支(表A1)。然而，根据设定的阈值，45种病毒中只有38种具有可接受的全基因组覆盖范围和覆盖广度。一些样品(n= 7)追溯到2019年，覆盖率非常低，没有用于下游系统发育和系统动力学分析(表A1)。样本病毒的谱系名称，包括2018-2019年基线研究中的病毒，是基于MADDOG数据库中的全局参考集(表A2)。Rambaut等人的命名系统已应用。来自达沃地区的大多数序列被归入谱系亚洲SEA4_A1.1.1 (69.4%的旧序列和63.2%的新序列)，其出现于SEA4_A1.1亲代谱系(图2)。一些剩余的序列被怀疑属于从A1.1.1谱系下降的两个新谱系:SEA4_B1 (30.6%的旧序列和7.9%的新序列)和SEA4_B1.1 (26.3%的新序列)，而其中一个新序列(2.6%)被归入SEA4_A1.1.2谱系。请注意，A1.1.2和B1.1谱系仅在新样本中识别。有趣的是，MADDOG在新病毒中也发现了一个潜在的新出现但目前样本不足的谱系(亚洲SEA4_B1.1_E1)。

图2. 2018-2021年菲律宾达沃市狂犬病全基因组序列与MADDOG全球参考集的谱系分配。新谱系和现有谱系的等级关系显示在最左侧面板的旭日图中，从主要谱系(Asian_A1)向外扩展至次要谱系(Asian SEA4_B1.1)。中间的图显示了输入序列的系统发生树(n= 87)加上来自相关谱系的所有序列，根据旭日图的配色方案对尖端进行着色。来自基线研究的序列(n= 49)在最右边的红色圆圈中表示，来自该研究的新序列(n= 38)是绿色的，而来自MADDOG储存库中相关现有谱系的序列是灰色的(表A3)。3.3.RABV的遗传多样性和系统地理分布

进行了系统发育分析，以评估在达沃市及其邻近省份的不同地区内流行的狂犬病病毒株的遗传关系。2018-2019年基线研究的样本来监测病毒进化枝和传播的变化。通过贝叶斯推理生成了系统发育树，绘制了2018年6月至2021年6月共87个报告样本(图3)。观察到三个主要分支，每个分支由来自基线(以红色突出显示)和新序列数据(以绿色突出显示)的样本组成。在所有三个分支中，2019-2021年新取样的病毒相对于2018-2019年的基线序列形成了一个分歧的亚分支，在系统发育树中具有祖先位置。对于分支1、2和3，每个分支的最近共同祖先分别追溯到2014年、2013年和2010年。贝叶斯后验概率(≥0.95)显示在分支标签上。

图3.来自菲律宾达沃市和邻近省份的87份报告的RABV样本的贝叶斯系统发育。样品根据样品ID_地区/省_镇_取样日期贴上标签。来自菲律宾(KM148260)、巴西(KM594034、KM594033、KM594041和KM594028)和中国(JN234411)的参考序列的基因组数据被用作构建系统发育树的外群。主要分支在括号中表示，而新旧样本分别用红色和绿色文本表示。

为了直观显示样本在达沃地区各省的地理分布，推断了时间分辨的系统发育树。类似于贝叶斯推断的系统发育，三个主要分支是明显的。分支1由来自托里尔达沃市的40个样品组成(n= 6)、Tugbok(n= 8)、塔洛莫(n= 20)、卡兰(n= 1)、碧瑶(n= 1)、布汉金(n= 2)、马赛丽(n= 1)，而Poblacion(n= 1);4个来自北达沃的样本:塔古姆(n= 1)、帕纳波(n= 1)，而Samal(n= 2);3个样品来自达沃德奥罗:新巴丹(n= 1)和孔波斯特拉山谷(n= 2);和1个来自南达沃迪戈斯的样品(图4).这个分支从其父分支分化出来的推断日期是2014年9月30日，表明这些样本最近的共同祖先可以追溯到2014年。样品的时间分辨率与取样日期一致。

图4.菲律宾达沃市及邻近省份报告样本的树木时间系统发育(进化枝1)。样品根据样品ID_地区/省_镇_取样日期贴上标签，并根据城市或省进行颜色编码。

分支2中绘制的大多数样本来自达沃市的不同地区:Buhangin(n= 8)、Poblacion(n= 4)、Tugbok(n= 7)、塔洛莫(n= 8)，而卡兰(n = 1) (图5).来自邻近省份的新采样病毒，如苏洛普(n= 1)南达沃的萨马尔(n= 3)的北达沃，和新巴丹(n= 1)的达沃德奥罗也被发现在这个进化枝。进化枝2中所有样本的最近共同祖先可追溯到2014年，样本的时间分辨率与采样日期一致。

图5.菲律宾达沃市及邻近省份报告样本的树木时间系统发育(进化枝2)。样品根据样品ID_地区/省_镇_取样日期贴上标签，并根据城市或省进行颜色编码。

另一方面，分支3主要由来自达沃地区南部的南达沃的样品组成(n= 4)，尽管来自北达沃(Samal，n= 1)和达沃地区以外(哥打巴托市、n= 1)也被映射(图6)。达沃地区样本与哥打巴托毒株分离的推断日期是2010年5月25日，表明这些样本最常见的祖先可追溯到2010年。样品的时间分辨率也与取样日期一致。

图6.菲律宾达沃市及邻近省份报告样本的树木时间系统发育(进化枝3)。样品根据样品ID_地区/省_镇_取样日期贴上标签，并根据城市或省进行颜色编码。

3.4.达沃地区RABV的系统动力学

根据BEAST分析重建了系统发育的祖先状态，以阐明狂犬病在达沃地区的传播途径。如所示，共有16条路线被确定为具有强至决定性贝叶斯因子支持(BF > 10)表3和图7。在这16条路线中，有两条记录为达沃市内的市内传播，九条记录为达沃市和邻近省份之间的市间传播，一条记录为达沃地区省份内的省内传播，四条记录为跨省传播。重建的传播途径仅对本研究和Bacus等人捕获的RABV分离株具有特异性且不包括研究期间其他可能传播的变异体。

4.讨论

4.1.RABV多样性和传播模式的新趋势

我们通过牛津纳米孔技术使用NGS获得了2019年6月至2021年6月期间从菲律宾达沃市收集的38个RABV分离株的全基因组序列。将我们的新序列整合到2018年至2019年收集的达沃市49个RABV分离株的基线数据中显示，2019年至2021年期间收集的较新病毒从2018年至2019年的旧变种中分化出来，这表明出现了新的亚分支。这一假设得到了MADDOG谱系分析的支持，其中28.9%的2019-2021序列属于命名为SEA4_A1.1.2和SEA4_B1.1的新亚谱系；在达沃地区的基线研究中，以前没有在分离株中检测到这些。相反，基线研究中的分离株都由SEA4 B1谱系组成。2013年首次在吕宋发现A1.1.2谱系，而在2018-2021年检测到的达沃地区的B1和B1.1谱系在该国其他地方从未报道过。请注意，在新样本中检测到了A1.1.2和B1.1谱系，而不是基线样本，这表明它们在2019年后的某个时候出现在达沃地区。这一发现与我们的树时系统发育数据一致。然而，也有可能这些变异只是在基线研究中没有被捕获，因此这些谱系中的序列先前不可用。

我们还报告了一个新兴的但目前采样不足的谱系，SEA4_B1.1_E1。只有在五年内至少鉴定出10个样本时，才能确定新的谱系。只有五个样本与B1.1_E1血统相关，所有这些样本都是在2019年至2021年期间在菲律宾达沃地区分离出来的。因此，建议持续监控。据我们所知，我们是菲律宾第一个报道B1.1_E1血统的人，尽管目前还很难说这是仅在菲律宾达沃地区出现还是在菲律宾更广泛的地区出现。

用基线序列绘制我们的数据也揭示了最近的样本形成了新的亚支，其从更古老的祖先菌株中分化出来，而不是从后面的谱系中分支出来。这意味着可能存在至少一种未被发现的、未被确认的、更古老的变种在达沃市传播。新的和以前未检测到的变异可能在传播性、毒力和疫苗耐药性等方面具有先进的特征。因此，这些新分支的流行病学适合度以及它们是否会在达沃地区胜过其他分支将会引起人们极大的兴趣。

除了出现新的谱系和分支，还观察到传播途径的显著变化。在2018年至2019年的基线监测中，RABV变异体根据地理位置进行聚类。我们在2019年至2021年的后续几年中的基因组监测显示，以前与某个地理区域相关的变异已经被引入到其他位置。例如，分支1主要由来自达沃市和达沃地区北部邻近地区的样本组成。然而，我们的分析显示，这些变种现在已经到达了南达沃的迪戈斯市。另一方面，进化枝2最初仅在达沃市发现，但现在有来自其他省份的新病毒，如北达沃、德奥罗和南达沃。

这些观察结果得到了种系地理传播分析的验证，在种系地理传播分析中，跨区和跨省传播是明显的。与基线研究中主要的城市间传播相比，当前研究中的大多数移动发生在达沃市及其邻近省份之间，这表明长距离传输增强了。达沃市是一个主要的城市中心，位于棉兰老岛南部达沃地区的中心27].人类社会经济和人口密度之间的相互联系因素可能决定了RABV的传播动态。因此，达沃地区内频繁的农村到城市的病毒传播是非常可能的。然而，值得注意的变化也可能是由于达沃市以外地区样本提交范围的扩大(当前研究中约占样本收集总量的25%，而Bacus等人2021年的基线研究中为16%)。

在我们对北达沃岛的岛屿花园城市萨马尔的分析中也观察到了类似的情况。所有三个分支都包含来自Samal的样本，这表明达沃地区的各个地区都参与了该岛的传播事件。尽管水体可以成为狂犬病传播的强大屏障29]，萨马尔岛是棉兰老岛南部的热门旅游目的地，距离大陆只有几公里。因此，从达沃地区的几个地方到萨马尔岛的大量人口流动，使该岛非常容易受到人介导的岛间传播的影响。菲律宾的其他地区以前也报告过人为媒介的岛屿间传播而在泰国等其他国家和北非。

RABV的长距离传播很可能是由于人类介导的感染狂犬病的犬的转移，这是RABV在无狂犬病地区反复和成功再引入的重要驱动因素。多项研究报告称，疫苗接种活动后频繁的重新引入会抵消当地狂犬病的消除，从而使狂犬病在菲律宾等流行地区持续传播。在加纳最近的RABV变异株中也观察到了长距离传播，这些变异株属于非洲1a和非洲1b亚系，这与早期报道的单一谱系(非洲2)在西非传播相反。以前已知非洲1a型在非洲北部和东部占主导地位，而非洲1b型在非洲东部、中部和南部被发现，这表明RABV可能从非洲其他地区反复传入加纳。随着在达沃市和萨马尔观察到的长距离和岛屿间传播事件，地方政府单位(LGU)应探索干预措施，严格管制此类流动，以限制传播，并防止病毒的传入或再传入。例如，Samal在2011年即将被宣布为无狂犬病(个人交流，达沃市兽医办公室)，但我们目前的发现表明，薄弱的边境控制可能阻止了这一点。

我们的分析还显示，达沃市的塔洛莫区报告的病例数最高，使其成为该市的RABV热点。这支持了基线研究的发现，该地区已被确定为达沃市RABV传播链中的捐赠者和接受者。人口结构是狂犬病维持的驱动力，特别是在人和犬的人口密度高的城市地区。因此，狂犬病病例在大城市地区的出现和复发是由不断增加的人和犬的数量引起的。

4.2.低资源环境下的基因组监测

我们关于以前未检测到的和新出现的谱系和新传播途径的新发现强调了通过基因组监测持续监测狂犬病的重要性。与传统的临床报告不同，基因组监测可以提供对病毒特征和行为的更深入的见解，这对指导狂犬病的最后阶段至关重要。扩大监测范围，包括目标地区以外的样本，可以捕捉长距离传播途径，这将回答有关病毒循环和潜在宿主转移的关键问题，并提供对持续谱系的见解。

在许多狂犬病仍然流行的第三世界国家，下一代测序的资金机会有限。因此，必须最大限度地利用可用资源，高效的基因组监测工作流程非常重要。相对于基线研究中使用的费力费时的常规测序，MinION工作流为WGS提供了一个快速而经济的平台。在我们的研究中，86%的基因组丰富的样本进行了测序。在MinION工作流程优化后，这些样本中有67%成功生成了全基因组数据，其中通过增加输入DNA、利用替代技术(如触地PCR)和优化清理步骤提高了输出。

在Brunker等人开发的测序工作流程中使用处方试剂盒以外的试剂可能导致了不同的效率。这可以解释为什么在病毒RNA提取后，在我们的靶向RT-PCR证实步骤中，只有81%的报告RABV阳性样品检测为阳性。还必须考虑样品的年龄和储存条件。总共有15个样本在靶向RT-PCR过程中反复测试为阴性，其中大多数(73%)在进行研究之前被搁置了两年多。在测序的样本中，旧样本的全基因组覆盖率(14%)明显低于阈值90% 。这表明组织样本和病毒RNA本身的质量可能在储存过程中受到了损害。通过将样本保存在适当的保存缓冲液中，可以减少遗传物质在长期储存或运输过程中的损失。或者，在收集后立即对样品进行测序可能会增加百分产率。

重复使用流动池时应小心谨慎。我们发现，ONT储存缓冲液如果没有从废液通道中完全清除，会堵塞废液口，因为残留的缓冲液稍后会在流动池中结晶。此外，我们发现，当第一次运行后剩余的孔数小于400时，重复使用流动池会导致较低的读取数(即，16小时运行后Fast5文件较少)。我们还注意到重复使用的流动池中存在污染，即使在第一次使用后流动池与ONT清洗试剂一起孵育30分钟，在解复用步骤中仍会出现之前运行的天然条形码。这得到了我们使用bammix工具的分析的支持，该工具揭示了在所讨论的样品的多个位置中的核苷酸混合物的痕迹。令人惊讶的是，甚至在独特的流动池中运行的干净样品中也检测到核苷酸混合物(<总基因组的1%)。这可能表明R9.4.1流动池中纳米孔的读数效率和准确度可能下降，或者可能是由于碱基调用过程中的错误，因为本研究中选择了“快速”模式，而不是“高准确度”选项。然而，重要的是要考虑bammix工具使用BAM文件作为输入来计算可能的核苷酸混合物。medaka工具(https://github.com/nanoporetech/medaka(2022年12月8日访问))在本研究中用于从相同的BAM文件产生共有序列，这保证了不管任何碱基调用或纳米孔错误都调用真实变体，只要在每个位置有足够的覆盖。因此，清洁样品(即在独特的R9.4.1流通池中处理并使用独特条形码的样品)被处理用于随后的系统发育分析，前提是它们通过了本研究中设定的所有标准，并且在文库制备之前，与样品一起处理的阴性对照的浓度< 1 ng/μL。我们建议将含脱氧核糖核酸酶的洗涤缓冲液的孵育时间延长一小时或更长时间，以完全消除污染样品。如果可能，我们还建议在同一个流动池中连续运行时使用不同的条形码组，以避免条形码重叠。

引入高效且成本较低的NGS工作流程将有力推动将基因组监测纳入狂犬病控制规划，尤其是在资源匮乏的环境中。与此同时，还需要增强当地实验室的能力，使其能够利用基因组分析方面的专业知识实施该技术。这可以通过学术界、公共卫生和兽医机构的合作来实现。

5.结论

我们的研究报告了菲律宾南部狂犬病病毒基因组监测的随访结果。在这里，我们使用牛津纳米孔高通量测序技术从传统测序转向下一代测序，这使我们能够在更短的时间内高效、经济地批量处理样品。在我们的测序工作流程中引入的修改验证了该技术在低资源环境中的实用性，并且可以在其他实验室中采用。我们的发现揭示了达沃地区RABV的遗传多样性和传播模式的新趋势，强调了通过基因组监测进行持续监测的价值。高效且廉价的下一代测序工作流程可能有助于将基因组监测纳入狂犬病控制计划，尤其是在资源匮乏的情况下。这将需要当地实验室的能力建设和基因组专家的协助，这可以通过跨部门合作来实现。

Capin JBG, Sanque AJC, Eng MNJ, Lagare A, Sepulveda MCB, Murao LAE. Emerging Genomic Trends on Rabies Virus in Davao Region, Philippines, 2018-2021. Viruses. 2023 Jul 30;15(8):1658. doi: 10.3390/v15081658. PMID: 37632001.