科研| Genome Biology（IF：13.583）：肠道微生物组对BCG诱导的训练免疫的影响

为了研究该微生物组对卡介苗疫苗诱导的训练免疫的影响，我们用卡介苗疫苗接种了321名西欧背景的健康个体[18-75岁，中位年龄23岁，其中有57%为女性，83%的受试者的体重指数（BMI）介于18.5和25之间]。在参与者接种卡介苗疫苗之前，收集粪便和血液样本。接种后在2周和3个月后采集血液，如图1A所示。从321份（见方法）粪便样本中提取DNA，进行扩增，并进行全基因组测序。随后用金黄色葡萄球菌或结核分枝杆菌体外刺激每个时间点从参与者分离的外周血单个核细胞（PBMC）；24h后测定IL-6、IL-1β和TNF-α，7d后测定IFN-γ。此外，在起始点检查时对外周血血浆中的代谢物进行了分析。我们的方法包括：（1）所有微生物物种的组装，（2）测试微生物对细胞因子反应的调节作用，（3）评估免疫调节微生物是否比其他微生物对循环代谢物具有更强的作用，以及（4）确定哪些代谢物与免疫调节微生物物种相关（图1B）。

图1 研究设计和粪便微生物组样本的从头组装分析概述。A研究设计概述：在采集粪便和血液样本进行宏基因组、代谢组和体外细胞因子产生分析后，健康成年志愿者（N=321）在起始点时接种卡介苗疫苗。在起始点、2周后和3个月后检查时，用金黄色葡萄球菌和结核分枝杆菌体外刺激从受试者分离的PBMC。B关联发掘策略：将单个宏基因组物种泛基因组（metagenomic species pangenomes, MSPs）的丰度与2周和3个月后细胞因子表达的倍数变化进行比较。随后，对2种模式中显著的MSP进行相关循环代谢物富集评估，同时对宏基因组组装进行搜索，以寻找解释差异代谢潜能的酶。C微生物组图谱：用MSPminer软件将量化为345个MSP丰度。MSP与门（全部345个）、属（315个）和种（208个）水平上的已知分类群相匹配。条形图显示了MSP的注释水平，按基因数量分为四分位。基于相对MSP丰度之间Jensen-Shannon距离的样本；D树状图显示了3个聚类簇：截止点（虚线）是簇间距离和簇内距离的最大平均差（见图S1）。E基于Jensen-Shannon距离的样本微生物丰度多维标度（MDS）图：样本由簇（C1–C3，如D中所示）着色，簇大小在括号中注明。突出显示了三个丰度最高的门（厚壁菌门、拟杆菌门和放线菌门）。F每个样本检测到的门水平的相对丰度：括号中的百分比报告了队列范围内的丰度。G MSP的数量与问卷变量显著相关，按门分层。颜色编码如F所示。通过Maaslin2广义线性模型（FDR<0.25）评估每个个体变量对每个MSP相对丰度的变化。

基于参考的MetaPhlAn 2.2分析检测到202种微生物，而使用MSPminer重新组装基因组则产生了170多万个微生物基因。345个核心宏基因组物种泛基因组（MSP）与一个门匹配，其中315个（91.3%）和208个（60.3%）分别与属和种水平匹配，更多的基因导致更精确的注释（图1C）。随后在所有下游分析中使用从头组装；除了检索所有202个初始物种外，该方法还使所分析的宏基因组空间的大小增加了1.7倍。在下文中，使用MSP将一个分类单元称为与一门、属或种相匹配的基因集。

使用核心MSP丰度进行的群落水平剖面分析产生了三组不同的样本（图1D；附件1：图S1A）。多维比例图中的平均聚类组成和样本定位反映了三个最常见门的丰度（图1E）。厚壁菌门在簇1中占优势，平均相对丰度为72.6%，其次是拟杆菌门（16.64%，簇2）和放线菌门（9.32%，簇3）（图1F）。簇1还显示疣状结肠炎和Euryarchaeota的相对丰度增加（附件1：图S1B-E）。总体构成反映了发达的、西方类型的群落，这些群落往往以厚壁菌为主，而拟杆菌的存在往往反映了饮食驱动的变化。这一点得到了与473名健康参与者组成的独立荷兰500FG队列的物种水平特征无法区分的证实（附件1：图S1F-G）。

我们研究的所有参与者都被要求完成一份包含一般健康相关问题的问卷调查。185个变量中有18个与三个微生物组相关（P<0.05，未经调整的Cauchy检验），包括饮食（水果摄入量）、性别、吸烟状况、接种史和胃肠道状况相关变量（如：附件1：图S2；附件2：表S1）。

