全球大流行霍乱源头的演变

霍乱是由霍乱弧菌引起的急性腹泻感染，已造成七次全球大流行。前六次大流行由经典生物型引起，据信起源于恒河三角洲 ² 。上世纪初，以埃及埃尔托检疫站命名的埃尔托生物型开始取代经典生物型，成为流行性疾病的主要病原体。第七次霍乱大流行始于 1960 年的印度尼西亚，是所有病原体中持续时间最长的大流行 ³ ，近年来记录到的最大规模疫情包括 2010 年的海地(参考文献 ⁴ )和 2016 年(参考文献 ⁵ )及 2019 年(参考文献 ⁶ )的也门疫情。自 2022 年以来，霍乱疫情再次出现。黎巴嫩 ⁷ 和叙利亚 ⁸ 等国首次报告了霍乱病例。马拉维曾出现过非季节性疫情爆发 ⁹ ，孟加拉国在 2022 年经历了有史以来最大规模的疫情爆发（参考文献 ¹⁰ ）。

第七次大流行可归因于一个名为第七次大流行埃尔托谱系（7PET）的单一独立遗传谱系 ¹¹ 。尽管 7PET 在表型上与经典生物型（O1）属于同一血清群，但它是在获得 CTXφ噬菌体和两个致病岛（ 弧菌第七次大流行岛（VSP）-I 和 VSP-II）后，独立于非流行性埃尔托生物型祖先而出现的 ¹² 。自 1960 年以来，7PET 已在全球范围内传播，并经历了三次重叠的疫情浪潮，均起源于恒河流域 ¹ 。

7PET 的全球传播表现为反复的克隆扩增，已观察到 11 次传入非洲的事件（分别称为 T1 和 T3-T12） ¹³ ，3 次传入拉丁美洲（LAT1-LAT3） ¹⁴ ，8 次传入欧洲（EUR1-EUR8） ¹⁵ 。尽管霍乱弧菌在特定疫情或特定地理区域的动态变化已被研究 ^6,7,13 ，但尚未对恒河流域（横跨印度和孟加拉国，被誉为“全球霍乱疫源地”）的整体进化动态进行纵向研究。已知的是，过去 15 年间，在印度和孟加拉国持续观察到两个 7PET 谱系，此前分别命名为“BD1”和“BD2”（参考文献 ^{10,16,17,18,19,20} ）。尽管它们的基因组差异小于 150 个单核苷酸多态性（SNP），但它们在可移动遗传元件（MGE）方面表现出显著差异，从而导致噬菌体防御、抗生素耐药性和毒力方面的差异。这些包括基因组岛（例如 CTXφ、VSP-II 和超级整合子 ²¹ ）、SXT 整合和接合元件 (SXT-ICE)（例如 ICE ^TET 和 ICE ^GEN (参考文献 ²² )）、质粒（例如 pCNRVC190243；参考文献 ⁶ ）以及噬菌体诱导的染色体岛样元件 (PLE) ²³ 。尽管许多此类遗传元件在疾病严重程度中所起的作用尚不清楚，但最近的宏基因组分析表明， 霍乱弧菌特异性噬菌体 ICP1 的存在可以减轻疾病的严重程度 ²⁴ 。

我们收集并测序了来自孟加拉国（ n = 1,516；2014–2023 年）和印度北部（ n = 794；2002–2023 年）的霍乱分离株，据我们所知，这是迄今为止恒河流域最全面的霍乱纵向数据集。多年来，人们一直认为霍乱的主要全球来源是恒河三角洲和孟加拉湾的咸水。然而，通过追踪恒河流域的霍乱弧菌 ，我们观察到霍乱的传播受到国界限制，因此可能受到人口流动的影响，而不是像临床疾病主要通过环境传播途径所预期的那样。我们推断，尽管霍乱在孟加拉国的流行率很高，但其在全球的传播很可能是由印度的输出造成的。在孟加拉国境内， 霍乱弧菌展现出独特的进化模式，特别是抗噬菌体元件的快速丢失/获得，这与米汤样便或严重脱水的风险增加有关。此外，这些遗传元件的存在可能会阻碍霍乱弧菌从孟加拉国传播到其他国家。

我们对 2014 年至 2018 年孟加拉国全国系统性霍乱监测研究中收集的菌株（ n = 1453）以及在达卡孟加拉国国际腹泻病研究中心（icddr,b）医院开展的全肠道疾病筛查活动中收集的菌株（ n = 63）进行了测序。在 icddr,b 医院的筛查活动中，每 50 名就诊患者中就有一人被纳入研究并接受肠道病原体检测，最终获得了 1477 个高质量基因组。此外，我们还对来自印度的 794 个基因组进行了测序，这些基因组分布于印度北部八个邦，并通过昌迪加尔医学教育与研究研究生院（PGIMER）提供的转诊、临床和监测服务收集。将我们的分离株置于已发表的 3112 个全球基因组（补充数据 1 和补充图 1 ）的背景下，我们发现来自孟加拉国和印度的分离株主要属于两个系统发育分支，分别对应于先前描述的谱系 BD1 和 BD2（参考文献 ¹⁶ ）。随着来自恒河流域和三角洲的更多数据，我们清楚地认识到，先前命名的 BD1 分支不仅包含先前定义的 BD1 谱系，还包含传播谱系 T11-13 和 LAT3。我们将该谱系称为超 BD1（sBD1）。在 sBD1 和 BD2 内部，样本对之间的 SNP 中位数分别为 37 和 19 个；而在 sBD1 和 BD2 之间，样本之间的 SNP 中位数为 112 个。为了提供追踪霍乱弧菌在恒河流域和三角洲传播所需的分辨率，我们根据 SNP 距离阈值 20 将 sBD1 和 BD2 细分为离散的亚谱系，并根据它们在全球系统发育树中的位置进行编号（补充图 2 ）。

