李生强 陈赛楠 谢冰颖 陈娟 谢丽华 叶云金 黄景文 葛继荣

福建省中医药科学院骨质疏松证候基因组学研究室，福建 福州 350003

骨质疏松症(osteoporosis，OP)可造成骨折及骨折并发症引起的残疾甚至死亡，骨质疏松的防治已成为老年人健康的重大课题。近年来，中药治疗骨质疏松症得到广泛认可。续苓健骨方是以中医药理论遣方用药形成的治疗OP的经验方。前期的动物实验表明，续苓健骨方可提高骨质疏松模型大鼠的骨密度[1]，改善骨生物力学性能[2]，提升血钙含量，降低血磷含量[3]，并对其机制进行了一定探索[4-5]。miRNA参与了骨的生长、发育、代谢及药物治疗等各个过程[6-7]。本研究采用二代测序方法，对骨质疏松大鼠模型在续苓健骨方治疗后进行腰椎miRNA测序，探讨小分子非编码RNA在治疗过程中的变化，为临床应用该方治疗OP提供实验室依据。

1 材料与方法

1.1 材料

1.1.1实验动物：6月龄SPF级SD雌性大鼠18只，购自上海斯莱克实验动物有限责任公司，实验动物许可证号：SCXK(沪)2007-0005，合格证编号：2007000530078。饲养于福建省中医药科学院比较医学实验中心，实验动物许可证：SYXK(闽)2016-0005。

1.1.2实验药物与试剂：续苓健骨方(国家发明专利号ZL201510784227.4)由川续断、白术、陈皮、赤芍、红花、甘草等12味中药组成，共122 g/剂。中药饮片购自福建省医药公司，按传统中药煎熬方法，每剂取汁 81.33 mL，生药含量为1.5 g/mL，4 ℃冰箱保存备用。大鼠灌胃量计算方法：按成人50 kg：122/50=2.44 g/(kg·d)，大鼠灌胃量：2.44×6.25倍=15 g/kg。

