基于转录组测序分析壮骨止痛方治疗去势大鼠骨质疏松症的基因表达差异

陈红琼，杨珊珊，陈瑶，郁洁，雷晓明

（湖南中医药大学，湖南 长沙 410208）

骨质疏松症是一种以骨量减少，骨组织微结构损坏，导致骨脆性增加，易发生骨折为特征的全身代谢性骨骼疾病［1］。其主要临床表现为疼痛、腰膝酸软、脊柱变形、易于骨折［2］。研究表明，骨骼代谢受到许多因素调节，包括激素、年龄、遗传以及环境等［3］。其主要发病机制是成骨细胞与破骨细胞之间骨重建过程的失衡，骨形成与骨吸收两者之间的失衡［4］，但目前骨质疏松症发病的分子机制尚不明确。骨质疏松症没有根治的方法，现代临床医学治疗原发性骨质疏松症普遍应用抑制骨吸收类药物、促骨形成类药物和骨矿化类药物，这些药物在临床上取得了显著的治疗效果，但长期应用会带来较多不良反应［5］。

壮骨止痛方是湖南省名中医莫新民教授的经验方，以古方补骨脂丸为基础化裁而来，由淫羊藿、补骨脂、女贞子、枸杞、骨碎补、牛膝和狗脊组成，具有补益肝肾、壮骨止痛功效。前期研究已证明壮骨止痛方可以多途径抑制骨吸收，促进骨代谢，提高雌激素水平，对于绝经后骨质疏松症具有良好的防治作用［6-7］。但其治疗骨质疏松症的具体机制仍不清楚，因此，本研究基于转录组学测序分析壮骨止痛方治疗去势大鼠骨质疏松症的差异表达基因以及对相关通路的影响，寻找相关的基因靶点，以期为临床应用壮骨止痛方提供实验依据。

1 材料

1.1 主要仪器

KL05A 型医用离心机（湖南凯达科学仪器有限公司）；991 型超低温冰箱（赛默飞世尔科技有限公司）；Illumina 测序平台（广州市唯誉智合科技有限公司）。

1.2 药品与试剂

壮骨止痛胶囊（四川美大康药业股份有限公司，批准文号：Z20050118）；戊巴比妥钠（YaErbio公司，批号：20200422）；青霉素钠（华北制药股份有限公司，批号：F9102109）；红细胞裂解液（合肥志宏生物技术有限公司，批号：700577116）；大鼠淋巴细胞分离液试剂盒（北京索莱宝科技有限公司，批号：20200616）。

1.3 动物

SPF 级雌性SD 大鼠12 只，体质量180～220 g，购自湖南斯莱克景达实验动物有限公司，动物合格证号（1107272011003527），于湖南中医药大学动物实验中心饲养，许可证号：SCXK（湘）2019-0009，饲养室温度24～26 ℃，相对湿度50%～70%。人工12 h 昼/夜循环照明，分笼饲养。每日定时清洗笼舍，大鼠能自由摄食及饮水。实验前适应性喂养1 周，造模前禁食12 h，自由饮水。

2 方法

2.1 分组、造模与给药

将12 只SD 大鼠随机分为实验组、模型组，每组6 只。两组均采用国内外公认的雌性大鼠去卵巢3 个月造成绝经后骨质疏松模型［8-9］。方法：大鼠用2%戊巴比妥钠（0.2 mL/100 g）腹腔麻醉，在无菌条件下，从大鼠腰椎两侧摘除双侧卵巢。术后连续3 d 大腿肌内注射青霉素纳，每只大鼠4 万U/天，术后5 d 拆线。术后1 周开始灌胃给药，实验组给予壮骨止痛胶囊（567 mg/kg），模型组给予同体积的蒸馏水（1 mL/150 g），每日1次，连续13周。

2.2 样本制备

13 周后，所有大鼠腹腔注射2% 戊巴比妥钠（0.2 mL/100 g）麻醉，经腹腔主静脉取血后，2 h内用大鼠淋巴细胞分离液试剂盒按说明收集PBMC，红细胞裂解液裂解PBMC 中残余红细胞，随机将各组每2 只大鼠的PBMC 合并为1 个样本，每组3 个样本，于液氮速冻后放-80 ℃保存，备用。