根据两项孟加拉国监测研究，BD2 在 2014 年至 2017 年间在孟加拉国占主导地位（图 1a ），之后在 2018 年被 sBD1 取代。根据祖先重建推断，sBD1 是从印度重新引入的（扩展数据图 1 ）。值得注意的是，这是本研究中观察到的孟加拉国唯一一次主要的病毒引入事件；在孟加拉国流行的所有其他主要亚系均与前一年的亚系相同或源自前一年的亚系。各亚系的持续时间各不相同（图 1a ），在 2014 年至 2018 年的系统监测研究中，平均持续时间为 16 个月（四分位距：8-24 个月）。在监测研究期间于孟加拉国观察到的 16 个亚系中，有 15 个在多个行政区划中被检测到。考虑到病原体迁移速率，从首次检测到亚谱系扩散至孟加拉国所有八个行政区之间存在 12 个月的滞后（线性模型；95%置信区间：9-16 个月）（扩展数据图 2 ）。为了推断迁移方向，我们进行了重复子抽样，以消除各地区急性水样腹泻病例抽样差异的影响。对于每个子样本，我们构建了 sBD1 和 BD2 的时间尺度系统发育树，并推断了每个节点的祖先状态。推断的区域间传播事件主要来自达卡（占 49%），其次是吉大港（占 24%；补充图 3 ）。

The alternative text for this image may have been generated using AI.

Full size image

a ，2014-2018 年全国监测研究和 2022-2023 年 2%研究中各亚系丰度随时间的变化，按亚系着色（见图例）。2014-2018 年监测研究和 2%研究均未开展的时间段以白色阴影显示。b ， 2014-2018 年全国监测研究和 2022-2023 年 2%研究中各基因/移动遗传元件（MGE）谱的丰度随时间的变化，按谱型着色（见图例）。2014 – 2018 年监测研究和 2%研究均未开展的时间段以白色阴影显示。c、 d ，从达卡（ c ）和吉大港（ d ）传播的不同霍乱弧菌亚群的系统地理学。为避免抽样偏差，仅纳入 2014-2018 年全国监测研究中的样本。对于每个行政区划，箭头指示了各亚群体首次引入该区的来源。箭头颜色代表所关注的亚群体（见图例）；箭头大小表示传播事件的相对数量，透明度表示相对于首次传播事件的月份（见图例）。图 c 和图 d 中的底图（包含行政区划边界）来自 R 语言的 bangladesh 软件包。

除了核心基因组中的单核苷酸多态性（SNP）变异外， 霍乱弧菌基因组还发生了单个基因和整个移动遗传元件（MGE）的丢失和获得（图 1b 和扩展数据图 3 ）。在 2014 年监测研究的第一年，由于 17 个碱基对（bp）的缺失和由此导致的移码突变，81%的 BD2 基因组中完整的 hlyA 外毒素基因（VCA0219）已经丢失（扩展数据图 4a ）。此外，VSP-II 内的操纵子 ddmABC 也发生了连续降解，该操纵子在噬菌体感染或质粒摄取过程中会触发细胞凋亡 ²⁵ 。这包括 ddmC 中的移码突变，这可能导致整个 ddmABC 操纵子失去功能，随后是 ddmB 中的转座酶插入，以及 ddmAB 和邻近基因 VC0493 – VC0494 的缺失（扩展数据图 4b ）。随后， wbeT 基因（ VC0258 ；扩展数据图 4c ）发生破坏，导致 2016 年表型血清型转换；2017 年，两个不同的亚系分别独立丢失了抗噬菌体岛 PLE1；同年，SXT-ICE ^TET 元件也丢失（图 1b ）。此外，在 2016 年，一个单系群（其超整合子中 VCA0455 – VCA0459 缺失，扩展数据图 4d ）逐渐增多，对应于亚系 BD2.028 和 BD2.029，占 2016 年测序的 BD2 分离株的 51%，但在 2017-2018 年下降至 27%。

2018 年，sBD1 亚系 sBD1.070 在孟加拉国完全取代了 BD2。与 BD2 不同，sBD1.070 携带 SXT-ICE ^GEN 元件、 ctxB7 等位基因、完整的 hlyA 外毒素基因、完整的 ddmABC 操纵子和完整的 wbeT 基因 。然而，由于无义突变，sBD1.070 的定植因子 acfC （ VC0841 ；扩展数据图 4e ）存在移码突变，并且最初缺乏 PLE 抗噬菌体元件。到 2022 年，达卡开始进行 2%监测研究时，已有三个新的亚系取代了 sBD1.070，其中两个（sBD1.067 和 sBD1.069）源自 sBD1.070，另一个（sBD1.060）可能来自印度。 sBD1.060 亚系与近期在巴基斯坦、黎巴嫩和马拉维爆发的疫情的亚系亲缘关系更为密切。在 sBD1.067 中， wbeT 基因被插入序列破坏，使其具有稻叶血清型。此外，sBD1.067 和 sBD1.069 均获得了一个新的 PLE 岛，本文将其命名为 PLE11，其基因组位置与 BD2 中的 PLE1 不同（扩展数据图 4f ）。

接下来，为了探究孟加拉国基因功能丧失和获得的时间和空间特征，我们仅分析了 2014-2018 年系统性监测中的基因组，以避免抽样偏差。图 1b 和扩展数据图 3 显示，基因丢失存在明显的时空模式：亚谱系 BD2.032 中的 PLE1 和 ICE ^TET ，BD2.029 中的 PLE1，BD2.026 中的 ddmAB ，以及 BD2.028 中的 VCA0455 – VCA0459 。这些变化在所有衍生亚谱系中均已固定（补充图 2 ）。利用 TreeTime，可以明显看出两种地理传播模式：sBD1.3.1、PLE1 ⁻ 和 ICE ^TET− BD2.2.2 和 Δ ddmAB BD2.2.1 从达卡传播（图 1c ），而 PLE1 ⁻ BD2.4.3 和 Δ VCA0455 – VCA0459 BD2.4.1 从吉大港传播（图 1d ）。

