骨质疏松症(osteoporosis, OP)是一种骨骼疾病，其特征是骨密度降低、骨微结构受损、骨强度受损和骨折风险增加[1]。在英国，每年大约有549 000例新的脆性骨折是由OP引起的，其中包括105 000例髋部骨折和86 000例椎体骨折[2]。最新的流行病学数据表示，目前我国约有9 000万人患有OP[3]。

近年来，随着“骨免疫学”概念的提出，免疫系统已被证实与OP密切相关，从免疫角度维持骨形成与骨吸收之间的平衡是维持骨稳态的重要策略之一，也为OP的药物研发提供了新的思路[4-5]。“骨免疫学”一词是由Aaron和Choi在2000年提出的，突出强调了骨骼和免疫系统之间的双向作用，特别是在各种炎症性骨骼疾病中[6]。因此，正确理解炎症与骨骼系统之间的关系就显得至关重要。炎症是OP发生发展的核心机制，炎症通过促进骨吸收和抑制骨形成在OP中发挥作用[7-9]。炎症反应的主要介导分子包括肿瘤坏死因子(TNF)、白细胞介素(IL)和干扰素(IFN)等[10]。这些分子通过影响各种免疫细胞和信号通路，打破骨形成与骨吸收之间的微妙平衡，导致异常的骨重建，从而可能减缓或加速OP的进展[11]。此前的动物实验和临床试验发现并证实了在去卵巢啮齿动物的血液和骨髓中以及绝经后女性体内，发现炎症介质(包括IL-1β、IL-6和TNF-α)的存在增加[12-13]。同时，绝经后女性的CD19+B淋巴细胞数量较低，而这些细胞能分泌更多的粒细胞巨噬细胞集落刺激因子(GM-CSF)[14]。总之，OP患者会表现出慢性低度炎症以及细胞因子和免疫细胞谱的改变。然而，大多数关于炎症介质与OP关系的研究均为观察性研究。这些研究易受偏倚、混杂因素以及潜在逆因果关系的影响。因此，需要更好的证据来了解炎症因子在OP中的作用。

孟德尔随机化(mendelian randomization, MR)是一种基于孟德尔独立分布规律的流行病学病因推断分析方法。MR基于孟德尔第二遗传定律，利用遗传变异作为工具变量来识别暴露因素与疾病之间的因果关系[15]。其独立分配法则表明遗传变异是独立的，不会受到其他混杂因素的影响。因此，本研究汇总了两样本MR分析结果，以评估炎症因子与OP之间的因果关系，并进一步通过药物预测和分子对接研究，验证潜在药物靶点的药理活性，同时拓宽了预测药物的临床应用范围。为进一步研究骨免疫学提供新的生物诊断标志物、治疗策略和理论基础。

1、资料与方法

1.1 研究设计

本研究将41个炎症因子作为暴露因素，与这41个炎症因子显著相关的单核苷酸多态性(single nucleotide polymorphism, SNPs)作为工具变量，OP作为结局变量。同时，双能X线吸收测定法检测骨密度(bone mineral density, BMD)是诊断骨质疏松的主要定量依据，被称为诊断OP的“金标准”,故将BMD(全年龄段和各年龄段)也纳入结局变量。使用两样本MR分析的方法对41个炎症因子与OP和BMD进行因果关联分析，利用Cochran's Q检验进行异质性检验，MR Egger法评估其水平多效性，最后对结果进行敏感性分析来确定结果的可靠性。

为了筛选出可靠的工具变量，得到可靠的结果，本研究进行的MR分析需要满足以下3个关键假设：①最终纳入的SNPs需要与炎症因子之间存在密切的关联性；②SNPs要独立于41个炎症因子与OP和BMD之间的所有混杂因素；③SNPs仅通过41个炎症因子影响OP和BMD。见图1。

图1孟德尔随机化核心假设示意图

1.2 数据来源

作为暴露因素的41个炎症因子的数据来自于Ahola-Olli等[16]的一项研究，其中包括8 293例参与者。OP的GWAS汇总统计的数据来自芬兰数据库,其中包括3 203例患病人群和209 575例对照人群。全年龄和各年龄段BMD的GWAS数据全来自于Medina等[17]的一项研究，其中包括66 628例参与者。本研究关于OP和BMD的GWAS数据都是从公共网站数据库中获取，并且全部为欧洲人群。OP的ID为：finn-b-M13_OSTEOPOROSIS,BMD的ID分别为：ebi-a-GCST005345(0～15岁)、ebi-a-GCST005344(15～30岁)、ebi-a-GCST005346(30～45岁)、ebi-a-GCST005350(45～60岁)、ebi-a-GCST005349(60岁以上)、ebi-a-GCST005348(全年龄段)。

1.3 工具变量

