生酮饮食导致小鼠骨质疏松的转录组学

吴秀华，范应静，叶永浓，李萍，朱青安，陈泽森，李博，王文，郑磊

南方医科大学南方医院1检验医学科，2脊柱骨科，广东 广州 510515；3广州市番禺区中医院药学部，广东 广州511400

生酮饮食（KD）是一种高脂肪、充足蛋白质和极低碳水化合物的配方饮食方案，它模仿机体禁食状态时的代谢来诱导酮体的产生，可提高血液中β-羟丁酸水平［1］。近年来KD治疗在神经退行性疾病、内分泌和代谢性疾病、心脏疾病、炎症性疾病以及肿瘤中取得了一定成效［2］。然而其不良反应如尿液改变、骨密度降低、神经病变、心血管疾病和肾结石的发生等引起关注［3,4］。早在1979年报道，KD疗法与骨质疏松相关［5］，KD治疗难治性癫痫时，患者骨矿物质含量进行性丢失［6］，这些临床研究提示，KD治疗癫痫患者骨折的风险不仅来之癫痫发作，而且因骨密度（BMD）降低而增加骨折风险［7］；儿童正处于生长阶段，KD用于治疗儿童耐药性癫痫尤其要关注其对骨骼的负面影响。我们针对KD对骨代谢的影响进行研究，在小鼠模型中发现KD降低了松质骨的骨密度、骨体积分数、骨小梁的连接程度，减少皮质骨的横截面积，降低骨力学强度，导致骨质疏松，骨质疏松的程度接近经典去势模型；且KD小鼠的骨髓腔中成骨细胞减少，破骨细胞增多［8,9］。骨髓中细胞包括骨髓间充质干细胞、成骨细胞、破骨细胞和血管内皮细胞等多种类型细胞。骨髓细胞可产生可溶性因子进行细胞间通讯，形成骨髓微环境［10］，任何骨髓微环境的干扰都可能导致骨稳态的破坏［11］。我们前期根据高脂饮食及热量限制饮食的文献报道总结了KD导致骨质疏松的原因，骨髓微环境改变是其中原因之一［12］，但具体机制尚不清楚。

转录组学测序可对特定组织或细胞在特定状态下转录的所有mRNA进行测序，比对参考基因组，进而实现mRNA序列、丰度、基因结构和新转录本等的分析，探寻差异表达基因及信号通路，已在疾病机制探讨、疾病分型、疾病生物标志物研究等领域发挥重要作用［13］，但未见转录组学用于KD导致骨质疏松的研究报道。

本研究通过提取KD小鼠骨髓细胞的RNA，行转录组学测序和生物信息学的分析，筛选与正常饮食组的差异基因及相关通路，并进一步通过RT-qPCR验证差异基因，为探究KD导致骨质疏松的分子机制奠定基础，为预防KD引起骨质疏松提供思路。

1 材料和方法

1.1 动物模型与分组

16只8周龄，体质量16～18 g的雌性C57BL/6J小鼠由南方医科大学动物中心提供（伦理编号为：NFYY-2021-0969），随机分为生酮饮食（KD）组和正常饮食（Sham）组，每组8 只。KD 组小鼠喂养生酮饲料（Zeneca），其碳水化合物与脂肪的比例为1：3；Sham组给予普通饲料喂养，所有动物均喂养3个月。

普通饲料购于南方医科大学实验动物中心。每100 g普通饲料热量为1338 kJ，含蛋白质14.5 g，脂肪4 g，碳水化合物55.5 g，膳食纤维4.5 g，钙720 mg，磷600 mg，维生素D 2.5 μg；每100 g 生酮饲料热量为2804 kJ，含蛋白质18.2 g，脂肪65.1 g，碳水化合物2.7 g，膳食纤维7.4 g，钙500 mg，磷300 mg，维生素D 2.5 μg。

1.2 体质量、血糖和血酮监测

每14 d（2周）每组随机选取6只小鼠测量其体质量，每次每组选取3只小鼠通过剪尾采血，并用血糖仪（JPS-5怡成）、血酮仪（MeterT-1 Sentest Inc.）测量血糖和血酮。

1.3 骨微结构分析

喂养3个月后，用异氟烷气体麻醉过量麻醉使小鼠安乐死。取小鼠双侧股骨，左侧用于提取骨髓细胞中的RNA，行转录组学检测。右侧股骨置于4%多聚甲醛中浸泡固定24 h，用于骨微结构分析。