此外，还观察到零星获得的移动遗传元件（MGE），包括质粒 pSA7G1（扩展数据图 5a ），该质粒在分布于六个行政区的五个 BD2 亚系中多次获得。pSA7G1 的功能相关性尚不明确。溶原性噬菌体 K139 可与霍乱弧菌 O1 抗原结合，并携带与小鼠模型毒力相关的 glo 基因 ^26,27 ，在七个行政区的十个 BD2 亚系中均有发现，其中在吉大港最为常见（扩展数据图 5b ）。值得注意的是，从 15 个霍乱弧菌样本中检测并测序了裂解性噬菌体 ICP1，与 K139 类似，该噬菌体在吉大港采集的样本中最为常见（扩展数据图 5c ）。

为了解孟加拉国与恒河流域其他国家霍乱弧菌之间的关系，我们分析了在印度北部收集的 794 株分离株基因组（扩展数据图 6 ）。尽管孟加拉国和印度都存在属于 sBD1 和 BD2 亚系的菌株，但两国亚系更替的时间模式截然不同（图 2a ）。2004 年至 2011 年间，sBD1 和 BD2 在孟加拉国和印度共同流行，之后 sBD1 在印度占据主导地位，到 2011 年占样本总数的 94%；而 BD2 在孟加拉国占据主导地位，到 2013 年占样本总数的 95%。2013 年至 2018 年间，孟加拉国共存在 11 个 BD2 亚系，其中 10 个亚系仅存在于孟加拉国境内，未在世界其他任何国家发现，这表明在此期间，BD2 在孟加拉国独立进化。在两国，Tajima’s D 值均持续低于 -4，提示存在选择性清除，且核苷酸多样性随时间变化，在多个 7PET 谱系同时存在的时期有所增加（补充图 4 ）。在孟加拉国流行的亚谱系分离株所特有的 PLE 在该国以外地区较为罕见， ddmABC 基因的丢失和 Inaba 血清型的完全替换也同样罕见（补充图 5 ）。

The alternative text for this image may have been generated using AI.

Full size image

a ，孟加拉国、北印度和加尔各答各亚谱系随时间变化的丰度，包括本研究采集的样本和先前发表的背景样本。仅显示样本数在十个或以上的亚谱系。b，南亚地区推断的传播事件和关键基因/移动遗传元件（MGE）获得事件。箭头颜色代表亚谱系，大小代表我们全球样本库中源自该传播事件的样本数量。c ， 2003 年至 2023 年间推断的从南亚到其他地区的传播事件。箭头颜色代表亚谱系，大小代表我们全球样本库中源自该传播事件的样本数量，不仅包括所示国家，还包括向其他国家的进一步传播。b 和 c 中的底图来自 Natural Earth（ https://www.naturalearthdata.com/ ），遵循通用公共领域 CC0 许可协议。

加尔各答距离孟加拉国库尔纳地区仅 60 公里，且与孙德尔本斯湿地毗邻，罗伯特·科赫于 1884 年提出，孙德尔本斯湿地可能是霍乱的起源地 ²⁸ 。然而，在加尔各答发现的霍乱弧菌的时间模式与印度北部的情况一致（图 2a ）。相反，2015 年至 2018 年间在库尔纳地区采集的 167 个基因组样本中，没有一个属于与加尔各答共享的亚系，而全部 167 个都属于与孟加拉国其他地区共享的亚系（扩展数据图 7a ），这表明霍乱的传播遵循的是边界而非水文特征。尽管阿萨姆邦地理位置靠近孟加拉国，但 2002 年至 2005 年间在阿萨姆邦采集的 6 个样本形成了一个系统发育簇，该簇与其他印度样本的亲缘关系最为密切（扩展数据图 7b ）。 2010 年在尼泊尔采集的 12 个样本中，12 个属于当时印度存在的亚谱系，没有一个属于孟加拉国存在的亚谱系（扩展数据图 7c ）。

接下来，为了根据全球 7PET 的时间尺度系统发育树（扩展数据图 1 ）重建全球传播事件的历史，我们推断了每个节点的祖先状态和日期（补充数据 2 ）。南亚国家之间的传播率很高，但孟加拉国和印度之间的传播事件（408 起中的 40 起）比其他任何国家对都多（图 2b ）。更详细地分析这些事件，我们发现 sBD1 和 BD2 在 20 世纪 90 年代末期分别从 T9 谱系中分化出来。此外，sBD1 似乎在印度获得了 ctxB7 等位基因，随后传播到尼日利亚、肯尼亚和海地，并导致了 2010 年海地疫情的爆发（图 2c ）。自那以后，大多数来自南亚的导致疫情爆发的全球传播事件似乎都源于在印度进化的 ctxB7⁺ sBD1 基因库（图 2c ）。从 2003 年到 2023 年，我们全球系统发育树中来自孟加拉国的样本数量是来自印度的两倍，然而，南亚以外地区近期源自印度传播事件的样本数量（ n = 359；2003–2023）是源自孟加拉国传播事件的样本数量（ n = 11；2003–2023）的 32.6 倍。为了检验这是否是由于不同时期采样强度的差异造成的，我们进行了重复子采样，使每年来自孟加拉国、印度和世界其他地区的样本数量保持一致（补充图 6 ）。结果再次表明，自 2000 年以来，源自印度的传播事件数量更多。