首先，为了在标准范围内筛选出更多与炎症因子相关的SNPs, 本研究以P<5.0×10-6为筛选条件从GWAS数据库中筛选出与暴露显著相关的SNPs作为工具变量。随后，为排除连锁不平衡(linkage disequilibrium, LD)的干扰，本研究将其参数设置为R2<0.001,等位基因距离为10 000 kb, 以确保工具变量的独立性。同时，为满足独立性假设，使用了PhenoScanner v2数据库,将筛选出的SNPs进行逐个检索，去除与混杂因素相关的SNPs[18]。此外，通过计算F统计值评估工具变量的强度，使用了公式F=[(N-K-1)/K]×R2/(1-R2)计算F统计值，其中，N是样本数量，K为工具变量的数量，R2是每个工具变量所解释的暴露变异所占的比例，计算公式为R2=2×EAF×(1-EAF)×(β/SD)2[19]。若F>10表示不存在弱工具变量偏差，MR分析结果可以免受弱工具变量偏差的影响[20]。

1.4 孟德尔随机化分析

本研究采用了5种不同的孟德尔随机化分析方法来评估41个炎症因子对OP和BMD的影响，其中主要的3种分析方法分别为逆方差加权法(inverse variance weighted, IVW)、MR-Egger回归法(MR-egger regression, MER)和加权中位数法(weighted median, WME)。假设所有SNPs都是有效的工具变量时，IVW方法能得出最准确的效应估计值[21]。MR-Egger方法可以检测多效性，但其提供的评估精度较低[22]。WME方法可以在提供至少一半的工具变量都是有效的前提下，给出准确的评估结果[23]。综上，因为IVW方法比MR-Egger方法和WME方法检验的效能更高，所以本研究选用IVW方法作为最主要的分析方法。

1.5 敏感性分析

为了评估41个炎症因子对OP和BMD的因果效应的稳健性，本研究进行了一系列敏感性分析。首先，使用留一法逐个剔除工具变量来评估是否存在某一个SNP驱动因果关联。其次，使用MR-Egger回归法和漏斗图评估工具变量与其他潜在混杂因素之间的多效性，以保证选择的工具变量不会通过除暴露因素以外的因素影响结局。另外，使用了Cochran’s Q统计法来评估SNPs之间的异质性，在该统计法中，若P>0.05则表示没有异质性。最后，还提供了更加直观的散点图来展示MR分析结果。

1.6 候选药物预测

评估蛋白质与药物的相互作用对于理解目标蛋白是否可以用作实际药物靶点至关重要。本研究将使用DSigDB数据库来完成此任务[24]。DSigDB拥有22 527个基因集和17 389个不同的化合物，同时涵盖了19 531个基因，是一个庞大的数据库，它能够将药物和其他化学物质与其靶基因联系起来。具体而言，将MR结果中确定的炎症因子的基因形式结果上传至DSigDB中，预测出能够与其结合的候选药物，从而达到定向靶基因治疗的目的。

1.7 分子对接

为进一步理解药物候选物对药物靶蛋白的影响以及靶蛋白的药物可用性，本研究在原子水平上进行了分子对接，以评估药物候选物与其靶标之间的结合能和相互作用模式。分子对接模拟能够分析配体与药物靶标之间的结合亲和力和相互作用方式。通过识别具有高结合亲和力和良好相互作用模式的配体，可以优先考虑将药物靶标用于进一步的实验验证，并优化潜在药物候选物的设计。本研究使用计算机辅助的蛋白质配体对接软件Autodock Vina 1.2.2[25],将筛选出的前5位重要药物与其对应的2位靶蛋白进行分子对接，因Sphingosine 1-phosphate和sulfasalazine在PubChem中无3D结构，所以最终有3种药物与其靶点进行了分子对接。药物结构数据来源于PubChem Compound Database[26],蛋白质结构数据来源于AlphaFold数据库。首先，从蛋白质和配体文件中移除了所有水分子，并添加了极性氢原子。网格框架的中心被设定为包围每种蛋白质的结构域，并且允许无限制的分子运动。采用了网格点距离为0.05纳米和口袋尺寸为30 Å×30 Å×30 Å来建立对接口袋。整个分子对接过程通过Autodock Vina 1.2.2在模型中进行了可视化。

1.8 统计学分析

本研究使用了R4.3.1软件中的TwoSample MR包进行MR分析。所得结果均以OR值以及相对应的95%置信区间(confidence interval,CI)来解释41个炎症因子对OP和BMD的影响，P<0.05表示差异具有统计学意义。[LL]

2、结果

2.1 孟德尔随机化分析结果