单个物种相对丰度的变化与25个变量相关（MaAsLin2广义线性模型，P_FDR<0.25），其中10个与更多物种相关（图1G；附件2：表S1）。尤其是阴道酵母菌感染、性别和日常运动与大多数微生物种类有关（图1G）。聚类关联遗漏的其他物种特异性变量包括高血压、青年期花粉热和抗生素使用。

我们经过测定受刺激PBMC上清液中细胞因子的产生，以评估训练免疫应答（24小时刺激后的IL-6、IL-1β、TNF-α）和特异性免疫应答（7天刺激后的IFN-γ）。与起始点相比，卡介苗疫苗接种后2周和3个月，金黄色葡萄球菌刺激产生的IL-1β和TNF-α以及结核分枝杆菌刺激产生的IFN-γ逐渐增加（图2A；附件1：图S3A）。训练后的免疫（通过TNF-α、IL-1β、IL-6产生的增加表示）以及特异性（结核分枝杆菌诱导的IFN-γ）和异源淋巴细胞衍生的（金黄色葡萄球菌诱导的IFN-γ）细胞因子反应相互关联（图2B，C）。在体外刺激PBMC后测得的任何细胞因子反应中，微生物组均未显示出显著差异（对于金黄色葡萄球菌IFN-γ，最大P=0.16，附件1：图S3B），这促使我们寻找潜在的个体分类群或功能性免疫调节因子。

图2 微生物物种丰度随着训练和刺激特异性免疫反应的不同而变化。在卡介苗疫苗接种后的起始点（b）、2周（2w）和3个月（3m）用结核分枝杆菌或金黄色葡萄球菌刺激时的体外细胞因子表达。P值通过配对Kruskal-Wallis检验确定。在2周（B）和3个月（C）时，细胞因子表达变化之间存在显著相关性，测量为与起始点水平的倍数变化（Spearman相关性，P<0.05）。对角线细胞显示每个单个细胞因子相对于起始点的显著变化（T检验，*P<0.05，***P<0.001）。D与FDR<0.2且患病率>20%的训练免疫应答相关的显著MSP。每个MSP的效果通过TNF-α、IL-1β和IL-6反应与起始点的倍数变化的线性模型进行估计，并根据年龄和性别进行调整。E与特异性反应（结核分枝杆菌刺激的IFN-γ）相关的显著MSP，FDR<0.2，患病率>20%，使用与训练和特异性反应模型中每个MSP未调整P值的D.F分布相同形式的模型。均匀分布的偏差表明MSP对细胞因子表达的影响更为显著。G每个MSP和时间点对细胞因子折叠变化的个体效应，量化为解释方差乘以效应符号。左侧报告了每个MSP的关联数和物种水平注释的置信度。在随后的2个时间点，H MSP 112（罗斯氏菌属）丰度降低了三种训练免疫表型中的2种。将MSP的相对丰度转换为秩，并通过等频分块将秩分为若干个单元。标题中的P值指经过训练免疫线性模型。中的系数为Spearman 相关性（*P<0.1，***P<0.001）。I MSP 091（瘤胃球菌）和MSP 181（缓慢爱格士氏菌）将丰度与特异性反应表型进行排序。将MSP的相对丰度转换为秩，并通过等频分块将秩分为若干个单元。标题中的P值指具体响应线性模型中的系数为Spearman相关性（*P<0.1，**P<0.01）。

微生物与2种免疫系统模式之间的关联通过2个线性模型进行评估，根据年龄和性别进行调整，测试对训练免疫和特异性T细胞反应的影响。经过训练免疫力与40种重要物种相关（患病率>20%，FDR<0.2，图2D；附件3：表S2），而特异性免疫力与三种物种相关：瘤胃球菌（MSP 091）、Eggerthella lenta（MSP 181）和嗜热链球菌（MSP 076）（图2E；附件3：表S2）。当FDR<0.1时，27个物种与训练免疫相关。由于训练免疫反应的模型结合了比特异性反应（一种细胞因子）更多的生物标记物（三种细胞因子），因此它在相同的显著性阈值下检索到更多的物种。尽管如此，这2种方法都引起了P值的不均匀分布（图2F），支持了微生物组对卡介苗疫苗反应的影响。