尽管印度和孟加拉国之间发生了大量本地传播事件，但 40 起事件中有 29 起未能产生超过 10 个后代样本。2010 年代，这些数据表明发生了 8 起印度到孟加拉国的传播事件，其中 7 起未能持续传播，导致孟加拉国的患者样本不足 10 个。虽然 BD2 在 2013 年至 2017 年间成为孟加拉国的主要毒株，但自 2010 年以来仅观察到一起输出事件，导致在加尔各答检测到一个 BD2 样本。BD2 和 sBD1 均在传入或重新传入孟加拉国后获得了 PLE。BD2 大约在 2001 年传入孟加拉国，并在两年后在孟加拉国获得了 PLE1。同样，sBD1 在 2020 年左右获得了 PLE11，此时距离 sBD1 在孟加拉国的样本中占比超过 50%已有两年。除孟加拉国外，未在其他任何国家检测到 PLE11。

鉴于基因/MGE 丢失的快速动态，但直接源自孟加拉国的全球传播事件数量相对较少，我们研究了基因和 MGE 的变化如何影响适应性（感染潜力、对噬菌体的易感性和抗菌敏感性），以及不同谱系向孟加拉国和恒河三角洲以外地区扩散的潜力。

我们首先进行了多变量逻辑回归分析，以确定基因或移动遗传元件（MGE）的缺失是否与临床严重程度相关 （ 图 3a ），并校正了其他基因和 MGE 的存在。ddmA（ P = 0.01）、 wbeT （ P = 5.8 × 10⁻⁵）和 ICE（ P = 0.005）的存在与米汤样便风险降低相关，而 VCA0455 – VCA0459 （ P = 0.0001）和 PLE1（ P = 0.04）与严重脱水风险降低相关，提示这些基因的连续缺失可能加重疾病严重程度。我们认为，这些遗传元件的缺失可能促进其更快的复制，并在肠道内赋予其选择优势。我们对来自 21 个霍乱流行国家的 13105 个肠道宏基因组进行了系统性荟萃分析（补充表 1 和补充数据 3 ），发现当 ICP1 缺失时，携带 ctxB1 ⁺ （因此可能属于 BD2 谱系；图 3b）的霍乱弧菌丢失 ICE TET 后，粪便宏基因组中霍乱弧菌的 reads 比例从 16.8%中位数增加到 40.5%。即使在 ICP1 存在的情况下，霍乱弧菌丢失 ICE TET 后，粪便宏基因组中霍乱弧菌的比例也从 1.5%增加到 25.8%。

The alternative text for this image may have been generated using AI.

Full size image

a ，不同基因和移动遗传元件（MGE）与米汤样便和严重脱水关联的对数转换优势比，校正了其他基因和 MGE 的存在与否以及采集日期和地点。中心点表示计算的对数转换优势比，误差线代表 95%置信区间； n = 1617 个生物学独立样本。b ， Kraken 在粪便宏基因组数据中分配给霍乱弧菌的 reads 百分比。箱线图表示中位数和四分位距，须线延伸至 1.5 倍四分位距内的最极端值。样本根据 ctxB7 和 ctxB1 的存在情况进行分类，在该时间段和地点，ctxB7 和 ctxB1 分别代表 sBD1 和 BD2； ctxB1 ⁺ 样本根据 PLE1 和 ICE ^TET 的检测结果进行亚分类； n = 230 个生物学独立样本。 c ，孟加拉国样本中检测到的 ICP1 抗防御系统的频率，按年份和共测序霍乱弧菌的谱系划分 。d ，基因组中存在溶原性噬菌体 K139 的样本比例，根据 ddmABC 抗病毒系统中存在的基因进行分类。P 值表示多元逻辑回归检验 ddmABC 基因座变化与 K139 之间关联性（Wald 检验；双侧）的结果，并根据图 3a 所示的协变量（ICP1、pSA7G1、SXT-ICE、PLE、 wbeT 、 VCA0455 – VCA0459 、 hlyA 、谱系、日期和地点）进行调整 。 未进行多重比较校正。e，PLE ⁺ 和 wbeT ⁺ 霍乱弧菌的代表。 2003 年至 2023 年，孟加拉国样本总数中霍乱的推断出口事件以及由此出口事件导致的孟加拉国境外样本数量。

我们的数据表明，ICP1 类型与循环霍乱弧菌谱系之间存在显著的时间关联。在 BD2 被 sBD1 取代的过程中，我们观察到与其共测序的 ICP1 噬菌体类型也发生了同步转变：从携带 CRISPR-Cas 抗防御系统的 ICP1（2018 年 4 月之前的 13 个样本中有 11 个）转变为携带 odn 抗防御基因的 ICP1（2018 年 4 月之后的 6 个样本中有 5 个）（图 3c ）。为了更好地了解 ICP1 的全球分布和特异性，我们对上述分析的 13,105 个肠道宏基因组进行了 ICP1 筛查。我们在孟加拉国的 87 个宏基因组和印度的 13 个宏基因组中鉴定出了 ICP1，其中大部分存在于同时含有霍乱弧菌的宏基因组中（88%；扩展数据图 8a ）。我们从这些宏基因组中生成了 85 个覆盖率超过 50%的 ICP1 基因组，并将其与（1）公开可用的 ICP1 组装序列和（2）与我们全球霍乱弧菌基因组库共测序的 ICP1 整合，构建了全球 ICP1 系统发育树（扩展数据图 8b 和补充数据 4 ）。我们发现来自孟加拉国、印度、刚果民主共和国和也门的 ICP1 形成了不同的系统发育簇。携带 CRISPR-Cas 抗防御系统的 ICP1 在系统发育上与仅在孟加拉国发现且主要与 BD2 共存的 odn ⁺ ICP1 明显不同。Odn ⁺ ICP1 在所有四个国家均有发现，并且似乎已适应于 sBD1 谱系。 2018 年在孟加拉国出现的 odn ⁺ ICP1 与 2017 年在加尔各答发现的 odn ^{+ ICP1 的亲缘关系比之前在孟加拉国分离的菌株更为密切。我们还发现了与 K139 噬菌体易感性相关的基因改变。 值得注意的是， ddmABC 降解与 K139 流行率的下降显著相关（ ddmB 和 ddmA 的 P 值分别为 3.3 × 10 −6 和 0.01；逻辑回归调整了日期、研究地点以及其他基因的存在与否；图 3d ）。}