首先，本研究确定了与41个炎症因子相关的SNPs, 经过一系列的质控步骤，最终选出强相关的SNPs进行孟德尔随机化研究。结局为OP的IVW结果表明，TRAIL(OR=1.15,95%CI:1.03～1.29,P=0.012),MCP-1/MCAF(OR=1.22,95%CI:1.04～1.43,P=0.014)与OP发生风险呈正相关。结局为BMD的IVW结果表明，TNF-β(β=-0.076,95%CI:0.87～0.98,P=0.011)与0～15岁人群的BMD呈负相关；IL-7(β=-0.099,95%CI:0.82～0.99,P=0.043)与15～30岁人群的BMD呈负相关；HGF(β=-0.131,95%CI:0.77～0.99,P=0.040)与30～45岁人群的BMD呈负相关；MIP-1α(β=-0.102,95%CI:0.84～0.97,P=0.006)、M-CSF(β=-0.058, 95%CI:0.89～0.99,P=0.048)与60岁以上人群的BMD呈负相关；MIP-1α(β=-0.056,95%CI:0.90～0.99,P=0.025)与全年龄段人群的BMD呈负相关，详见图2。其次，本研究还用了MR-Rgger回归法和Weighted median法对上述8个炎症因子与OP和BMD之间的因果关系进行了分析，详见表1。

图2与OP和BMD有因果关联的炎症因子效应值散点图

表1与OP和BMD具有显著相关性的炎症因子MR效应值

2.2 敏感性分析

Cochran's Q检验和MR-Egger回归法的结果均表明本研究不存在显著的异质性和水平多效性，详见表2。留一法的结果表明，在逐个剔除SNP后对结局无影响，表明本研究具有稳健性。同时漏斗图结果是逐个以SNP作为工具变量分析因果关联效应，其结果显示因果关联效应的散点基本呈对称分布，表明本研究不存在水平多效性。

表2敏感性分析

2.3 候选药物预测结果

本研究利用DSigDB数据库进行潜在干预药物的预测。结果表明，Sphingosine 1-phosphate和trolox是两种最显著的药物，与HGF和TNF相关。大部分药物(例如Sphingosine 1-phosphate, trolox, (+)-chelidonine和sulfasalazine等)参与免疫系统调节，表现出抗炎抗氧化作用(表3)。鉴于慢性炎症和氧化应激是导致OP的三大病理机制，因此这些药物的抗炎抗氧化和激素调节作用对于OP的治疗至关重要。该研究结果不仅预测了药物的潜在靶点，展示了靶点的药物性质，还拓宽了已有药物的临床应用范围，为未来这些药物在临床治疗OP中提供了理论支持。

表3DSigDB预测候选药物

2.4 分子对接结果

为了评估候选药物与其靶点的亲和力，并由此了解药物靶点的成药性，本研究进行了分子对接。使用Autodock Vina v.1.2.2获得前5种候选药物与相应的靶蛋白的结合位点和相互作用，并生成每次相互作用的结合能，因Sphingosine 1-phosphate和sulfasalazine在PubChem中无3D结构，所以最终总共产生6种蛋白质与药物的有效对接结果(见表4和图3)。每一种候选药物都与它的靶蛋白形成可见的氢键和强大的静电相互作用。在该部分研究中，发现HGF与(+)-Chelidonine的结合能最低(-7.747 kcal/mol),结合非常稳定。

表4分子对接结果

图3分子对接结果图

3、讨论

本研究通过孟德尔随机化分析方法全面评估了41种炎症因子对OP和各年龄段BMD的因果效应，结果表明某些炎症因子确实与OP和BMD之间存在因果关联，这些结果从基因遗传角度提示炎症可能是OP和BMD的重要影响因素，同时也意味着炎症在骨质疏松症发病过程中发挥着重要作用。此外，本研究对与这些炎症因子相对应的药物进行了预测和分子对接，进一步证明了这些炎症因子作为治疗靶点的药用价值。