将获取的小鼠右侧股骨垂直放入Micro-CT 标本杯中，行Micro-CT（μ80，瑞士Scanco）扫描分析。预扫完毕后，选择股骨远端1/3的区域进行扫描。设置扫描参数：电压为55 kV，电流为145 μA，层厚为12 μm。扫描完毕选取图层中没有股骨髁的层面，数100层进行松质骨显微结构分析，包括骨体积分数（BV/TV）、骨小梁数目（Tb.N）、骨小梁厚度（Tb.Th）、骨小梁间隙（Tb.Sp）、各向异性程度（DA）、骨小梁连接密度（Conn.D）、结构模型指数（SMI）、骨密度（BMD of TV）、骨组织密度（BMD of BV）等参数。

1.4 转录组学分析

提取小鼠股骨骨髓细胞总RNA并进行浓度、纯度及完整性检测，质检满足要求后，RNA样本质检结果见表1。遵照试剂盒（TruSeq Stranded mRNA LT Sample Prep Kit）说明书进行操作。纯化提取总RNA中的信使RNA（mRNA）并打断mRNA成为200-300 bp片段。进一步合成第一链cDNA及第二链cDNA，构建链特异性文库。对双链cDNA进行末端修复补平，在3’端加个A，再在链接酶的作用下加上测序接头，然后进行片段的选择，去除多余的接头序列，PCR扩增富集文库片段，文库大小在300～400 bp。检测文库大小（Agilent2100 Bioanalyzer），荧光定量检测文库总浓度（Quantifluor-ST fluorometer，Promega；Quant-iT PicoGreen dsDNA Assay Kit，Invitrogen），qPCR 定量检测有效文库浓度（Thermo Scientific StepOnePlus Real-Time PCR Systems），根据文库的有效浓度、文库需要得到的数据量，将含有不同index序列的文库按比例进行混合。混合文库统一稀释到2 nmol/L，通过碱变性，形成单链文库。Illumina HiSeq平台上机测序，以单链文库为模板进行桥式PCR扩增、测序引物退火、边合成边测序。样本处理、上机测序及生物信息分析由南京派森诺基因科技有限公司完成。

表1 RNA样本质检表Tab.1 Quality check table for RNAsamples

1.5 生物信息分析

采用Excel软件t检验对得到的转录组学原始数据进行P值计算，再按照倍数变化值的计算公式计算，过滤掉接头以及低质量序列并筛选出P＜0.05且log2FC绝对值＞1的基因。在DAVID网站进行基因本体（GO）富集和京都基因与基因组百科全书（KEGG）通路分析。

1.6 实时荧光定量PCR

按RNAiso Plus（Takara）说明书提取总RNA，检测小鼠骨髓细胞RNA浓度，采用逆转录和扩增两步法对DEGs法进行RT-qPCR验证。引物由上海生工生物工程技术服务有限公司合成。将目的基因扩增结果经内参校正后，以空白组的目的基因扩增结果作为对照，2-△△Ct法比较不同样品间mRNA的表达量差异。引物信息见表2。

表2 逆转录PCR 引物序列Tab.2 Primer sequences for reverse transcription PCR

1.7 统计分析

实验结果采用SPSS Statistics软件进行统计学分析，计量资料以均数±标准差形式表示，不同时间点两组动物体质量、血糖、酮体水平变化采用重复测量的方差分析进行评价，两组间骨结构参数采用独立样本t检验分析。设置检验水准α=0.05，以P＜0.05 为差异有统计学意义。

2 结果

2.1 小鼠的体质量、血糖、血酮结果

两组小鼠在2周、4周、6周、8周、10周和12周测量的体质量、血糖水平、血酮水平如图1所示，两组间体质量、血糖水平在各个时间点的差异均无统计学意义。两组间血酮体水平在各个时间点的差异均有统计学意义（P＜0.05）。

图1 生酮饮食和正常饮食小鼠的体质量、血糖、血酮水平Fig. 1 Body weight(A),blood glucose(B)and blood ketone(C)of the mice in the Sham and KD groups.*P＜0.05 vs Sham group.

2.2 小鼠股骨骨微结构

Micro-CT的2D和3D图可见，KD小鼠松质骨的骨小梁间隙变大，骨小梁数量减少，骨密度降低（图2）。

图2 远端股骨骨松质Micro-CT 图Fig. 2 Micro-CT images of the distal femurs of the two groups.

KD 组小鼠股骨远端松质骨骨体积分数（P=0.0039）、骨小梁数量（P=0.0001）、骨小梁连接密度（P=0.0021）、骨密度（P=0.0050）低于Sham组，且骨小梁分离度升高（P=0.0158），差异均有统计学意义。两组的骨小梁厚度（P=0.8895）、骨小梁结构模型指数（P=0.5306）、骨组织密度（P=0.9927）没有明显差异（表3）。

表3 生酮饮食和正常饮食小鼠股骨远端松质骨的结构参数Tab.3 Structural parameters of the cancellous bone in the distal femur in the two group mice

2.3 DEGs分析