文库构建、测序相关试剂：NEB Next®Poly(A) mRNA 磁性分离模块、NEB Multiplex Small RNA Library Prep Set(美国New England Biolabs公司)；RiboZero Magnetic Gold Kit (Human/Mouse/Rat) (美国Epicentre公司)；TruSeq Rapid SR Cluster Kit (#GD-402-4001，美国Illumina公司)。

1.1.3仪器：Discovery W 双能 X 线骨密度仪(变异系数1.0 CV%，精度0.25%，美国Hologic 公司)；NanoDrop ND-1000(美国Thermal Scientific公司)；生物芯片分析系统Agilent 2100 Bioanalyzer(美国Agilent公司)；Illumina NextSeq 500 测序仪(美国Illumina公司)。

1.2 方法

1.2.1动物分组、造模及给药：大鼠按随机数字表分成3组：假手术组、模型组及续苓健骨方组，每组6只。根据参考文献[1]造模。术后 30 d 开始灌胃给药，续苓健骨方 15 g/kg，上午灌胃；假手术组和模型组给予生理盐水 15 g/kg，连续给药12周。各组动物均普通饲料饲养，自由饮水和活动。

1.2.2取材：最后一次灌胃完成2 h后，以10%水合氯醛麻醉后，取左侧胫骨和第一、二腰椎。胫骨用于骨密度检测；每组随机选3只大鼠腰椎用于提取RNA及测序。

1.2.3骨密度检测：采用双能X线骨密度仪测定各组大鼠左侧胫骨骨密度，操作过程严格按仪器使用说明进行。

1.2.4样品RNA提取及质检：从大鼠第一腰椎提取总 RNA，Nanodrop ND-1000测定RNA 260/280的浓度及蛋白质污染情况，采用变性琼脂糖凝胶电泳检测RNA纯度及完整性。

1.2.5文库构建及测序：构建好文库，用Agilent 2100 Bioanalyzer进行文库质量测定，混合好的不同样品的测序文库，通过0.1 M NaOH变性生成单链DNA，然后在Illumina NextSeq 500测序仪上进行51循环测序。文库构建及测序均委托上海康成生物工程有限公司完成。

1.2.6测序数据分析：采用 miRDeep2 0.0.8 reads 软件把trimmed miRNA 进行已知 miRNA定量和新miRNA预测。使用R软件edgeR 3.18.1进行差异表达计算，筛选差异miRNAs，在线VENNY2.1分析共同差异miRNAs。利用miRNA靶基因数据库统计top10差异miRNA的靶基因，并进行靶基因GO功能显著性富集分析、靶基因Pathway显著性富集分析。

1.3 统计学方法

使用edgeR筛选差异表达miRNA时，设定阈值为1.5倍差异，P<0.05且组内CPM均值≥1。采用 SPSS 17.0 软件进行数据分析，计量资料用均数±标准差表示，根据数据是否正态分布，选择非参数检验或t检验比较组间差异，P<0.05 表示差异有统计学意义。

2 结果

2.1 各组大鼠胫骨骨密度比较

与假手术组相比，模型组大鼠骨密度显著下降(P<0.01)；续苓健骨方组骨密度较模型组显著提高(P<0.01)；续苓健骨方组与假手术组相比无显著性差异(P>0.05)，见表1。

表1 各组大鼠胫骨骨密度值比较

2.2 miRNA表达谱组间差异分析

各组之间差异miRNAs数见表2。VENNY 2.1在线进一步分析模型组与假手术组的差异表达基因在续苓健骨方组中表达情况。模型组表达上调，在续苓健骨方组表达下调的miRNAs有50条；模型组表达下调，在续苓健骨方组表达上调的miRNAs有60条，即模型组共有110条异常表达的miRNAs在续苓健骨方组得到纠正。表3列出了这110条差异miRNAs中，续苓健骨方治疗之后上调或下调的top 10 miRNAs，它们在假手术组、模型组、续苓健骨方组的表达变化及其已知功能。

表2 组间差异miRNA数比较

表3 前10位差异表达 miRNAsTable 3 Top 10 of differential expressed miRNAs

2.3 实时荧光定量PCR验证miRNAs表达

根据miRNA测序的结果，从表3中随机选取了4个miRNA进行定量PCR验证，内参选择 snRNA U6。结果表明，rno-miR-136-3p、rno-miR-26b-5p在模型组中显著提高(P<0.01)，在续苓健骨方组中显著下降(P<0.01)；rno-miR-485-5p、rno-miR-138-5p在模型组中显著降低(P<0.01)，在续苓健骨方组中显著上升(P<0.01)。定量PCR的结果与测序结果一致(图1)。

图1 定量PCR验证部分测序结果Fig.1 Partial sequencing results verified by real-time PCR注：与假手术组比较，*P<0.01；与模型组比较，†P<0.01

2.4 差异miRNAs靶基因预测及GO富集分析

基于miRNA数据库，从110条miRNA中选择经过续苓健骨方治疗后上调或及下调的top 10 miRNAs，进行靶基因预测，然后对靶基因进行Gene Ontology(GO)分析，包括分子功能(molecular function，MF)、细胞组成(cellular component，CC)及参与的生物学过程(biological process，BP)。图2分别列出了差异最显著的前10个GO条目。这些靶基因参与了蛋白结合、多肽结合、转录因子调节、RNA聚合酶转录因子活性等过程。

图2 GO分析柱形图 A：表达上调miRNAs靶基因参与的GO分析；B：表达下调miRNAs靶基因参与的GO分析Fig.2 Column map of enrichment analysis. A: Go analysis of target genes involved in up regulated miRNAs; B: Go analysis of target genes involved in down regulated miRNAs.

2.5 差异miRNAs靶基因信号通路分析

根据预测的靶基因，用其在KEGG数据库进行信号通路分析。分别选取前10个信号通路，根据P值绘制成柱状图(图3)。表达上调的miRNA的靶基因参与的信号通路包括甘油磷脂代谢通路、磷脂酰肌醇信号通路、HIF-1信号通路、甲状腺激素信号通路、Hedgehog信号系统、细胞内吞信号通路等；表达下调的miRNA的靶基因参与的信号通路主要包括内质网蛋白质加工信号通路、谷氨酸突触信号通路、GnRH信号通路、自噬调节、CAMP及MAPK等信号通路等。这些信号通路很可能参与了续苓健骨方治疗骨质疏松模型大鼠的治疗过程。

2.6 预测新的miRNAs

通过使用Dicer对miRNA前体处理的简单模型，miRDeep2预测新的miRNAs。表4列出了miRDeep2得分最高的10个新miRNAs。

图3 KEGG信号通路分析 A：表达上调miRNAs靶基因参与的信号通路；B：表达下调miRNAs靶基因参与的信号通路Fig.3 Pathway analysis on KEGG. A: signaling pathways involved in up regulated miRNAs; B: signaling pathways involved in down regulated miRNAs

表4 新miRNA表达数据展示Table 4 Display of expression data of new miRNAs

3 讨论

中医药治疗骨质疏松的miRNA表达谱的研究较少。赵清等[20]应用Illumina高通量测序技术检测中药骨康治疗绝经后骨质疏松症(postmenopausal osteoporosis，PMOP)SD大鼠血浆miRNAs，认为差异miRNAs可能在PMOP发生发展以及中药治疗中起重要的调控作用。严芳娜[21]用miRNA芯片检测黄精多糖对大鼠骨髓间充质干细胞诱导分化为成骨细胞和骨髓源性单核巨噬细胞诱导分化为破骨细胞过程中miRNAs表达谱的变化，认为miRNA的表达对诱导成骨细胞和破骨细胞均有重要作用。陈哲[22]运用高通量测序技术检测并qPCR验证滋肾降糖丸可能通过miRNAs对Wnt通路进行调控，进而影响成骨分化、实现防治糖尿病骨质疏松的作用。

本动物实验中，续苓健骨方组大鼠骨密度则显著提升。从差异miRNAs来看，续苓健骨方对在模型组出现异常表达的miRNAs起了较好的恢复作用。经12周的治疗后，共有110条miRNAs的异常表达得到纠正。与假手术组相比，续苓健骨方组仅有8个异常表达的miRNAs，这表明无论是骨密度还是miRNA层面，续苓健骨方对骨质疏松症均具有良好的治疗作用。

在筛选到的差异miRNAs中，许多miRNAs均与骨质疏松症或骨代谢密切相关，如Chen等[11]发现miRNA-136-3p在酒精诱导的小鼠骨量减少模型中的表达显著降低，而且它能靶向PTEN调节血管生成和骨形成并改善酒精诱导的骨量减少。Panach等[23]从临床患者血液中筛选到miR-122-5p，定量PCR证明其与骨质疏松性骨折相关；Mandourah等[13]采用芯片从临床患者外周血筛选到差异表达miR-122-5p，定量PCR验证它与低骨密度及脆性骨折相关，认为是一个潜在的诊断标记物；Liao等[24]也发现miR-122-5p负性调控SPRY2，升高受体酪氨酸激酶(RTK)活性，从而促进成骨细胞的增殖和分化。Zhang等[14]发现在骨质疏松症患者中miR-485-5p表达水平上调，而在成骨分化过程中，miR-485-5p受到抑制，后又通过数据库预测及荧光素酶报告基因证实miR-485-5p靶向Ostex促进骨质疏松。此外，本文在测序的基础上，通过软件预测得到了高信度的新miRNAs，不过其具体功能还有待进一步深入研究。

miRNA是机体重要的调控体系，对相关信号通路的调控是miRNA发挥作用的重要方式。本研究对top10差异miRNA的靶基因进行信号通路分析，发现靶基因参与HIF-1、甲状腺激素、Hedgehog、GnRH、自噬、cAMP及MAPK等信号通路。这些信号通路均参与了骨代谢的过程。

综上所述，本研究采用Illumina测序结合生物信息学分析，建立了卵巢切除大鼠模型骨组织miRNA表达谱，分析了续苓健骨方对骨质疏松模型大鼠腰椎miRNAs作用，筛选出差异表达miRNAs。top10 miRNAs靶基因参与了蛋白结合、多肽结合、转录调控、RNA聚合酶转录因子活性等生物学过程，及HIF-1、甲状腺激素、Hedgehog、GnRH、自噬、cAMP及MAPK等信号通路，对于从miRNA角度探讨骨质疏松症的形成及中药治疗骨质疏松症的分子机制具有重要参考价值。