鉴于在孟加拉国以外地区未观察到 PLE11，我们随后研究了 PLE 是否普遍影响全球传播潜力。与在孟加拉国采集的 PLE ⁺ 基因组数量相比，PLE 在 2003 年至 2023 年间传播到其他国家的事件中显著偏低（图 3e ； P = 1.8 × 10⁻¹；Fisher 精确检验），并且在这些事件产生的霍乱弧菌基因组中偏低更为显著（ P = 2.2 × 10⁻²；Fisher 精确检验）。2003 年至 2023 年间从印度输出的 26 株霍乱弧菌中，没有一株（0/26）为 PLE ⁺ 。同样，缺乏完整 wbeT 基因的稻叶型霍乱弧菌也很少从孟加拉国传播出去（ P = 8.4 × 10⁻⁴；Fisher 精确检验）。

最后，正如预期的那样，表型抗生素敏感性测试表明，ICE ^TET 的丢失与对甲氧苄啶/磺胺甲噁唑（ P = 9.3 × 10 ⁻⁸ ；Fisher 检验；扩展数据图 9 ）、四环素（ P = 9 × 10 ⁻⁴ ）和多西环素（ P = 2.5 × 10 ⁻⁵ ）的敏感性增加相关，而携带 ICE ^GEN 元件的 sBD1 谱系霍乱弧菌对多西环素的敏感性显著更高（ P = 3.7 × 10 ⁻²² ）。

孟加拉国的孙德尔本斯和恒河三角洲人口稠密地区毗邻孟加拉湾，被许多人认为是霍乱疫情的中心地带 ² 。基因组研究进一步证实了这一观点，研究表明，一种名为 7PET 的霍乱弧菌克隆株从该地区扩散，引发了数波疫情，导致了持续至今的第七次霍乱大流行 ¹ 。7PET 是目前全球霍乱病例激增的罪魁祸首，也是世界卫生组织将霍乱疫情复燃列为三级紧急事件的原因 ²⁹ 。尽管霍乱的全球蔓延显而易见，2023 年全球霍乱病例达 66.7 万例，世界卫生组织六个区域中的五个都报告了疫情暴发（参考文献 ³⁰ ），但人们对霍乱弧菌在全球起源地的进化动态知之甚少。我们利用基因组学方法，研究了霍乱弧菌在大恒河流域（包括孟加拉国和印度境内的三角洲地区以及覆盖印度北部八个邦的恒河上游流域）的动态变化。我们追踪了该时期该地区主要的 7PET 循环谱系 sBD1 和 BD2 的演化过程。

我们的研究结果表明，尽管恒河三角洲和流域横跨孟加拉国和印度，但霍乱的演化似乎更倾向于沿着国界线而非水文特征和水流方向发展。加尔各答和阿萨姆邦都存在北印度特有的霍乱亚系，而孟加拉国则没有。此外，孟加拉国的 BD2 型霍乱亚系演化几乎与世界其他地区完全隔绝。而且，BD2 型霍乱的许多基因变异似乎是从人​​口密度最高的达卡向全国扩散的。虽然东印度样本可能存在偏倚，但这进一步证实了 7PET 型霍乱主要通过短周期人际传播。

在孟加拉国， 霍乱弧菌的进化以基因和移动遗传元件（MGE）的快速变化为特征，尤其是与噬菌体防御相关的基因和 MGE。携带这些变化的亚系迅速成为优势菌株，表明霍乱弧菌的进化受到强烈的选择压力驱动。2018 年霍乱弧菌噬菌体防御系统和 ICP1 抗防御系统几乎同时发生的转变就印证了这一点。ICP1 可能在 BD2 中与噬菌体防御系统共同进化后，保留这些系统不再具有进化优势，导致它们依次丢失。PLE1 在 BD2 进化过程中多次独立丢失，表明维持 PLE1 基因存在沉重的负担。抗噬菌体 ddmABC 操纵子也观察到了类似的模式；一旦 ddmC 中的移码突变导致该操纵子失效，其余部分就会被删除。我们发现 ddmABC 操纵子的降解与溶原性噬菌体 K139 的流行率降低有关。虽然从本文提供的数据来看，致病机制尚不明确，但对该复合物的结构和功能的持续研究可能会提供潜在的解释 ^32,33,34 。

BD2 谱系失去防御能力后，迅速被 sBD1 取代。sBD1 携带一个 ICE ^GEN 元件，其中包含一个抗噬菌体 BREX 系统 ¹⁰ ，该系统最初可能对本地 ICP1 有效。然而，携带 odn 抗防御系统的 ICP1（可能与 sBD1 一同从印度引入）随后迅速占据主导地位。最终，sBD1 通过获得抗噬菌体岛 PLE11 克服了这一挑战。PLE11 通过破坏 ICP1 尾部组装来阻止其增殖，从而成为 ICP1 进化的强大驱动力 ³⁵ 。此外，我们的数据表明，携带这些噬菌体防御系统与较低的疾病严重程度和霍乱弧菌载量相关，这表明维持这些防御系统会显著降低致病性。然而，由于我们对疾病严重程度的分析仅限于孟加拉国的样本，因此尚不清楚其他地区是否也存在同样的模式。赋予四环素耐药性的 ICE ^TET 的丢失表明，抗生素治疗并未施加如此强大的进化压力。