采用转录组学技术对Sham组和KD组的基因表达进行了分析比较，测量基因总数16 667 个。以差异倍数log2FC＞1 或log2FC＜-1 为变化阈值，Pval值＜0.05 作为筛选差异基因的标准，结果如图3A所示，共有165个差异基因，其中91个基因上调，74个基因下调，按Pval值大小排序，前20位差异基因如图3B所示。筛选KD对比Sham小鼠差异基因Padj＜0.05且FC＞2的基因作为进一步研究的基因，共筛得5个差异基因，分别为Acot1、Mpig6b、Gp9、Ppbp、Slc2a9（表4）。

图3 Sham组与KD组显著差异基因表达分布及前20位基因聚类图Fig. 3 Distribution of the significantly differentially expressed genes(A)and cluster map of the top 20 differential genes (B) between Sham and KD groups.In A,the gray dashed line represents the threshold line of the differential gene screening criteria,the blue dot represents the down-regulated gene with significant difference,the red dot represents the up-regulated gene with significant difference,and the gray dot represents the non-significant differential gene.

表4 生酮饮食对比正常饮食小鼠显著差异基因Tab.4 Significantly differentially expressed genes in ketogenic diet group compared with normal diet group

2.4 GO富集分析

以log2FC＞1 或log2FC＜-1 为变化阈值，Pval值＜0.05 作为标准筛选差异基因进行GO富集分析，分析结果显示，KD治疗组相比Sham组的差异基因主要富集在黏膜的免疫反应、器官和组织特异性免疫反应、细胞对有机物质的反应、运动、对外界刺激的反应等生物过程（图4）。细胞组分的富集结果显示，差异基因主要分布在细胞外区部分、细胞外间隙、细胞外区、纤维胶原三聚体、带状胶原纤维等部位。

图4 Sham组和KD组差异基因GO功能富集的气泡图Fig. 4 Bubble diagram of GO functional enrichment of the differential genes between Sham and KD groups.The color and size of the bubbles indicate the size of the P value of each item and the number of differential genes that the item contains,respectively.

2.5 KEGG Pathway分析

对log2FC＞1 或log2FC＜-1 为变化阈值，Pval值＜0.05 作为标准筛选的差异表达基因结果进行KEGG通路富集分析，结果如图5所示，差异基因主要参与了人类疾病、机体系统、环境信息处理、细胞过程等。在环境信息处理信号通路中，Apelin信号通路、PI3K-Akt信号通路、ECM受体互作通路、MAPK信号通路、病毒蛋白与细胞因子和细胞因子受体的相互作用通路的富集程度和差异基因数均较高，尤以Apelin信号通路、PI3K-Akt信号通路、ECM受体互作通路最明显。

图5 Sham组与KD组差异基因的KEGG Pathway分析Fig. 5 KEGG pathway analysis of the different genes between Sham and KD groups.

2.6 RT-qPCR验证显著差异基因表达

基于转录组学结果，为了进一步验证和明确KD组相较于Sham组基因表达差异，筛选Padj＜0.05且FC＞2的显著上调差异基因，采用RT-qPCR技术验证。结果显示，各基因均明显上调，Acot1上调倍数为2.49±0.665（P=0.01），Mpig6b上调2.58±0.474倍（P=0.02），Gp9上调1.58±0.199倍（P=0.01），Ppbp上调2.59±0.611倍（P＜0.01），Slc2a9上调1.45±0.162倍（P=0.00），其中上调倍数最明显的基因为Acot1、Mpig6b和Ppbp（图6）。

图6 KD组关键基因验证结果Fig. 6 Verification the differentially expressed genes in KD-fed mice using RT-PCR (Mean±SD).*P＜0.05,**P＜0.01 vs Sham group.

3 讨论

本研究通过转录组学分析KD小鼠骨髓细胞的差异基因及通路富集，发现骨髓细胞的Acot1、Mpig6b、Gp9、Ppbp、Slc2a9基因显著上调，Apelin 信号通路、PI3K-Akt信号通路、ECM受体互作通路、MAPK信号通路、病毒蛋白与细胞因子和细胞因子受体的相互作用通路的富集程度和差异基因数均较高，前三者信号转导通路富集程度显著。进一步RT-qPCR实验表明KD小鼠骨髓细胞上调最明显的基因为Acot1、Mpig6b和Ppbp。

酰基辅酶A硫酯酶1（Acot1）作为一种胞质酶，位于细胞质中，启用酰基辅酶A水解酶活性，参与酰基辅酶A代谢过程、长链脂肪酸代谢过程和非长链的脂肪酸代谢过程。已有研究表明，在禁食期间Acot1 调节PPARα（过氧化物酶体增殖物激活受体）耦合脂肪酸通量与氧化能力，且Acot1敲低能降低PPARa基因的表达，Acot1 过表达会驱动PPARα活性［14］。PPARα是PPAR三种亚型之一，PPAR另一亚型PPAR-γ能介导骨髓间充质干细胞向成脂分化，减少了骨形成［15］，PPARα在骨代谢中的报道比较少。有学者推测PPARα的激活参与骨髓间充质干细胞增殖、分化以及骨代谢，抑制氧化应激缓解骨质疏松［16］。本研究骨髓细胞中Acot1基因显著上调，Acot1上调可能促进PPARα活性，但是否导致骨质疏松，需要进一步研究。