2.3 RNA提取及测序

将两组的PBMC 样本，共6 个样本，送至广州市唯誉智合科技有限公司进行RNA 提取、建库及转录组测序，测序平台为Illumina测序平台。

3 信息分析

3.1 差异表达基因分析

采用DEGSeq 分析差异表达基因，差异基因筛选主要参考差异倍数（Fold change 值）以及q值（Padj 值，矫正之后的Pvalue 值）作为相关指标，通常选取|log2 Fold change|≥1 和q＜0.05 的差异基因作为显著差异基因。

3.2 GO和KEGG富集分析

基于筛选出的差异表达基因（|log2 Fold change|≥1，q＜0.05），应用超几何检验对差异表达基因进行GO和KEGG 富集分析，以校正后的q＜0.05 为阈值，满足此条件的GO词条和KEGG通路被认为具有显著性。

4 结果

4.1 差异表达基因分析结果

实验共获得20 276 个转录组mRNA，与模型组相比，实验组共有10 166 个表达上调，10 110 个表达下调，见图1。为差异表达基因及样本垂直聚类结果的热图。DEGs 分析显示，差异明显的基因有149 个（|log2 Fold change|≥1，q＜0.05），67个基因上调，82个基因下调（见图2）；进一步筛选显著差异基因（|log2 Fold change| ≥2，q＜0.01），共选出18 个显著差异基因，其中上调基因15 个，下调基因3 个，具体信息见表1。

表1 实验组与模型组外周血单个核细胞的差异表达基因

图1 实验组与模型组差异基因的热图

图2 实验组与模型组差异表达基因的火山图

4.2 差异表达基因的GO富集分析

对差异基因（|log2 Fold change|≥1，q＜0.05）进行GO 富集分析，结果显示有45 个富集项（见图3），差异表达基因主要富集在细胞部分、膜部分、细胞器、细胞过程、生物调节、代谢过程、绑定分子类、催化活性等；进一步校正检验后的P值之后，以q＜0.05 为阈值，发现满足此条件的差异表达基因中显著富集的GO 条目有316 条，功能富集到生物过程的GO 条目有303 条（见表2，取显著性最高的前10条）、细胞组分有4条（见表3）、分子功能有9条（见表4）。

表2 差异表达基因显著性GO富集结果-生物过程

表3 差异表达基因显著性GO富集结果-细胞组分

表4 差异表达基因显著性GO富集结果-分子功能

图3 实验组与模型组差异表达基因的GO统计柱状图

4.3 差异表达基因的KEGG富集分析

对差异基因（|log2 Fold change|≥1，q＜0.05）进行KEGG富集分析，共富集到170条生物通路，结果显示：上调基因主要参与细胞因子与细胞因子受体的相互作用、MAPK 信号通路、TNP 信号通路、吞噬体、趋化因子信号通路、癌症通路、癌症转录调控失调等通路（见图4）；下调基因主要参与嘌呤代谢、内质网蛋白加工、cGMP-PKG信号通路、黏着斑、IL-17信号通路、轴突引导等通路（见图5）；进一步校正检验P值之后，以q＜0.05 为阈值，得到TNP 信号通路、IL-17 信号通路、沙门氏菌感染三条显著性富集的通路见表5。

表5 差异表达基因的显著性富集通路

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

图5 下调差异表达基因的KEGG分类

5 讨论