尽管在孟加拉国以外地区，PLE 抗噬菌体岛较为罕见，但 BD2 和 sBD1 在孟加拉国流行约两年后均获得了 PLE。此外，PLE 在孟加拉国以外的传播事件中出现频率较低，这表明它们可能仅在恒河三角洲地区具有优势，并会阻碍南亚地区内外的远距离传播。这或许可以部分解释为何推断的大多数全球传播事件的源头似乎是印度而非孟加拉国，从而修正了之前关于恒河三角洲拥有所有流行性 7PET 霍乱弧菌亚系全球多样性的观点。相反，包括印度在内的周边国家才是全球传播的跳板。

目前尚不清楚为何噬菌体噬菌体（PLE）似乎主要局限于孟加拉国，因为在印度、也门和刚果民主共和国也发现了 ICP1 ^6,36,37 。我们的假设是，孟加拉国是少数几个人口密度和霍乱流行率都足够高的地区之一，同时由于缺乏预防性的水、环境卫生和个人卫生（WASH）基础设施，导致暴露率很高，这使得噬菌体能够在人群中有效传播，同时维持昂贵的抗噬菌体防御系统。此外，我们发现稻叶型霍乱弧菌很少在孟加拉国以外传播，这表明，虽然血清型转换可能使小川特异性血清阳性率高的人群获得免疫逃逸 ^38,39 ，但它也可能给免疫力低下的人群带来适应性代价。

总之，我们发现孟加拉国霍乱弧菌的进化独特之处在于其移动遗传元件（MGE）的快速获得和丢失。这一过程遵循相对可预测的模式：霍乱弧菌在孟加拉国定殖时获得噬菌体相关元件（PLE），随后随着时间的推移而丢失，这可能是由于维持针对 ICP1 的高成本抗噬菌体防御机制与致病性之间的权衡所致——经典的基因对基因理论，但在此是在区域尺度上体现 ^22,40,41 。此外，在孟加拉国， 霍乱弧菌似乎从人口密度高的地区（例如达卡和吉大港）辐射扩散。对这些地区的霍乱弧菌进行实时纵向基因组监测，可以作为全国的早期预警系统，识别携带预计会增加重症疾病或抗生素耐药性风险的变异的新兴亚系。这有助于迅速调动干预策略，例如疫苗接种或改善水、环境卫生和个人卫生（WASH），以防止高风险亚系的传播。鉴于抗生素耐药性威胁日益加剧，噬菌体疗法作为一种新的干预策略，目前备受关注。尽管我们的研究结果应该与将 ICP1 用作治疗剂联系起来考虑，但其他噬菌体的共同进化能力可能有所不同。

我们还发现，移动遗传元件（MGEs），例如 PLEs，可能会损害霍乱弧菌从孟加拉国向恒河上游地区的全球传播潜力和输出。因此，尽管孟加拉国霍乱流行率高且种类繁多，但恒河三角洲可能并非当今全球霍乱的主要来源地，印度和整个恒河流域才是。我们的数据表明，仅将霍乱监测重点放在非洲、中东和海地的疫情上不足以终结第七次霍乱大流行。必须建立类似于目前流感监测和应对系统的全球年度霍乱监测机制，才能利用 WASH 干预措施和有限的口服霍乱疫苗储备，有效锁定传播源。

伦理声明

2014-2018 年孟加拉国监测研究和 2% icddr,b 达卡医院监测研究的研究方案已获得 icddr,b 研究审查委员会和伦理审查委员会的批准。所有参与者均签署了知情同意书，同意收集数据和样本。在印度北部开展的研究已获得昌迪加尔医学教育与研究研究生院（PGIMER）伦理委员会和合作委员会的批准。所有参与者在收集数据和样本前均签署了知情同意书。

孟加拉国的样本采集和测序

2014-2018 年监测研究的具体流程详见参考文献 ⁴² ，并已获得国际腹泻病研究中心（icddr,b）研究审查委员会和伦理审查委员会的批准。2014 年至 2016 年，在 10 个监测点开展了哨点监测。由于资金缺口，2016 年 1 月至 5 月期间监测工作暂停，随后自 2016 年 5 月起，监测范围扩大至 21 个地区的 22 个监测点。所有成年参与者均签署了知情同意书，18 岁以下儿童的知情同意书由其法定监护人签署。采集粪便样本用于检测霍乱弧菌 O1/O139。在 2014 年建立监测的 10 个监测点，还对样本进行了产肠毒素大肠杆菌 、 沙门氏菌和志贺氏菌的检测。在 26,221 例急性水样腹泻患者中，6.2%（ n = 1,604）确诊为霍乱弧菌 O1 型感染，其中 1,526 例进行了全基因组测序（补充图 7 ）。同样，在孟加拉国达卡国际腹泻病研究中心（icddr,b）开展的 2%达卡医院监测项目中，每 50 名到访该中心医院的患者均被纳入研究并接受肠道病原体检测，2022 年至 2023 年间收集的 63 份霍乱弧菌样本也进行了测序。