接下来，我们检查了疫苗接种后2个时间点的物种相对丰度和个体细胞因子产生之间的相关性（图2G）。罗斯氏菌属基因组（MSP 112）仅在属水平上可识别，在2个时间点与IL-1ß和TNF-α产生的增加呈负相关（P=2.67e−05），FDR=0.0011，图2H）。在特异性反应模型中，瘤胃球菌（MSP 091，P=0.0012，FDR=0.17）和轮状大肠杆菌（MSP 181，P=0.0030，FDR=0.17）在2个时间点分别对IFN-γ表现出一致的负效应和正效应（图2I）。特异性反应模型揭示了单个细胞因子和微生物种类之间的相关性，尽管对于所有细胞因子-刺激物组合，未调整的相关性数量在10到23个微生物之间（如：附件4：表S3）。

细胞因子产生的可检测变化产生了潜在的免疫调节物种，其确切数量随选定的统计阈值而变化。为了证实细胞因子的相关性并缩小潜在的致病机制，我们比较了免疫调节微生物对血浆代谢组的影响。共有1607个质量/电荷（m/z）峰与已知的化合物特征相匹配，并分为20组，针对8个主要分子类别进行了浓缩。

典型相关分析证实了微生物组对化合物强度的强烈影响，因为前25个典型组分（CCs）中的15个显示出显著的解释差异（如：附件1：图S4D-E）。前25个CCs的t-SNE预测显示，按分子类别进行了强有力的分组，在含有羧酸和衍生物、甘油磷脂和有机硫酸的簇中，微生物明显共存（图3A）。三种免疫调节物种罗斯氏菌属（MSP 112）、瘤胃球菌属（MSP 091）和缓慢爱格士氏菌（MSP 181）强烈影响代谢物空间的组成，如图3A中的节点所示。

 图3 免疫调节物种强烈影响血清代谢组。对所有321个样本中345个MSP和1607个注释代谢物的丰度进行典型相关分析（Canonical correlation analysis, CCA）。使用t-SNE在2个维度上预测了25个解释方差显著增加的典型成分。突出显示了三种免疫调节MSP（左）和代谢物类别（右）的位置。点与原点（坐标0，0）的联合偏差反映了代谢物和MSP之间增加的协方差和分组。B对细胞因子反应有影响的MSP对代谢物强度有较大影响。所有代谢物中每个MSP的最大绝对Spearman相关性分布。至少20%的样本中存在的MSP以及2个对数线性模型中P<0.05（左）、<0.1（中）和<0.2（右）的显著MSP以蓝色显示，其余非显著MSP以白色显示。垂直虚线显示了显著相关物种和代谢物的阈值（P_Bonferroni<0.05）。C所有32种显著相关的代谢物均为34 MSPs。彩色热图细胞显示出与P_Bonferroni<0.05的Spearman相关性。列（MSP）按2个对数线性模型中任何一个的最小P值排序。星号表示与细胞因子反应相关的MSP，如图2所示。 

接下来，我们研究了与其他微生物相比，免疫调节微生物是否对个体代谢物具有更强的影响。所有成对Spearman相关性检验显示在34个MSP和32个代谢物强度之间存在198个显著关联（P_Bonferroni<0.05）（如：附件5：表S4）。虽然大多数相关性为阳性（183/198），但由此产生的代谢物强度增加可能是微生物介导的合成代谢或化合物分解所致。为了测试在2个线性模型中具有较大影响的物种是否更有可能对代谢组产生显著影响，我们使用了三种不同的阈值，选择了至少40种物种，并计算了对任何一种代谢产物的最大绝对影响。所有这三种设置都在免疫调节微生物物种中产生了显著的富集效应（最大P<10⁻¹¹，卡方检验，图3B）。这证实了对细胞因子产生影响的微生物种类也对循环代谢物产生强烈影响。

在198个重要的代谢物关联中，18个与罗斯氏菌属（MSP 112）的存在相关，这是所有MSP中最多的，一个与缓慢爱格士氏菌（MSP 181）相关（图3C）。由于注释的不确定性和途径介导的对代谢物丰度的依赖性，我们试图通过使用KEGG途径对代谢物进行分组来确定一般机制。苯丙氨酸代谢受罗斯氏菌属（MSP 112）的影响最为强烈，该途径中有四种代谢产物关联。罗斯氏菌属（MSP 112）还影响色氨酸代谢和次级代谢产物生物合成中的三种代谢产物（图3C）。

非靶向血浆代谢组学检测到2种SCFA，乙酸盐和丙酸盐，它们是膳食纤维代谢的微生物产物和宿主微生物信号的候选介质。这2种分子与训练免疫呈显著正相关，但与特异性反应无关（附件1：图S5）。虽然免疫调节微生物与乙酸盐或丙酸盐均无显著相关性，缓慢爱格士氏菌与丙酸盐的正相关性最高。

