Illumina高通量测序技术研究肌少-骨质疏松症与骨质疏松症患者骨组织microRNAs差异表达谱及差异分析

陈锦成 朱国涛 秦晓飞 陈彦丞 罗骏 刘洪文 余博飞 徐杰*

1.福建中医药大学中医骨伤及运动康复教育部重点实验室，福建 福州 350122 2.福建医科大学省立临床医学院(福建省立医院)骨科，福建 福州 350001

肌少症是指与自然增龄相关的，其特征是肌肉强度的下降和(或)全身肌量的减少或人体肌肉的生理活动功能衰退状况呈进行性发展[1-2]。骨质疏松症(osteoporosis，OP)属于全身性骨病范畴，其特征是骨量逐步减少、骨的微观结构退化及脆性增加，其系统性退化伴随骨折发生率的增加[2-3]。肌少-骨质疏松症(sarco-osteoporosis，SO)的概念是由Binkley 等[4]提出的，目前临床上尚未有行业专家一致公认的诊治标准，本课题组前期的临床实验研究规定的SO患者纳入标准是合并骨质疏松(腰椎BMD T评分<-2.5)及肌肉功能异常(肌力0～3 级或出现明显肌肉萎缩)的全髋关节置换患者；OP患者纳入标准是合并骨质疏松但不伴肌肉功能异常的全髋关节置换患者。排除标准：已知体内存在影响骨代谢的疾病，包括但不限于：糖尿病、甲状腺疾病、肿瘤性疾病患者、免疫系统疾病等；近12个月内服用影响骨代谢的药物，包括但不限于：双膦酸盐制剂、利尿剂、肾上腺和甲状旁腺药物、糖皮质激素、维生素D和钙剂。符合以上任1条件可排除试验病例。根据前期课题组研究了SO与OP的肌肉组织与骨组织差异蛋白质谱分析结果[5]，证实了SO组与OP组存在相互作用的显著差异蛋白，那么SO与OP在miRNA水平是否也同样存在差异基因调控其mRNA及下游靶基因。目前研究miRNA调控下游靶基因以揭示其潜在的调控关系，可以通过单纯的生物信息学多网站和数据库进行预测分析靶基因[6]，但往往会得到庞大的靶向调控关系结果的笼统数据，其中不可避免存在大量的假阳性结果[7]。因此，本课题组从临床收集病例，术中分别取SO组与OP组的骨组织样品，并运用Illumina高通量测序技术鉴定、筛选肌少-骨质疏松症与骨质疏松症患者骨组织microRNAs的差异表达谱[8]，以期筛查出在肌少-骨质疏松症疾病演变过程中发挥重要的调控作用的候选miRNA。

1 材料和方法

1.1 一般资料与取材

选取2017年至2019年期间在福建省老年医院和福建省立医院骨二科行手术治疗的股骨颈骨折患者中，随机选择并自愿参与本临床研究6例患者的股骨头组织(2～5 g)，SO组与OP组各选取3例；术中取一侧开路器切除的粗隆部位松质骨组织并置于液氮保存。两个组别的纳入对象均为汉族女性。本课题已通过福建省立医院的伦理委员会批准，明确患者及其近亲家属知情同意且签署了知情同意书，伦理审批号为K2019-03-034。

1.2 Illumina高通量测序建库

首先提取SO组与OP组骨组织样品的total RNA，利用small RNA的3'及5’端特殊结构直接在small RNA两端加上接头，然后反转录为c-DNA并进行PCR扩增，再通过凝胶电泳分离目的片段，验证此次实验样品RNA完整性后构建small RNA文库，随后可进行二代Illumina高通量测序。根据诺禾致源SE50测序策略进行，使用TruSeq SR Cluster Kit v3-cBot-HS(Illumia)在cBot群集生成系统上对索引编码的样本进行群集，产生簇后将文库制备物在Illumina Hiseq 2500平台和50 bp上测序。进行已知的miRNA比对检测，使用miRBase 20.0作为参考，使用改良软件mirdeep2和srna-tools-cli获取潜在的miRNA并绘制二级结构；并运用miRNA前体发夹结构的特征来预测新型miRNA，将软件miREvo和mirdeep2集成在一起，通过探索二级结构、Dicer切割位点和未注释的小RNA标签的最小自由能来预测新型miRNA。预测miRNA的靶基因是在miRanda软件中进行。对于SO组与OP组具有生物学重复的骨组织样品进行筛选miRNA的差异表达谱，使用DESeq R软件(1.8.3)对这两组进行差异表达分析，并使用Benjamini&Hochberg方法将校正的P值0.05设置为显著差异表达的阈值。对差异表达的miRNA的目标基因候选物进行基因本体论(GO)富集分析，并对目标基因候选物进行KEGG富集分析(http:// www. genome. jp/kegg/)，旨在从分子水平理解生物系统的高级功能和信号通路信息[9]。