本研究采用高通量测序技术对去势大鼠骨质疏松模型组与使用壮骨止痛方干预的实验组的外周血单个核细胞进行转录组学测序，差异表达基因分析结果表明，实验组与模型组差异表达基因有149 个，以|log2 Fold change| ≥1 和q＜0.05 作为筛选条件，上调基因有67 个，下调基因有82 个。GO 富集分析涉及生物过程、细胞组分以及分子功能三大类，差异表达基因主要富集在细胞部分、膜部分、细胞器、细胞过程、生物调节、代谢过程、绑定分子类和催化活性等；校正P值后，在生物过程分类中，细胞迁移、免疫系统过程、防御反应所占比例最多，其次是应激反应、调节免疫系统的过程、细胞运动性、调节对外界刺激的反应、炎症反应等；在细胞组分分类中，细胞膜部分所占比列最多；在分子功能分类中，细胞因子受体活性所占比列最多，其次是整联蛋白结合及受体结合类分子功能。KEGG 富集分析结果显示，显著性富集的通路为TNF 信号通路、IL-17 信号通路、沙门氏菌感染三条通路，涉及相关差异基因有Cebpb、Cxcl2、Icam1、Jun、Mapk10、Ppbp、LOC100359515。研究发现，肿瘤坏死因子成员之一（TNF-α）可以通过NF-κB 信号转导通路促进成骨细胞的分化和细胞外基质的矿化［10］，白介素17（IL-17）可以显著抑制骨保护素（OPG）和增强NF-κB 配体受体活化剂（RANKL），并且IL-17 可以通过TNF 受体相关因子6（TRAF6）、TANK 结合激酶1（TBK1）、c-Jun N 端激酶（JNK）或NF-κB 调节RANKL 和OPG 的表达［11］，而RANKL/RANK/OPG 系统是调节破骨细胞分化、成熟和凋亡最重要的信号通路［12］。但是暂无沙门氏菌感染在骨质疏松领域的相关报道。

在差异表达基因分析中，差异表达基因主要与脂质代谢相关疾病、骨质疏松疾病、炎症、肿瘤等疾病有关。如显著性最高的G0S2，其参与机体内脂质代谢过程，与肥胖、糖尿病胰岛素抵抗等脂质代谢相关疾病密切相关［13］；Dgat2在脂质代谢旺盛的组织中的表达水平最高，如肝脏、肾脏、白色脂肪组织等［14］；Fpr3参与多种生理和病理过程，包括防御反应、炎症、神经退行性疾病、肿瘤、糖脂代谢紊乱相关疾病等［15］；Nefh 与肌肉萎缩，阿尔茨海默病、帕金森综合征有关；Niacr1 可使小鼠结肠炎和结肠癌风险增加。但目前以上几种基因在骨质疏松领域还没有相关报道，其表达是否会影响骨质疏松症的发生还需进一步验证。

已有骨质疏松相关报道的差异基因如Ccrl2 参与了chemerin 信号通路调节，研究发现，chemerin 与其受体CMKLR1、Ccrl2 在骨髓充质干细胞和前成骨细胞的成脂化或成骨化的过程中起就决定性作用［16］，但目前有关Ccrl2 在骨质疏松领域的报道很少；Csf3r 是一种蛋白质编码基因，其在破骨细胞上表达，在成骨细胞上不表达，与骨质疏松症的形成有一定关系［17］，但在慢性中性粒细胞白血病领域的研究更常见；LOC100359515 是诱导型一氧化氮合酶（iNOS）编码基因，在骨质疏松症的发病机制中起着重要作用。iNOS 产生一氧化氮（NO），可导致骨骼形成和雌激素缺失引起的骨吸收之间的不平衡［18］，这可能是治疗骨质疏松症的重要靶点；Nlrp12 通过抑制NF-κB 诱导的破骨形成在骨中发挥保护作用，并在破骨刺激下下调，调节Nlrp12 水平可控制OC 前体中NF-κB 信号，改变骨稳态和骨溶解反应［19］；破骨细胞相关受体基因Oscar 在破骨细胞分化中起关键作用，可以使骨密度降低，是调控骨密度的重要候选基因［20］。在本研究中Nlrp12、Oscar 是上调基因，其可能是壮骨止痛方治疗骨质疏松症的关键靶点。

目前，关于骨质疏松症的具体发病机制尚不明确，有关骨质疏松症的分子机制研究也比不多，并且没有根治的方法。近年来随着科学技术的迅速发展，高通量测序技术已经广泛应用于基因表达转录水平研究，是研究疾病基因功能的重要手段之一。本研究通过基因差异表达分析和富集分析，筛选了壮骨止痛方干预去势大鼠骨质疏松症相关的基因和通路，但由于样本量不多，随机误差和假阳性难以排除，需要进一步研究，且本研究是通过高通量测序结合基因数据库进行生物信息学的预测，还需后续的实验去验证。