目前，多数研究已经表明，炎症因子调控着OP的发生发展过程，这些研究与本研究的结果是一致的[27-28]。最近的研究表明，骨和免疫反应之间存在串扰和通讯[29]。同时临床数据表明，绝经后妇女体内炎症细胞因子水平升高，包括IL-1β、IL-6和TNF-α等[30-31]。此外，另一项研究发现，在绝经后妇女体内能够检测出最高的TRAIL[32]。T细胞在PMOP患者中也更有可能表达TNF-α[33]。MCP-1通过影响破骨细胞的发育和激活，在骨重塑过程中发挥着重要作用，MCP-1的血清水平与骨密度呈负相关，是早期发现和预防OP的重要生物标志物[34]。此外，MIP-1α通过阻断ERK激活的Dickkopf相关蛋白1上调介导的β-catenin活性来抑制BMSCs的成骨分化潜力。通过中和抗体在体内阻断MIP-1α可改善老年小鼠的小梁骨丢失和骨髓肥胖[35]。除了更高水平的炎症介质外，绝经后妇女的免疫细胞数量也发生了改变。根据免疫学临床研究数据表明，在绝经后女性患者中，CD19+B淋巴细胞数量较低，而这些细胞能分泌更多的粒细胞巨噬细胞集落刺激因子(GM-CSF)[36]。除此之外，临床前数据也证实了在去卵巢啮齿动物的血液和骨髓中，发现炎症介质(包括IL-1β、IL-6和TNF-α)的存在增加[37-38]。除老年女性绝经外，老年男性和女性都会由于增龄导致骨重建失衡，骨吸收与骨形成的比值升高，最终导致进行性骨丢失；另一方面，增龄也会使免疫系统持续低度活化，最终形成促炎症状态[3]。一系列炎症因子，如TNF-α、IL-1、IL-6、IL-7及IL-17等均能诱导M-CSF和RANKL的表达，刺激破骨细胞，最终造成骨量减少，导致OP。并且雌激素和雄激素在体内均具有抗氧化应激的作用，老年男性性激素结合球蛋白持续增加，使睾酮和雌二醇的生物利用度下降，体内的活性氧类堆积，促使间充质干细胞、成骨细胞和骨细胞凋亡，最终导致骨形成减少[39]。总之，老年男性和女性均会表现出炎症表型，尽管目前炎症调控OP的机制十分复杂，但上述结果充分说明了炎症状态的改变很大可能是引发OP的关键因素。

本研究使用的两样本MR的研究方法，从遗传学角度评估了炎症对OP和BMD的因果效应，具有一定的优势。首先，本研究自身能够克服混杂因素和反方向因果的干扰，极大地提高了研究结果的准确性，同时还避免了人力、物力、财力等的浪费，极大提高了工作效率。其次，双能X线吸收测定法检测BMD被誉为OP诊断的“金标准”,因此本研究将各年龄段和全年龄段的BMD均纳入结局进行MR分析，不仅增加了骨质疏松这一结局的可靠性，还将BMD分为各个年龄段进行分析，使得结局更加全面和系统，最终使得本研究的结果更加准确可靠。但该研究依旧存在一些局限性。首先，孟德尔随机化研究通常仅考虑单个或两个变量，但实际情况往往十分复杂，可能会有其他变量会影响实验结果。此外，由于本研究中使用的主要是欧洲人群的数据，这样虽然能够极大程度地减少环境因素、种族特征和饮食结构等差异造成的影响，但这些结果可能并不适用于其他种族群体。最后，虽然孟德尔随机化研究推断了暴露的遗传易感性对结局的潜在因果影响，但它没有提供证据表明治疗干预引起的暴露变化会对结果产生的重大影响。因此，在实际情况中运用孟德尔随机化分析问题时，要综合考虑多种因素以保证实验结果的准确性。

综上所述，本研究使用了孟德尔随机化分析方法、药物预测和分子对接方法来探讨炎症因子与OP之间的因果关系。研究结果表明，炎症确实是OP和BMD的影响因素，提示未来在治疗OP时，需要关注对炎症的调节，进而更好地指导临床。同时这些发现还为设计针对这些靶点的骨质疏松治疗方法提供了可行性，为未来的治疗策略开发提供了重要的理论和实践基础。

参考文献:

[3]中华医学会骨质疏松和骨矿盐疾病分会,章振林.原发性骨质疏松症诊疗指南(2022)[J].中国全科医学,2023,26(14):1671-1691.

[4]何姣姣,陈以发,陈玉林,等.绝经后骨质疏松的骨免疫学机制[J].中国骨质疏松杂志,2023,29(7):1032-1036.

[5]郭曹培,程飘涛,杨成兵,等.骨免疫与骨代谢[J].中国组织工程研究,2024,28(14):2261-2266.

[7]宋昕,孙莹,张琦,等.基于RANKL/RANK/OPG通路研究绝经后老年骨质疏松患者骨代谢与肠道菌群失调的关系[J].河北医科大学学报,2022,43(3):302-306.

[15]刘明,高亚,杨珂璐,等.孟德尔随机化研究的报告规范(STROBE-MR)解读[J].中国循证医学杂志,2022,22(8):978-987.

[27]韩学明,孙忠良,徐建华,等.老年人血清炎症因子与骨密度变化及骨吸收相关性研究[J].中国骨质疏松杂志,2020,26(2):221-226.

基金资助:江苏省科技计划专项资金(基础研究计划自然科学基金)(BK20220468);江苏省研究生科研与实践创新计划项目(KYCX24_2332,KYCX24_2251,KYCX23_2208);

文章来源:王亮,郭翔云,秦金然,等.两样本孟德尔随机化研究炎症因子与骨质疏松症的因果关联[J].中国骨质疏松杂志,2025,31(01):62-69.