1.3 统计学处理

2 结果

2.1 SO组与OP组比较组合的miRNA差异基因表达水平分析

利用DESeq软件对差异表达基因进行筛选，统计SO组与OP组各样本中已知的和新筛选得到miRNA的表达量，并对两组的骨组织样本间表达量相关性进行分析。SO组与OP组比较组合的差异基因统计结果显示，共筛选到3 064(包括上调和下调)个差异基因，其中上调的有1 432个，下调的有1 632，设置Pvalue <0.05且 ∣log2 Fold Change∣>0.0。①火山图(volcano plot)分析结果：通过火山图的分析可以更直观的观察到SO组与OP组之间的差异基因分布情况，横坐标X轴表示log2 Fold Change值，纵坐标Y轴则表示为-log10(Pvalue)，绿色点表示有差异且下调miRNA，红色的点则表示为有差异且上调miRNA。差异基因火山图结果提示，肌少-骨质疏松症骨组织存在差异显著表达的miRNA(图1)。②差异miRNA层次聚类分析分析结果：通过聚类分析将表达模式相近的miRNA聚集成类，进而可以识别未知miRNA的功能或已知miRNA的未知功能；此外，通过聚类分析还可以更好的认识这些同类的miRNA共同参与同一代谢过程或细胞信号通路。层次聚类分析横坐标为样品名，纵坐标为差异miRNA的每百万个映射读取的每千克碱基中的碎片归一化后的数值，分析结果为SO组与OP组组间的层次聚类分析(图2)。③GO富集结果分析：对两组差异表达 miRNA的集合进行gene ontology富集分析。其中对差异显著的miRNA进行GO富集分析，分析结果提示差异miRNA不仅在参与生物学过程(biological process，BP)、还参与了细胞组分(cellular component，CC)和分子功能(molecular function，MF)调控机体生命功能的过程。研究miRNA差异基因在 gene ontology 中的分布状况将阐明实验中SO组与OP组骨组织样本差异在基因功能上的体现。在GO富集分析结果中，将选取最显著的30个Term(keg-20)绘制柱状图(图3)、散点图(图4)。根据分析结果得到差异表达miRNA，预测其靶向的mRNA，并对差异表达miRNA的靶基因进行KEGG富集功能分析(图5)，有助于构建mRNA与miRNA调控网络，为深入研究miRNA对肌少-骨质疏松症的调控机制以及后续研究何种治疗方案对肌少-骨质疏松症患者机体细胞生物功能产生何种影响提供更多研究思路。

图1 横坐标为log2 Fold Change值，纵坐标为-log10(P value)Fig.1 The abscissa is log2FoldChange value, the ordinate is -log10 (P value)

图3 横坐标是GO Term，纵坐标则是GO Term富集的显著性水平Fig.3 The abscissa is GO Term, and the ordinate is the significance level of GO Term enrichment

图4 横坐标为注释到GO Term上的差异基因数与差异基因总数的比值，纵坐标为GO TermFig.4 The abscissa is the ratio of the number of differential genes annotated to GO Term to the total number of differential genes, and the ordinate is GO Term

图5 KEGG富集散点图Fig.5 KEGG enrichment scatter plot注：Gene Ratio 这一列表示是差异基因集；Padj 多重假设检验校正后的P值；Count注释到KEGG通路编号上的差异基因数。

2.2 qRT-PCR验证结果

内参基因U6和目的基因miR-382-3p溶解曲线呈单峰状曲线，表明样品无污染并且设计的引物特异性足够稳定，实验数据与预期结果较为吻合，miR-382-3p的实验组与对照组相对表达量见图6所示。