从哨点采集的粪便样本在采集后 15 天内用 Cary-Blair 培养基运送至 icddr,b ⁴² 。样本直接划线接种于牛磺胆酸-碲酸盐明胶琼脂平板上，并在碱性蛋白胨水（pH 8.6）中富集培养 18 小时，之后再接种于牛磺胆酸-碲酸盐明胶琼脂平板上。这些平板在 37℃下培养过夜。疑似霍乱弧菌菌落使用 O1 Ogawa 特异性抗体、O1 Inaba 特异性抗体和 O139 特异性抗体进行血清分型。每五个霍乱弧菌阳性培养物中抽取一个子样本进行抗生素敏感性试验（多西环素， n = 311；四环素， n = 195；红霉素， n = 366；阿奇霉素， n = 350；环丙沙星， n = 350；复方新诺明， n = 195；萘啶酸， n = 140），使用市售抗生素药敏纸片（Oxoid），并遵循临床和实验室标准协会（CLSI）指南 ⁴³ 。2017 年，新增了氨苄西林（ n = 187）、头孢曲松（ n = 171）和头孢克肟（ n = 171）的敏感性试验。对所有抗菌药物均敏感的大肠杆菌 （美国典型培养物保藏中心 25922）用作对照菌株。

使用 Wizard 基因组 DNA 提取试剂盒（Promega），从 37°C 下于 Luria-Bertani 培养基中培养过夜的 5 ml 霍乱弧菌培养物中提取基因组 DNA。通过琼脂糖凝胶电泳确认 DNA 完整性，并使用 NanoDrop 2000 分光光度计（Thermo Fisher Scientific）评估其纯度。在 2014-2018 年的监测研究中，使用 HiSeq 2500 平台（Illumina）在 Wellcome Sanger 研究所进行了 150 bp 双端测序。测序数据已提交至欧洲核苷酸数据库 (ENA)，研究编号为 ERP112767 。在 2% icddr,b 达卡医院监测研究中，使用 NextSeq 2000 进行了双端测序。测序数据已提交至 ENA 数据库，研究编号为 ERP167534 。使用 fastp v.0.23.4 对 FASTQ 文件进行修剪，从 reads 的 5′ 端和 3′ 端滑动窗口，并修剪平均质量低于 20 的碱基。使用 FastQC v.0.11.8 和 MultiQC v.1.8 验证 read 质量。

印度北部样本采集和测序

本研究已获得昌迪加尔医学教育与研究研究生院（PGIMER）伦理委员会的批准。在印度北部昌迪加尔及其邻近邦的疫情调查期间，采集了粪便和水样。此外，还从包括河流和池塘在内的淡水地点采集了水样。粪便和水样分别置于 Cary-Blair 培养基中，并在冷链条件下运送至昌迪加尔医学教育与研究研究生院（PGIMER）进行处理。

水样经 0.22 μm 硝酸纤维素滤膜过滤后，在磷酸盐缓冲液中涡旋振荡以释放黏附细胞。过滤后的水样（100 μl）和粪便样本均在 37 °C 的碱性蛋白胨水中富集培养 6-8 小时。随后将富集后的样本接种于血琼脂和硫代硫酸盐-柠檬酸盐-胆盐-蔗糖琼脂平板上，并在 37 °C 下继续培养 18-24 小时。收集的样本还检测了其他肠道病原体，例如志贺氏菌和沙门氏菌 。使用基质辅助激光解吸/电离飞行时间质谱法鉴定疑似霍乱弧菌或其他肠道病原体的菌落。基因组 DNA 的提取方法如上所述。采用 Illumina 平台进行高通量基因组测序，生成 150 bp 双端测序序列，并按上述方法进行测序数据的质量控制。测序数据已存入 ENA 数据库，研究编号为 ERP188886 和 ERP188887 。

系统发育分析

我们将这些基因组置于一个包含 7PET 基因组序列的全球集合中进行分析（补充数据 1 ）。使用 snippy v.4.6（参考文献 ⁴⁴ ）将 reads 比对至霍乱弧菌 N16961，构建伪基因组比对。使用 snp-sites v.2.5.1（参考文献 ⁴⁵ ）构建仅包含 SNP 的比对，并使用 snp-dists v.0.7（参考文献 ⁴⁶ ）计算基因组间的成对距离。使用 Kraken v.1.1.1（参考文献 ⁴⁷ ）计算归属于霍乱弧菌和 ICP1 的 reads 比例。与参考菌株 N16961 的 SNP 距离大于 400 或归属于霍乱弧菌的 reads 比例低于 90%的样本被排除。使用 rhierBAPS ⁴⁸ 将 7PET 样本分为主要谱系，并根据成对 SNP 距离小于 20 的标准将其分为亚谱系。使用 IQ-TREE v.1.6.12 和 HKY+F+I 替代模型 ⁴⁹ 构建了 (1) 全球 7PET 和 (2) sBD1 以及 (3) 2014-2018 年系统监测研究中的 BD2 的最大似然系统发育树。使用 TreeTime v.0.7.4 (参考文献 ⁵⁰ ) 对系统发育树进行重新定根，并构建最大似然时间尺度系统发育树。使用 TreeTime 的“迁移”功能估计每个节点的离散祖先状态，包括位置、亚谱系以及关键基因和移动遗传元件 (MGE) 的存在情况。利用全球 7PET 系统发育树中祖先位置之间的转换来估计国家间的传播，而利用 sBD1 和 BD2 系统性 2014-2018 年监测研究系统发育树中位置之间的转换来推断孟加拉国境内各行政区之间的传播事件。系统发育树使用 R 包 APE v.5.8（参考文献 ⁵¹ ）和 ggtree v.3.12.0（参考文献 ⁵² ）进行可视化。地图使用 R 包 sf v.1.0-21 生成（参考文献 ⁵³ ）。 ⁵³ ) 使用公开可用的空间矢量数据集，包括 Natural Earth（公共领域；R 包 rnaturalearth v.1.1.0；参考文献 54 ）、Bangladesh R 包 v.1.0.0（参考文献 55 ）、HydroRIVERS（HydroSHEDS 数据库 56 ）和世界银行主要流域数据集 57 。