丁酸盐，一种先前被证明能减少IL-6产生的SCFA，在循环中未被检测到，可能是由于其在肠道中的吸收。尽管如此，在一个罗斯氏菌属（MSP 209）和5个粪球菌（Coprococcus）物种中检测到丁酸激酶，它催化丁酸生物合成的最后一步，其中粪球菌（MSP 030）与经过训练免疫有最强的正相关性，特别是在接种疫苗后IL-6的增加（如：附件1：图S6A-C）。相反，粪球菌（MSP 213）与训练免疫反应降低相关，表现为3个月后TNF-α表达降低（附件1：图S6D）。免疫调节性罗斯氏菌属（MSP 112）中丁酸激酶的缺失可能源于组装基因组中的缺口，或允许在介导BCG诱导免疫时采用SCFA以外的代谢方式。

来自以下方面的综合证据：（i）微生物对细胞因子产生的模拟效应，（ii）疫苗接种后2个时间点的物种相对丰度和2种经过训练免疫相关细胞因子产生之间的关联，以及（iii）免疫调节微生物种类对血浆代谢组的影响表明，罗斯氏菌属介导的苯丙氨酸代谢影响训练免疫。

对苯丙氨酸代谢物的影响促使我们研究罗斯氏菌属基因组并确认相关酶的存在。我们通过更深入的粪便亚基因组分析，比较了罗斯氏菌属与系统发育和生物化学相关分类群的差异酶和代谢潜力。利用核心基因对所有345个已鉴定的MSP进行了系统发育分化分组（PhyloPhlAn 3.0，附件1：图S7），揭示了一个包含大多数罗斯氏菌属物种的分支（图4A）。系统发育最接近的厚壁菌包括5株瘤胃球菌、3株粪球菌和2株直肠真杆菌MSP。虽然大肠埃希菌在丁酸生产方面与罗斯氏菌属有关，但我们也包括了一株生物学上独特但丰度高的双歧放线杆菌（MSP 111）。组装的小组由59个MSP组成，这些MSP来自85208个组装的双歧杆菌（7个MSP）、罗斯氏菌属（13个MSP）、瘤胃球菌（17个MSP）、真杆菌（15个MSP）和粪球菌（7个MSP）基因，每个属在队列中的患病率都超过95%（图4B）。

图4 罗斯氏菌属编码的苯丙氨酸途径酶是一种普遍而高丰度的酶。A检测到的MSP之间系统发育关系的子树，包含大多数映射到罗斯氏菌属的MSP。B队列中每个属的总丰度。C在样本中组装编码苯丙氨酸代谢途径77种酶中20种酶的基因（KEGG直系同源ko00360），按每个属的MSP分类。热图显示了检测到基因的样本。D苯丙氨酸代谢途径中单个基因的总丰度。对于每个KEGG直系同源群，所有组装同源物的统计数据都进行了汇总。与小组中的其他四个属相比，由罗斯氏菌属编码的每种酶的丰度（每百万映射读取的读取数，RPKM）都有所增加。对于罗斯氏菌属丰度最高的情况，显示了最多四次比较的最小P值。黑点表示中值，误差条表示第25和75个百分位。E小组中MSP的Spearman相关性与苯丙氨酸代谢途径的16个检测到的代谢物特征（质量-电荷、m/z、质谱峰）对应。仅显示顶部匹配的复合名称，而右侧显示不明确的匹配。对于每个属，选择与任何相关MSP相关的最强效应，并标记为显著受试者，P<0.05。与MSP 112（罗斯氏菌属）相关的效果在黑框中突出显示。

共有330个基因与苯丙氨酸途径（KEGG ko00360）中77种酶中的20种酶对应。罗斯氏菌属贡献了最多的基因（92,45.1%平均患病率；图4C；附件6：表S5），而小组中的其他属贡献较少，瘤胃球菌（77,19.2%）、真杆菌（66,36.3%）、粪球菌（50,35.5%）和双歧杆菌（45,46.1%）。罗斯氏菌属还编码了苯丙氨酸酶（20种中有13种）。其中，hisC、mhpE、yhDR、paaH、paaF、paaI和enr在罗斯氏菌属中的含量明显高于其他四个属（图4D）。

共有15个m/z峰与苯丙氨酸途径（KEGG ko00360）中的化合物匹配，由于匹配度不明确，可能跨越多达19个化合物。罗斯氏菌属的物种与途径中15种化合物中的11种显著相关，其中MSP 112对其中6种化合物起作用（Spearman相关性P<0.05，图4E）。瘤胃球菌属、粪球菌属、双歧杆菌属和真杆菌属的关联性较小（图4E）。综上所述，与近缘物种相比，罗斯氏菌属对苯丙氨酸代谢途径的影响增加，这是通过增加酶的潜力和代谢产物的强度来证实的。

