绝经后女性骨质疏松相关基因的生物信息学分析

邓禹，黄禾，金鸿雁

(北京大学第一医院妇产科，北京 100034)

随着人口老龄化的到来，我国有1.2亿女性正在经历向更年期的过渡，这约占世界更年期女性人口的23%[1]。更年期女性在绝经时，由于雌激素缺乏引起许多健康问题，如骨密度(BMD)降低、关节疼痛、潮热、睡眠障碍等症状，严重影响生活质量[2]。据报道，绝经后骨质疏松(PMOP)的发生率为30%～55%[3]，预防和治疗骨质疏松(OP)对避免PMOP的发生至关重要[4]。尽管已有研究分析了影响OP进展的关键分子，但是我们对影响PMOP发生和进展相关基因的了解仍不详尽，需要深入研究[5]。

骨的发育和维持是一个动态平衡的过程，这一过程受到骨破骨细胞和成骨细胞的共同影响[6]。外周血单核细胞(PBMs)是一个研究OP的合适细胞模型[7]。对于成人外周骨骼(如股骨)来说，PBMs是破骨细胞前体的唯一来源[8]。PBMs被招募到骨重塑部位后，在接受特定刺激后形成破骨细胞[9]。研究表明，PBMs可以通过分泌一些细胞因子影响破骨细胞分化、活化和凋亡[10]。PBMs可以通过转化生长因子-β(TGF-β)促进破骨细胞的形成[11]。PBMs的异常与OP密切相关[12]。因此，本研究利用NCBI中的GEO基因芯片公共数据库进行数据分析，选择芯片数据(GSE56815)作为分析对象，筛选差异表达基因，并对上调差异表达基因进行GO和KEGG富集分析，利用String数据库分析关键上调差异表达基因的蛋白质-蛋白质相互作用(protein-protein interaction，PPI)网络，深入了解关键差异表达基因的生物学功能，以期筛选出影响PMOP发生和进展的关键分子。

材料与方法

一、研究对象及分组