Mpig6b在人和小鼠中具有直系同源生理作用，调节血小板的产生和功能，鼠中Mpig6b位点的缺失可导致大血小板减少和骨髓纤维化［17］。另有文献［18］报道，Mpig6b基因缺乏的雌性小鼠从8周龄开始出现股骨皮质骨厚度增加，骨髓面积减少，到16周龄时髓腔被骨小梁阻塞；Mpig6b缺乏雄性小鼠仅在骨干近端的髓腔发育少量的骨小梁，并且在16周时表现出骨体积分数和小梁厚度的暂时降低。由此可见Mpig6b基因与骨代谢相关，Mpig6b基因缺乏导致骨硬化，那么其过表达可能导致骨质疏松。在本研究中，我们发现Mpig6b基因上调，我们推测Mpig6b基因上调可能是KD导致骨质疏松机制之一，需要进一步研究。

Ppbp（CXCL7）编码的蛋白质是一种血小板衍生的生长因子，属于CXC趋化因子家族，是中性粒细胞的引诱剂和激活剂，能刺激各种细胞过程。Ppbp作为一种趋化因子，被鉴定为一种新的MMP-13底物和破骨细胞形成的调节因子［19］。有文献报道，RANKL和M-CSF存在下用阻断抗体处理RAW264.7细胞和小鼠BMC，阻断SDF-1，CXCL7和CX3CL1，强烈抑制小鼠BMC中的TRAP活性，表明Ppbp趋化因子是从小鼠BMC分泌的，并且与破骨形成有关［20］。在本研究中，我们发现KD小鼠Ppbp（CXCL7）基因上调，这可能促进了破骨细胞的形成及分化，从而导致骨吸收能力增加，破坏了骨平衡。

基于转录组学，我们发现差异基因在Apelin信号通路、PI3K-Akt信号通路、ECM受体互作通路等信号通路富集高。Apelin是一种参与骨稳态的内源性脂肪因子，Apelin-13激活了BMSCs中的线粒体自噬，改善氧化应激，恢复成骨功能［21］。PI3K-Akt通路是经典的骨代谢通路之一，由杜仲组成的混合物清氏丸通过PI3KAkt通路和ATM通路抑制铁死亡，促进细胞存活，对骨质疏松症有治疗作用［22］。也有研究表明镉通过P2X7/PI3K/AKT通路介导了成骨细胞和破骨细胞分化［23］，引起骨质疏松。ECM受体互作通路涉及一系列生物学过程，如细胞分化、增殖和凋亡，ECM受体互作通路与骨质疏松症有关［24］，此外，有研究报道了ECM受体互作通路调节人MSC的成骨分化［25］。由此可见，以上三条通路可能是KD导致骨质疏松的关键通路，需要进一步验证。

此外，虽然差异基因在MAPK信号通路富集不如前三条通路高，已有研究表明，MAPK信号通路可以通过激活多种骨标志物的表达调节成骨细胞的分化［26］，MAPKs及其下游转录因子如ERK1/2、p38和JNK已被发现与其他诸如碱性磷酸酶、骨钙素、骨桥蛋白等多种成骨细胞表型标记物协同激活［27］。MAPK通路也可能是KD治疗导致的骨生成和骨吸收不平衡进而形成骨质疏松的通路。

本研究有几个局限性，首先，我们仅使用雌性小鼠作为研究对象，虽然骨髓细胞中的雌激素水平有限，然而雌激素的干扰可能无法完全排除。如果能同时研究雄性小鼠和雌性小鼠在KD代谢中骨髓细胞基因的改变更加严谨。其次，本研究只用RT-qPCR验证差异基因，后续研究我们将用KD不同代谢产物干预骨相关细胞，并验证其差异基因及蛋白表达，深入研究KD引起骨质疏松的分子机制，为KD引起骨质疏松提供预防思路。

综上所述，本研究表明KD小鼠导致骨质疏松的差异基因在多个细胞外组分与信号通路中富集，其中Acot1、Mpig6b、Ppbp为差异显著基因，可能在KD引起的骨质疏松中发挥重要作用。通路富集中以Apelin信号通路、PI3K-Akt信号通路、ECM受体互作信号通路富集程度显著。本研究为KD引起骨质疏松的分子机制研究提供理论和实验依据，为预防KD导致骨质疏松提供思路。