图6 qRT-PCR验证miR-382-3p 的表达情况与测序结果一致Fig.6 qRT-PCR verifies that the expression of miR-382-3p is consistent with the sequencing results

3 讨论

随着全球老龄化进程的加速，世界各国面临老年人群的健康问题所带来的压力越来越大，其中SO的发病率逐年增加，引起了各国医疗卫生工作者的广泛关注，SO随着年龄的增加，患者往往出现全身肌肉的减少和骨密度的下降，并且伴随而来的跌倒和骨折风险显著增加。

近年来，随着转录调控测序技术的发展与广泛运用，特别是small RNA测序技术在采用先进的Illumina测序平台后，可针对各种疾病的样本中的miRNA、siRNA、piRNA展开全面分析，既能鉴定已知sRNA、也能预测新的sRNA并预测sRNA的靶基因，为深入研究miRNA在SO疾病过程中的功能及调控机制提供有力手段[10]。

miRNAs属于一类内源性的small RNA，其可以通过降解mRNA或者抑制翻译过程来抑制靶基因的表达[11]。在肌少-骨质疏松症的研究中miRNA发挥的作用在很大程度上还是未知的。本课题组利用高通量测序手段研究SO中起关键调节作用的 miRNA及其靶基因，有望从表观遗传学层面揭示肌肉减少症患者容易罹患OP的新机制。miRNA在骨病、骨性关节炎、肌少症和OP等老年性疾病的组织样品中存在差异表达，并且miRNA在老年性疾病中的调控作用已有广泛研究[12]。根据Illumina高通量测序结果筛选出可能参与调节SO疾病发生演变过程的候选miRNA，并充分利用GO和KEGG数据库的生物信息可以评估显著差异的miRNA在SO发生、发展中所发挥的调控作用。

本课题组在研究中利用Illumina高通量测序技术对SO组与对照组OP进行对比，差异表达显著的miRNA有22个，其中15个发生了上调(表1)，7个(表2)发生下调；与OP对照组相比，在SO组骨组织中存在显著高表达的miRNA，其中hsa-miR-382-3p、hsa-miR-27a-5p、hsa-miR-1226-5p和hsa-miR-3934-5p是发生上调较显著的miRNA，显著低表达的是hsa-miR-451a。研究结果提示在SO患者骨组织中存在部分差异表达的miRNA，其可能参与调控肌少症罹患OP的疾病过程发生发展过程。

表1 显著上调的miRNATable 1 Significantly up-regulated miRNAs

表2 显著下调的miRNATable 2 Significantly down-regulated miRNAs

本研究中发现的特征性候选基因hsa-miR-382-3p在膝骨性关节炎存在差异表达，目前国内外研究发现miR-382-3p可能通过直接靶向CX43通过TLR4/MyD88/NF-κB信号通路参与OA[13]；Heilmeier等[14]在糖尿病性骨病和绝经后骨质疏松症研究中发现miR-382-3p的差异表达对机体成骨、脂肪形成和细胞增殖的存在调控作用。本研究中，肌少-骨质疏松症患者的miR-27a-5p表达水平较显著上调，有研究表明miR-27a-5p通过抑制Atg7可减轻缺氧诱导的大鼠心肌细胞损伤[15]，其也可能在肌少-骨质疏松症的病理生理过程中发挥重要作用。hsa-miR-451a在肌少-骨质疏松症组中表达显著下调，Lu等[16]研究发现miR-451a在骨质疏松症中存在差异表达，其可能通过BMP6信号传导促进成骨分化并抑制骨量丢失进程；Karvande等[17]研究报道甲状旁腺素促进葡萄糖依赖性miR-451a表达，刺激成骨细胞分化，其机制是甲状旁腺素可以通过PI3K-mTOR-AKT轴抑制AMPK磷酸化，从而防止作用于miR-451a启动子区域的八聚体结合转录因子1(OCT-1)的磷酸化和失活，AMPK活性的调节会影响成骨细胞中miR-451a的水平。

本课题组在后续基础研究中将筛选出可能调控肌少-骨质疏松症的候选miRNA，并通过预测分析以确定候选miRNA所作用的靶基因，在细胞与动物水平上深入研究候选的miRNA与靶基因的表达水平及其调控机制。