糖尿病是一种由β细胞功能障碍和/或胰岛素抵抗引起胰岛素的相对或绝对缺乏,从而导致血糖水平高出正常范围的慢性代谢性疾病
[1]。作为最常见的慢性病之一,糖尿病的患病率和患病人数均呈持续增长的态势,给全球造成了沉重的公共卫生负担和经济负担。开展其可变病因研究对于实现疾病的三级预防和优化临床管理至关重要
[2-4]。肠道菌群不仅参与包括胆汁酸、蛋白质、维生素等多种物质在内的生理代谢过程,还从炎症与免疫反应的调节、肠道屏障完整性的维持、能量的提供等多方面对机体产生影响
[5-6]。肠道菌群失调表现为微生物组成和功能的改变
[6],可以通过不同的机制通路引起体内免疫、炎症反应和代谢失调,进而参与糖尿病的病理生理过程
[5, 7-10]。
孟德尔随机化(Mendelian randomization,MR)是一类利用基因型为工具变量(instrumental variables,IVs)来推断表型与结局(疾病)之间因果联系的分析方法
[11-12]。该方法在克服反向因果及混杂因素所带来的偏倚方面具有天然的优势
[13-14]。尽管国内外已有少量研究利用两样本MR方法来探究肠道菌群与不同类型糖尿病间的因果关联,但因存在选定的结局全基因组关联研究(genome-wide association study,GWAS)数据样本量较小、纳入的肠道微生物类别不全等问题,未能为肠道菌群与糖尿病间的关联提供可靠而系统的因果关联证据谱
[15-21]。本研究拟利用两样本MR方法,从不同分类水平系统探索肠道菌群与糖尿病[包括1型糖尿病(type 1 diabetes,T1D)、2型糖尿病(type 2 diabetes,T2D)和妊娠糖尿病(gestational diabetes,GDM)]的潜在因果关联,识别能显著降低或增加糖尿病发生风险的肠道菌群,为糖尿病病因研究提供新证据,同时也为糖尿病的预防和治疗提供新思路。
1 资料与方法
1.1 肠道菌群与糖尿病GWAS数据的选定
本研究所使用的肠道菌群与糖尿病的GWAS数据的基本信息汇总于
表1。
1.1.1 肠道菌群的GWAS数据
肠道菌群的汇总统计数据来源于MiBioGen联盟
[22],该联盟整合了24个队列的18 340例研究对象的16S rRNA基因测序图谱和基因分型数据,对人类常染色体遗传变异和肠道微生物群之间的关联进行了大规模、全基因组的荟萃分析,最终确定了影响肠道微生物类群相对丰度的显著位点
[23]。该联盟仅对至少存在于3个队列、有效样本量达3 000的肠道微生物进行微生物数量性状位点(microbiome quantitative trait loci,mbQTL)分析,最终共分析了211个肠道微生物类群,每个微生物类群都有1个与之对应的GWAS数据。由于211个微生物类群中包含有3个未知科和12个未知属,因此本研究仅纳入余下的196个微生物类群作为暴露,具体包括9个门、16个纲、20个目、32个科和119个属。
1.1.2 糖尿病的GWAS数据
本研究关注的糖尿病主要包括T1D、T2D和GDM。T1D发现样本的汇总数据来源于Chiou等
[24]对9个欧洲队列包含的18 942例T1D患者和501 638例对照进行的GWAS;T1D验证样本的汇总数据来源于Forgetta等
[25]对12个欧洲队列的荟萃分析,包括9 266例T1D患者和15 574例对照。
T2D发现样本的汇总数据来源于糖尿病遗传验证与荟萃分析(Diabetes Genetics Replication And Meta-analysis,DIAGRAM)联盟对已有的32个GWAS进行的一项荟萃分析
[26],共包含74 124例T2D患者和824 006例对照,研究对象均为欧洲人种;T2D验证样本的汇总数据来源于FinnGen联盟于2023年12月最新发布的GWAS数据集,研究对象包括65 085例T2D患者和335 112例对照
[27]。
GDM发现样本的汇总数据同样来源于FinnGen联盟于2023年12月最新发布的GWAS数据集,研究对象包括14 718例GDM患者和215 592例对照
[27]。GDM验证样本的汇总数据来源于Jiang等
[28]运用基于广义线性混合模型的全基因组关联方法对英国生物银行数据库中的2 989个二分类表型进行的二次GWAS。GDM的GWAS数据共包含864例患者和6 977例对照。与线性混合模型相比,该模型具有更好的统计性能。
1.2 IVs的筛选
以P<1×10-5为条件在肠道菌群的GWAS数据中筛选与暴露强相关的单核苷酸多态性位点(nucleotide polymorphisms,SNP)。使用PLINK clumping方法去除连锁不平衡,为保证各SNP之间的相互独立,具体参数设定为物理距离=10 000 kb,r2<0.001。根据前2步筛选得到的SNP位点,从肠道菌群与糖尿病的GWAS数据中分别提取效应等位基因、参考等位基因、效应等位基因频率、效应值、效应值的标准差和P值等数据列信息。对于无法在结局GWAS数据中找到的SNP,使用LDLink寻找目标SNP的代理SNP,即欧洲人种的参考基因组提示连锁不平衡的 r2>0.8的SNP。计算每个SNP的F值,剔除F<10的SNP以避免潜在的弱IVs偏倚。最后,对齐所有IVs的效应等位基因,并进一步剔除回文SNP(相同的等位基因出现在正负链)及与结局强相关的SNP(P<5×10-8)。
1.3 统计学处理
1.3.1 正向MR分析
在正向MR分析中,以人体肠道菌群作为暴露,糖尿病作为结局。采用4种不同的效应估计值模型进行因果推断,包括逆方差加权法(inverse-variance weighted,IVW)
[29]、MR-Egger回归法
[30]、加权中位数法(weighted median)
[31]和加权模式法(weighted mode)
[32]。其中,将IVW作为评估因果效应的主要分析模型,其他3种模型作为IVW的补充和参考,用于评估其结果的稳健性。当IVW模型分析结果提示因果效应显著且其他3种MR模型分析的因果效应估计值方向与其一致时,方能确定因果关联的存在。
1.3.2 敏感性分析
对于MR分析提示因果关联具有统计学意义的“肠道菌群-糖尿病”对,分别采用Cochran’s Q检验和MR-Egger截距检验以评估是否存在异质性和多效性。多效性的存在违背了MR分析的排他性假设,因此应排除存在多效性的因果关联。此外,为了防止肠道菌群与糖尿病之间的潜在混杂因素对因果效应估计的影响,研究进一步利用PhenoScanner(人类基因型-表型的数据库)查找每个IV的潜在相关表型。若IVs中存在与糖尿病或潜在混杂因素强相关的SNP(P<5×10-8),则说明其违反了MR分析的核心假设,予以剔除后再重新进行MR因果效应的估计。
1.3.3 反向MR分析
将正向MR分析中发现的对T1D、T2D或GDM有因果效应的肠道细菌类别再进行反向MR分析(以糖尿病作为暴露,以人体肠道菌群作为结局),以排除糖尿病对肠道菌群有因果效应的可能性。反向MR分析的研究样本、方法原理、分析步骤均与正向MR分析相同。
1.3.4 显著性估计值的验证
初步确定对结局(T1D、T2D和GDM)有因果效应的肠道细菌类别后,将这些肠道菌群类别的GWAS数据与验证样本结合进行二次MR分析,以IVW模型所得的因果效应估计值为参考,观察能否得到相同的结果。
1.3.5 统计软件
本研究所涉及的统计分析和图形可视化均通过R软件(版本4.2.3)实现,所使用的R包有TwoSampleMR、LDlink、MungeSumstats、forestploter等。
2 结 果
2.1 IVs选择
筛选得到2 179、2 176、2 166个SNPs,分别用于肠道菌群与T1D、T2D、GDM的MR分析。各IVs在肠道菌群不同分类水平(包括9个门、16个纲、20个目、32个科、119个属)的分布和
F值见
表2。经统计,所有IVs的
F值均大于10,表明不存在弱IVs偏倚。
2.2 肠道菌群与T1D的关联
IVW模型分析结果(
表3)显示:
Eggerthella属(
OR=1.15,95%
CI 1.03~1.27,
P=0.010)、
Lachnospira属(
OR=1.41,95%
CI 1.01~1.96,
P=0.043)、
Ruminococcaceae UCG004属(
OR=1.19,95%
CI 1.04~1.37,
P=0.013)与T1D存在正向因果关联;Bacillales目(
OR=0.90,95%
CI 0.83~0.99,
P=0.024)、Desulfovibrionales目(
OR=0.82,95%
CI 0.68~0.99,
P=0.043)、Victivallaceae科(
OR=0.91,95%
CI 0.84~0.99,
P=0.033)、
Parasutterella属(
OR=0.88,95%
CI 0.78~0.99,
P=0.035)、
Turicibacter属(
OR=0.87,95%
CI 0.76~1.00,
P=0.048)与T1D存在负向因果关联。除
Ruminococcaceae UCG004属外,其他7种肠道细菌的另3种MR模型分析的因果效应估计值方向与IVW模型一致。
在正向MR分析中确定与T1D存在因果关联的7种肠道微生物的异质性与多效性检验结果见
附表1(
https://doi.org/10.57760/sciencedb.xbyxb.00079)。Cochran’s
Q检验结果提示各肠道微生物的IVs间不存在异质性;MR-Egger截距检验结果提示Victivallaceae科的IVs与T1D间可能存在多效性(egger_intercept=0.058,
P=0.047),予以剔除。每个IV的潜在相关表型见
附表2(
https://doi.org/10.57760/sciencedb.xbyxb.00079)。暂未发现与T1D或潜在混杂因素相关的SNPs,因此无需进行剔除后的MR分析。反向MR分析结果(
附表3,
https://doi.org/10.57760/sciencedb.xbyxb.00079)显示T1D与发现的6种肠道菌群间均不存在因果效应。
验证阶段的IVW模型分析结果(
图1)显示:Desulfovibrionales目(
OR=0.70,95%
CI 0.52~0.93,
P=0.015)与T1D存在负向因果关联;
Eggerthella属(
OR=1.22,95%
CI 1.03~1.45,
P=0.020)与T1D存在正向因果关联。这与发现阶段的分析结果一致。
2.3 肠道菌群与T2D的关联
IVW模型分析结果(
表4)显示:Verrucomicrobia门(
OR=1.10,95%
CI 1.01~1.20,
P=0.037)、Deltaproteobacteria纲(
OR=1.10,95%
CI 1.04~1.17,
P=0.002)、Actinomycetales目(
OR=1.08,95%
CI 1.01~1.16,
P=0.030)、Desulfovibrionales目(
OR=1.09,95%
CI 1.02~1.16,
P=0.010)、Actinomycetaceae科(
OR=1.08,95%
CI 1.01~1.16,
P=0.030)、Desulfovibrionaceae科(
OR=1.09,95%
CI 1.02~1.16,
P=0.016)、
Actinomyces属(
OR=1.13,95%
CI 1.06~1.19,
P<0.001)、
Gordonibacter属(
OR=1.04,95%
CI 1.00~1.07,
P=0.047)与T2D存在正向因果关联;Bacillales目(
OR=0.95,95%
CI 0.90~1.00,
P=0.037),Alcaligenaceae科(
OR=0.90,95%
CI 0.82~0.99,
P=0.035),
Lachnospiraceae NC2004 group属(
OR=0.93,95%
CI 0.88~0.99,
P=0.014)与T2D存在负向因果关联。除Bacillales目和
Gordonibacter属外,其他9种肠道细菌的另3种MR模型分析的因果效应估计值方向与IVW模型一致。
在正向MR分析中确定与T2D存在因果关联的9种肠道微生物的异质性与多效性检验结果见
附表4(
https://doi.org/10.57760/sciencedb.xbyxb.00079)。Cochran’s
Q检验结果提示Verrucomicrobia门的IVs间存在显著异质性(
PMR-Egger=0.002,
PIVW=0.004);MR-Egger截距检验结果提示各肠道微生物的IVs与T2D间不存在多效性。每个IV的潜在相关表型见
附表5(
https://doi.org/10.57760/sciencedb.xbyxb.00079),rs2715439位点与臀围、体重强相关(
P<5×10
-8),2个表型可能是肠道菌群与T2D因果关系链上的潜在混杂因素,剔除该位点后重新进行对应肠道细菌(
Actinomyces属)与T2D的MR分析。MR分析结果显示:在剔除具有潜在多效性的IV后,
Actinomyces属与T2D间的因果关联仍具有显著性(
OR=1.12, 95%
CI 1.05~1.19,
P<0.001;附表6,
https://doi. org/10.57760/sciencedb.xbyxb.00079)。反向MR分析结果(附表7,
https://doi.org/10.57760/sciencedb.xbyxb. 00079)显示T2D与发现的9种肠道菌群间均不存在因果效应。
验证阶段的IVW模型分析结果(
图2)显示:
Actinomyces属(
OR=1.09,95%
CI 1.02~1.15,
P=0.007)与T2D存在正向因果关联。这与发现阶段的分析结果一致。
2.4 肠道菌群与GDM的关联
IVW模型分析结果(
表5)显示:Betaproteobacteria纲(
OR=1.19,95%
CI 1.00~1.40,
P=0.047)、Deltaproteobacteria纲(
OR=1.14,95%
CI 1.00~1.29,
P=0.046)、
Coprobacter属(
OR=1.10,95%
CI 1.01~1.21,
P=0.032)、
Ruminococcus2属(
OR=1.11,95%
CI 1.01~1.23,
P=0.038)与GDM存在正向因果关联;Euryarchaeota门(
OR=0.89,95%
CI 0.83~0.95,
P=0.001)、Tenericutes门(
OR=0.86,95%
CI 0.75~0.98,
P=0.025)、Clostridia纲(
OR=0.82,95%
CI 0.67~1.00,
P=0.046)、Methanobacteria纲(
OR=0.89,95%
CI 0.81~0.98,
P=0.015)、Mollicutes纲(
OR=0.86,95%
CI 0.75~0.98,
P=0.025)、Clostridiales目(
OR=0.83,95%
CI 0.69~1.00,
P=0.048)、Methanobacteriales目(
OR=0.89,95%
CI 0.81~0.98,
P=0.015)、Methanobacteriaceae科(
OR=0.89,95%
CI 0.81~0.98,
P=0.015)、
Methanobrevibacter属(
OR=0.84,95%
CI 0.75~0.93,
P=0.001)与GDM存在负向因果关联。除Deltaproteobacteria纲和Clostridiales目外,其他11种肠道细菌的另3种MR模型分析的因果效应估计值方向与IVW模型一致。
在正向MR分析中确定与T2D存在因果关联的11种肠道微生物的异质性与多效性检验结果见附表8(
https://doi.org/10.57760/sciencedb.xbyxb.00079)。Cochran’s
Q检验结果提示Clostridia纲的IVs间存在显著异质性(
PMR-Egger=0.043,
PIVW=0.022);MR-Egger截距检验结果提示Euryarchaeota门的IVs与GDM间可能存在多效性(egger_intercept=0.049,
P=0.040),予以剔除。每个IV的潜在相关表型见附表9 (
https://doi.org/10.57760/sciencedb.xbyxb.00079),暂未发现与GDM或潜在混杂因素相关的SNPs,因此无需进行剔除后的MR分析。反向MR分析结果(附表10,
https://doi.org/10.57760/sciencedb.xbyxb.00079)显示GDM与发现的10种肠道菌群间均不存在因果效应。
验证阶段的IVW模型分析结果(
图3)显示:
Coprobacter属(
OR=0.63,95%
CI 0.42~0.95,
P=0.029)与GDM存在负向因果关联。这与发现阶段的分析结果矛盾。
3 讨 论
本研究基于两样本MR研究设计,利用现有可获得的样本量最大的肠道菌群与糖尿病的汇总统计数据,从门、纲、目、科、属5个不同的分类水平系统探究了人类肠道菌群与T1D、T2D和GDM的潜在因果关联,从而为暴露与疾病结局提供了可靠且较全面的关联研究证据。
与其他流行病学方法相似,MR方法同样依赖于特定的假设,只有使用符合MR核心假设(具体包括关联性假设、独立性假设和排他性假设)的IVs进行分析方能得到正确的因果效应估计值
[33]。本MR研究通过设定严格的IVs筛选方案以确保IVs的有效性和结果的可靠性。首先,与大多数MR研究设定的相关性条件
P<5×10
-8不同,本研究采用更为宽松的
P< 1×10
-5作为条件筛选与肠道菌群强相关的SNPs位点。这样可解释更多的肠道微生物组特征变异,且已被广泛应用于不同的肠道菌群MR研究中
[15, 23];同时使用PLINK clumping方法去除连锁不平衡,从而保证IVs的强相关性与独立性,符合关联性假设。其次,本研究利用PhenoScanner查询每一个遗传位点的潜在相关表型,对与糖尿病或潜在混杂因素强相关的SNP予以剔除后再进行MR分析,以确保各IVs与混杂因素独立,符合独立性假设;尽管目前排他性假设无法被直接评估和严格验证,本研究采用特定的敏感性分析方法即MR-Egger截距检验以检测多效性存在与否,以
PMR-Egger>0.05作为因果关联发现的必要条件,从而在一定程度上排除多效性对结果的可能影响。
综合系列分析结果,共发现6种肠道微生物与T1D存在因果关联:
Eggerthella属和
Lachnospira属与T1D发生风险增加有关,而Bacillales目、Desulfovibrionales目、
Parasutterella属和
Turicibacter属与T1D发生风险降低有关。其中,
Parasutterella属与T1D的因果关联与既往同类MR研究报道的相一致
[34]。此外,
Eggerthella属、Desulfovibrionales目与T1D的关联研究结果在独立的验证样本中被完全复现。与后者相比,前者即
Eggerthella属在绝大部分自身免疫病患者中的丰度均显著高于健康人群
[35-38]。此外,Alexander等
[39]通过小鼠模型证实
Eggerthella lenta能够通过细胞和抗原非依赖性机制解除对辅助性T细胞17(T helper cell 17,Th17)转录因子Rorγt的抑制,从而诱导肠道Th17激活,最终加重自身免疫疾病模型小鼠的病情。由此可见,
Eggerthella属可能通过一定的免疫机制引起机体免疫稳态的破坏,从而导致免疫相关疾病(包括T1D)发生风险的增加。
本研究发现9种肠道微生物与T2D存在因果关联:Verrucomicrobia门、Deltaproteobacteria纲、Actinomycetales目、Desulfovibrionales目、Actinomycetaceae科、Desulfovibrionaceae科和
Actinomyces属与T2D发生风险增加有关,而Alcaligenaceae科和
Lachnospiraceae NC2004 group属与T2D发生风险降低有关。其中,
Actinomyces属与T2D的关联研究结果在独立的验证样本中被完全复现。
Actinomyces属是一种机会致病菌,在本MR研究中,无论是结合发现样本信息还是验证样本信息,均确证了其与T2D风险间的正向因果关联,这与已有的同类MR研究结果
[20, 40-41]一致。有研究
[42]报道
Actinomyces属在T2D组中的丰度显著高于健康对照组,且与血浆中参与炎症反应的肿瘤坏死因子(tumor necrosis factor,TNF-α)浓度呈正相关,提示其可能与机体内的炎症状态有关。
此外,基于MR分析,本研究发现10种肠道微生物与GDM存在因果关联:Betaproteobacteria纲、
Coprobacter属、
Ruminococcus2属与GDM发生风险增加有关,而Tenericutes门、Clostridia纲、Methanobacteria纲、Mollicutes纲、Methanobacteriales目、Methanobacteriaceae科、
Methanobrevibacter属与GDM发生风险降低有关。在独立的GDM验证样本中发现
Coprobacter属与GDM的关联方向与发现样本中的相反,这可能是用于外部验证的结局GWAS数据样本量(
n=7 841)过小所致,而目前已发表的相关文献报道的结论与从发现样本得出的结论相符。既往有研究
[43]利用两样本MR方法从属的水平探究了肠道菌群与GDM的因果关系,其报道的结果与本研究呈部分一致性,包括:
Coprobacter属、
Ruminococcus2属与GDM发生风险呈正相关;
Methanobrevibacter属与GDM发生风险呈负相关。此外,动物研究
[44]发现在肥胖小鼠模型中
Coprobacter属与空腹血糖、稳态模型评估胰岛素抵抗指数(homeostasis model assessment of insulin resistance,HOMA-IR)呈正相关,提示该菌属对糖代谢稳态的潜在有害效应。相较于T1D和T2D来说,目前有关肠道菌群与GDM的关联研究较为有限,较少有人群研究和动物实验报道特定肠道微生物类别与GDM的关联或解析其在疾病发生和发展中的潜在作用机制。尽管本MR研究为肠道微生物与GDM的关系提供了初步的关联研究线索,但仍需更多的研究予以证实并进一步解释其背后潜在的生物学意义。
本研究验证阶段中仅有少部分肠道菌群与糖尿病的关联研究结果予以完全复现,而其他大部分肠道菌群与糖尿病的因果关联未呈现显著意义,尽管它们大多数与糖尿病发现样本和验证样本分别所得的因果效应估计值具有方向上的一致性。样本量是MR研究中统计功效的重要决定因素之一,较大的样本量能够获得较高精度的因果估计值
[45-47]。研究
[13, 48]表明:导致遗传流行病学不可复制性的重要因素之一是统计功效不足(通常反映有限的样本量)。本研究为了保证验证样本的独立性,用于发现和验证的结局GWAS数据的样本来源不相重叠,使得验证阶段的汇总数据的样本量远小于发现阶段。因此,发现阶段中显著性结果在验证阶段的不可复现性可能归因于有限的样本量所导致的统计功效不足。未来随着全基因组层面上多中心、大样本的基因与各表型的关联研究的深入开展,基于MR研究设计的肠道菌群与糖尿病的因果关系探索需要在更大样本量、更高质量的GWAS数据中进行进一步确证。
本研究仍存在一定局限性:1)尽管结局GWAS数据的样本均来源于欧洲人种,而肠道菌群汇总数据的样本来源是多人种的(欧洲人种占比80%左右),因此可能存在潜在的人群分层偏倚;不同种族人群遗传背景的差异性也限制了本MR研究所得的因果关联证据推广至其他种族人群,即非欧洲人种。2)由于16S rRNA基因测序的分辨率有限,未能从更特定的种的水平来评估肠道菌群与糖尿病的因果关联。3)虽利用生物统计学的方法系统解析了不同分类水平的肠道菌群与各类型糖尿病的潜在因果关联,但是未能从生物功能的角度解释这些关联背后的确切机制,未来还需更多的人群研究和机制实验予以进一步阐明。
湖南省杰出青年基金(2022JJ10087┫。This work was supported by the National Natural Science Foundation ┣82073653)
Hunan Provincial Outstanding Youth Fund(2022JJ10087)