子抽样

我们采用重复子抽样来验证抽样强度是否会影响推断的区域和国家间传播事件。对于从 2014-2018 年孟加拉国全国监测研究中推断出的区域间传播事件，我们进行了重复子抽样（重复 10 次），以获得预期比例，即假设每年各区域急性水样腹泻病例的比例反映了居住在该区域的孟加拉国人口比例。我们每年都尽可能使用最大样本量来达到此比例。为了增加每个区域内的样本数量，我们排除了迈门辛格，并将拉杰沙希/朗布尔和库尔纳/巴里萨尔分别归入北孟加拉和南孟加拉。对于每个子样本，我们使用 IQ-TREE v.1.6.12 构建了 sBD1 和 sBD2 的系统发育树。我们使用 TreeTime v.0.7.4 对树进行重新定根，并推断每个节点的祖先状态。

对于国际传播事件，我们从霍乱弧菌基因组的全部样本中抽取子样本（重复十次），以确保每年来自孟加拉国、印度和世界其他地区的样本数量相等。我们每年都尽可能多地采集样本以实现这一目标。对于每个子样本，我们使用 IQ-TREE v.1.6.12 构建了 7PET 系统发育树。如上所述，我们使用 TreeTime v.0.7.4 对树进行重新定根，并推断每个节点的祖先状态。

由于核苷酸多样性 (π) 受分析基因组数量的影响，为了比较孟加拉国和印度的核苷酸多样性随时间的变化，我们对两国每年的 15 个样本进行了重复子样本采集（重复 10 次）。使用 Pegas v.1.3（参考文献 ⁵⁸ ）计算了每年每个子样本的核苷酸多样性和 Tajima’s D 值 。

毒力基因、抗生素耐药性和移动遗传元件

使用 SPAdes v.4.1.0 进行霍乱弧菌基因组组装。使用 ARIBA v.2.14.6 (ref. ⁵⁹ ) 检测抗生素耐药基因（针对综合抗生素耐药数据库）、毒力因子（针对毒力因子数据库）、 ctxB 类型、N16961 参考基因组中存在的所有基因以及 ICP1 抗防御基因（ csy1-4 、 cas1 和 cas3 基因的参考基因组为 MW794190.1 ； odn 基因的参考基因组为 MW794192.1 ）的完整阅读框的存在情况。使用 Mash v.2.1.1（参考文献 ⁶⁰ ）筛选包含在 MGE 序列中的样本读取集，包括 PLE1 ( KC152960.1 )、PLE2 ( KC152961.1 )、PLE 3-10（参考文献 ⁴¹ ）、PLE11（使用 Panaroo v.1.3.4 从头基因组组装中鉴定；参考文献 ⁶¹ ）、SXT-ICE ^TET ( MK165649.1 )、SXT-ICE ^GEN /ICEVchInd5 ( KY382507.1 )、ICEVchBan9 ( CP001485 )、ICEVchCHN143 ( KT151654 )、ICEVchInd4 ( GQ463141 ) 和质粒 pCNRVC190243（参考文献 ⁶ ）（ OW443149.1 ）。 ）和 pSAG71 ( CP053818.1 )。如果一个样本与另一个样本共享 1000 个哈希值中的 700 个，则认为该样本存在可移动遗传元件 (MGE)。使用 Kraken v.1.1.1，基于比对到 ICP1 的 reads 比例大于 0.1% 的条件，鉴定出与 ICP1 共测序的样本。使用 R 语言中实现的多变量逻辑回归模型评估关键基因和 MGE 与米汤样便、严重脱水和 K139 流行率之间的关联，并校正所有其他关键基因/MGE 的存在、时间和采样地点。

二次宏基因组分析

我们对来自 40 个霍乱流行国家的宏基因组进行了系统性筛查（补充表 1 ）。我们检索了美国国家生物技术信息中心序列读取档案库（NCBI Sequence Read Archive），寻找已进行 Illumina 全基因组 DNA 测序且分类 ID 分别为 408170（人类肠道宏基因组）、2705415（人类粪便宏基因组）或 749906（肠道宏基因组）的样本，其中宿主为智人（Homo sapiens ）。我们使用 Sylph v.0.8.1（参考文献 ⁶² ）快速筛查霍乱弧菌（V. cholerae） ，该软件使用包含参考基因组 N16961 的自定义数据库。我们使用 Kraken2 v.2.0.8 筛查 ICP1，该软件使用仅包含 RefSeq 病毒的数据库。我们使用 Mash v.2.1.1 检测移动遗传元件（MGEs）。使用 SPAdes v.4.1.0 对比对到 ICP1 的 reads 超过 1000 的样本进行组装，并使用 BLAST v.2.7.1 比对参考序列 GCA_000893175.1 来鉴定与 ICP1 对应的 contig。我们还从 NCBI 下载了 60 个 ICP1 组装结果，并从其中 74 个组装结果中鉴定出 ICP1 contig，这些组装结果中 ICP1 与霍乱弧菌共测序。使用 Prokka v.1.14.5 对参考序列 GCA_000893175.1 覆盖率超过 50% 的组装结果进行注释，并使用 panaroo v.1.3.4 生成核心基因比对。使用 Snp-sites v.2.5.1 创建仅包含 SNP 的比对，并使用 IQ-TREE v.1.6.12（HKY+F+I 替代模型）构建系统发育树。使用 phytools v.2.4-4（参考文献 ⁶³ ）对系统发育树进行中点定根。使用 BLAST v.2.7.1 鉴定基因组组装中特定基因的存在情况。

报告摘要

有关研究设计的更多信息，请参阅本文链接的《 自然投资组合报告摘要》 。