最后，我们更仔细地研究了罗斯氏菌属中发现的13种苯丙氨酸代谢酶及其对循环代谢物的不同影响。其中5种酶（hisC、mhpE、paaH、yhDR和paaF）为超过一半的MSP所共有，而其余8种则更具物种特异性（图5A）。aaaT（L-苯丙氨酸N-乙酰转移酶）和paaI（酰基辅酶A硫代酯酶）的共同出现是免疫调节MSP 112的特异性，因为它们在任何其他罗斯氏菌属MSP中都没有发现。这些酶在从苯丙氨酸到N-乙酰-L-苯丙氨酸和苯乙酰谷氨酰胺的途径中起作用，其血浆浓度与罗斯氏菌属（MSP 112）丰度增加有关（图4E）。

图5 罗斯氏菌属对苯丙氨酸代谢途径影响的综合分析（KEGG ko00360）。A罗斯氏菌属中MSP编码的酶的总丰度，黑框中突出显示MSP 112。B观察苯丙氨酸途径中的酶和化合物。蓝色和橙色基团中的化合物和酶相应地着色。罗斯氏菌属编码的酶与苯丙氨酸代谢复合物强度之间存在显著的Spearman相关性。仅显示顶部匹配的复合名称。复合酶组由橙色和蓝色条表示。黑线强调2组之间的分离。C酶的总丰度与B中蓝色和橙色条带分组的化合物的平均强度之间的相关性。D所选化合物与来自罗斯氏菌属基因组的六个（蓝色组）或5个基因（橙色组）的随机组之间的相关性的零分布。观察到的Spearman 相关性（如C）用黑色垂直线标记，零分布用灰色表示。苯丙氨酸代谢途径（KEGG ko00360）。检测到的化合物和酶呈蓝色和橙色。

 单个酶的丰度和化合物强度之间的相关性揭示了2种不同的途径模式。首先，L-酪氨酸、L-苯丙氨酸和苯丙酮酸（m/z 模糊匹配: 2- 羟基 -3- 苯基丙烯酸酯，反 -3- 羟基肉桂酸酯)伴随着5种常见的酶（hisC、mhpE、paaH、yhDR和paaF）以及hipO而增加（图5B中以蓝色条分组）。其次，HPD、paaK、aaaT和padE与苯丙酸盐、马尿酸盐（m/z模糊匹配：苯乙酰甘氨酸）和苯乙酰谷氨酰胺（图5B中以橙色条分组）的增加有关。后者的化合物也随着共有的5种酶的存在而强烈地减少，这表明了2组化合物和酶之间的对立机制。

我们评估了观察到的相关性的显著性，从而验证各组化合物强度的差异确实更可能是由于相关酶所引起的，而非是其他基因所引起。酶丰度和化合物强度之间的 Spearman 相关性显示了比任何单独的复合酶联合更强的影响（图5C，左上和右下）。相反，各组之间的相关性表现出弱到负相关性（图5C，右上和左下）。最后，将观察到的相关性与来自所有22,832个罗斯氏菌属基因的6个 (蓝色组) 或4个 (橙色组) 基因的10,000个随机抽取的相关性的零分布进行了比较。。事实上，与随机的基因组相比，所有四种观察到的效应都是极端的（图5D），证实了相反的化合物-酶关联不太可能是由其他罗斯氏菌属基因引起的。图5E显示了检测到的化合物和酶的全套，并显示了基团分离。这些结果揭示了2种机制，每种机制都以一组相互排斥的酶和化合物为特征，概述了罗斯氏菌属物种水平的差异如何影响循环中的苯丙氨酸代谢物。

一项研究利用在特定培养基上生长的单个微生物种类的上清液进行靶向代谢组学研究，提供了罗斯氏菌（包括了在物种水平上与罗斯氏菌MSP 112的邻近的物种：Roseburia inulinivorans）改变苯丙氨酸代谢途径中化合物的能力的额外证据（图4A）。Roseburia inulinivorans上清液中观察到2种化合物至少有2倍的变化：yhDR和hisC下游编码的乙醇酸增加，而在mphE上游编码的琥珀酸下降（如：附件1：图S8，图5D）。虽然我们在血浆代谢组方法中未检测到这些化合物，但变化方向与相关基因一致，并支持罗斯氏菌属在苯丙氨酸代谢途径中的积极作用。