利用NCBI(https：//www.ncbi.nlm.nih.gov/geo)中的高通量基因表达数据库(gene expression omnibus，GEO)公共数据库进行芯片数据搜索[13]。检索词为“Postmenopause”或“Bone mineral density”。选择40例样品芯片数据(GSE56815)为研究对象，其中包含了绝经后20例髋部高BMD和20例低BMD女性PBMs的RNA微阵列数据。基因表达谱中髋部高BMD组患者平均年龄为(57.2±1.9)岁，髋部BMD的平均Z值为(1.58±0.66)；髋部低BMD组患者平均年龄为(57.6±1.5)岁，髋部BMD的平均Z值为(-1.04±0.38)。为方便描述，我们将绝经后髋部高BMD组患者的PBMs数据命名为高BMD组，将绝经后髋部低BMD组患者的PBMs数据命名为低BMD组。

二、方法

1.差异表达基因的筛选：本研究使用R软件(3.6.1版)中的Limma程序包对GSE56815数据库中高BMD组和低BMD组的差异基因进行筛选，筛选标准为：|log2 FC|>0.379，错误发现率(false discovery rate，FDR)<0.05。

2.GO和KEGG富集分析：本研究使用R软件(3.6.1版)中的ClusterProfiler程序包对筛选出的上调差异基因进行GO和KEGG富集分析，得到P<0.05的功能富集和KEGG通路。为进一步探究关键差异基因的功能和参与的信号通路，我们使用Metascape数据库(http：//metascape.org/gp/index.html#/main/step1)对PPI网络分析中的关键差异表达基因进行GO功能富集分析[14]。

3.蛋白互作PPI网络分析：采用String数据库(http：//string-db.org)构建差异表达基因的PPI网络的分析图。首先，将所有上调的差异基因输入String数据库，在“Minimum required interaction score”的选项中选择“High confidence(0.700)”生成蛋白互作PPI网络图，然后使用Cytoscape软件(3.8.2版)进行可视化分析。

结 果

一、差异表达基因的筛选

通过对GSE56815数据集进行分析发现，数据库中40例样本的大部分基因表达量基本保持一致，说明该数据适合进行下一步分析(图1)。

图1 GSE56815基因表达箱线图

对GSE56815数据集进一步分析发现，共得到40个表达差异基因，其中上调基因34个(表1)和下调基因6个；并绘制火山图，其中红点示上调基因，蓝点示下调基因(图2)。

图2 火山图

表1 表达上调的差异基因

二、GO和KEGG富集分析

为了解显著表达上调的差异基因的功能，本研究对34个表达上调差异表达基因进行了GO和KEGG富集分析。

上调差异表达基因GO分析发现，主要富集在中性粒细胞脱颗粒、免疫反应的中性粒细胞激活、对细菌来源的分子反应、对脂多糖的反应和细胞对生物刺激的反应等功能(图3)。

图3 表达上调差异基因的GO富集分析

上调差异表达基因KEGG分析发现，主要富集在NOD样受体信号通路、金黄色葡萄球菌感染、成骨细胞的分化、糖类降解和精氨酸的生物合成等信号通路(图4)。

图4 表达上调差异基因的KEGG分析

三、上调差异表达基因的PPI网络及重要基因

采用String数据库构建了上调差异表达基因的PPI网络，其中CEACAM8、LTF、CRISP3、CAMP、LCN2、DEFA4、MMP8、ARG1、HP、BPI、CEACAM6、TCN1和ELANE存在着互作关系(图5)。

图5 表达上调差异表达基因PPI网络

四、关键差异表达基因GO富集分析

采用Metascape分析工具对PPI分析结果中的关键差异表达基因(CEACAM8、LTF、CRISP3、CAMP、LCN2、DEFA4、MMP8、ARG1、HP、BPI、CEACAM6、TCN1和ELANE)进行GO富集分析发现，关键差异表达基因主要富集在中性粒细胞脱颗粒、抗菌肽、细胞外基质组织以及对脂多糖的反应(表2)。

表2 关键差异表达基因的GO富集分析

讨 论

OP是老年人常见的全身性骨病，其特征是骨质减少和骨组织微结构破坏，骨脆性和骨折的风险增加[15]。OP可分为原发性骨质疏松症(POP)和继发性骨质疏松症(SOP)两大类。POP包括PMOP、老年性OP和特发性OP[16]。在绝经期妇女中，雌激素的缺乏会增加OP的发生风险[17]。据报道PMOP是最常见的POP，影响了超过一半的绝经后女性[18]。目前，我们尚未完全了解影响PMOP发生和进展的关键分子[5]。因此，进一步探究影响PMOP发生和进展的机制具有重要的社会意义。

PBMs作为破骨细胞的前体[8]，PBMs的异常与OP密切相关[12]。因此，本研究通过分析GSE56815基因表达谱中40例绝经后不同BMD患者的PBMs微阵列数据，筛选出了高BMD组表达上调差异基因34个和下调差异基因6个。PPI网络分析显示，CEACAM8、LTF、CRISP3、CAMP、LCN2、DEFA4、MMP8、ARG1、HP、BPI、CEACAM6、TCN1和ELANE存在着互作关系。进一步对上述关键差异基因进行GO分析发现，这些关键差异基因主要富集在中性粒细胞脱颗粒、抗菌肽、对脂多糖的反应、细胞外基质组织，其中LTF、CEACAM8、CAMP、LCN2、MMP8和ARG1与骨代谢的联系密切。

乳铁蛋白(LF)，曾称LTF，是一种铁结合的球状糖蛋白，属于转铁蛋白超家族，主要存在于中性粒细胞分泌小泡、上皮分泌物和乳汁中[19]。LF具有重要的免疫调节功能[20]，能减少肿瘤坏死因子-α(TNF-α)和白细胞介素-1β(IL-1β)等多种溶骨细胞因子的分泌[21]。LF存在于中性粒细胞的分泌颗粒中，它的产生受到炎症刺激的影响[22]。LF具有多种生物学效应，除了参与抗菌和免疫调节外，它可以诱导成骨细胞增殖和分化，同时抑制成骨细胞凋亡和破骨细胞的形成，在体内促进局部骨的形成[22-23]。有研究报道，LF能促进成骨细胞系(Saos-2)的成骨分化，以及增强骨再生[24-25]。Hou等[26]研究发现，切除卵巢可显著降低大鼠股骨以及腰2～4椎体的骨量、骨小梁厚度和骨小梁数，并增加骨小梁分离，LF或雌激素治疗可以保护去卵巢大鼠骨量、维持骨小梁厚度和防止骨小梁的丢失。骨保护素(OPG)和核因子激活因子-κB配体(RANKL)是调节骨形成和骨吸收所必需的细胞因子，RANKL/OPG比值是调节骨形成和骨吸收平衡的关键因素。如果缺乏雌激素，将导致这一比值失衡。Hou等[26]进一步研究发现，在切除了卵巢的大鼠的骨小梁中，LF可以下调RANKL/OPG比值，从而防止骨丢失。凌冰等[27]发现，大鼠去除卵巢后，血清LF含量会降低，而且血清LF含量与大鼠的左股骨和腰椎的BMD呈正相关。综上所述，LF在骨代谢中具有重要的作用，血清LF含量与BMD的变化有相关性，LF能否作为预测PMOP变化的指标需要进一步研究。

CEACAMs是免疫球蛋白超家族的成员，参与调节各种细胞功能，如增殖、凋亡和分化等[28]。CEACAM8是CEACAMs的成员之一，又称CD66b，编码糖基磷脂酰肌醇(GPI)连接的糖蛋白，这种糖蛋白仅在人类粒细胞中表达。当中性粒细胞被激活时，中性粒细胞中的CD66b表达上调[29]。有研究发现，当中性粒细胞受到强烈刺激活化后可能会引起颗粒细胞中部分特定基因表达增加，如CD66b和LF，且CD66b的表达水平与LF水平呈正相关[30]。但是，CD66b的功能很大程度上还是未知的，它与LF的相关性以及BMD的关系还需要进一步研究。

CAMP基因位于人3p21.31号染色体上，编码18-kDa的前体蛋白hCAP-18。CAMP基因的表达是由1α，25-二羟维生素D3[1，25(OH)2D3]通过Toll样受体通路调控[31]。有意思的是，1，25(OH)2D3会显著增加破骨细胞中CAMP的表达，但是CAMP的过表达是否会影响破骨细胞的功能需要深入研究[31]。

LCN2也称为中性粒细胞明胶酶相关脂钙蛋白(NGAL)，是脂蛋白家族的重要成员，是脂肪来源的细胞因子。在成骨样细胞MC3T3-E1细胞中白细胞介素-17和肿瘤坏死因子-α会增加LCN2表达[32]。有研究发现，LCN2水平升高可能与老年人OP松性骨折的风险增加有关[33]。

MMP8属于基质金属蛋白酶(MMPs)家族的一员，具有降解细胞外基质的能力，参与骨骼和软骨的发育[34]。Li等[34]发现，在去除了卵巢的大鼠中，MMP13在mRNA水平和蛋白表达水平均显著增加，而MMP8在mRNA和蛋白表达水平均无显著变化；MMP13蛋白表达与骨小梁分离度呈正相关，与骨小梁数目呈负相关。有研究发现，17β-雌二醇可以显著增加成骨细胞中雌激素受体和Sirtuin-1(Sirt1)的表达水平，降低NF-κB和MMP8的表达水平，即17β-雌二醇通过调节 Sirt1/NF-κB/MMP8通路改善成骨细胞功能[35]。

精氨酸酶(ARG)是参与鸟氨酸循环的重要酶之一，它包括定位于细胞质的ARG1和位于线粒体的ARG2两种亚型[36]。ARG1过表达时，抑制破骨细胞分化、骨吸收和一氧化氮(NO)生成[37]。有研究表明，NO可能会刺激破骨细胞分化和活性，进一步导致大鼠骨丢失[36]。ARG1可能通过调节NO来负向调节破骨细胞的分化[37]。骨代谢与CRISP3、DEFA4、HP、BPI、TCN1和ELANE之间的关系尚未见报道。

此外，本文也存在一些局限性，如本研究为了精确筛选PMOP与OP相关的差异基因，减少因数据集的不一致性造成的异质性，我们对纳入的GEO数据集进行了严格筛选，因此最终纳入的患者病例数较少，未来需要扩大样本量进行更深入的研究。

综上所述，本研究通过生物信息学分析，筛选出了与绝经后高BMD组患者表达上调差异表达基因34个，其中CEACAM8、LTF、CRISP3、CAMP、LCN2、DEFA4、MMP8、ARG1、HP、BPI、TCN1和ELANE存在相互作用，且CEACAM8、LTF、CAMP、LCN2、MMP8和ARG1与骨代谢的联系密切，但CRISP3、DEFA4、HP、BPI、TCN1和ELANE与骨代谢的联系还需要进一步研究。本研究筛选出的差异表达基因为研究绝经后女性骨量的变化提供了一系列的生物标记物，深入研究这些基因的生物学功能，将有助于我们了解PMOP发生和进展的分子机制。