运用高通量RNA测序探讨绝经后骨质疏松与免疫功能的关系

陶有娣 王 潇 陈 健

(厦门大学医学院，福建 厦门 361000)



运用高通量RNA测序探讨绝经后骨质疏松与免疫功能的关系

陶有娣 王 潇 陈 健1

(厦门大学医学院，福建 厦门 361000)

目的 探讨绝经后骨质疏松与免疫功能的关系。方法 选取6例符合要求的绝经后病人，3例骨质疏松患者为实验组，3例非骨质疏松者为对照组。双能X线吸收仪测定骨密度。空腹抽取肘静脉血适量，提取外周血单个核细胞RNA。 Nanodrop分光光度计及琼脂糖电泳对RNA进行质检，构建RNA文库。用 Illumina Hiseq2000平台测序，对数据进行统计，寻找差异表达基因，并且对这些差异表达基因进行功能富集分析。结果 根据差异分析P<0.05阈值筛选两组的差异基因，这些差异表达基因功能富集显示，大多集中在与免疫相关的信号通路上，如防御反应(IL-15,IL-1 RAP,TNFAIP6,TNFRSF4,IL-2RA,IL-8等)，免疫反应(CD8A,CD8B,IL-15,TNFRSF13C CD79B,CD14等)，炎症反应(IL-15,CD24,TLR10,TNFAIP6等)，白细胞活化(IL-8,CD8A,CD8B等)和淋巴细胞活化(CD3G,CD8A,CD8B,THEMIS,KIR2DS4,KIR3DL1等)。结论 绝经后骨质疏松的发生与患者的免疫功能改变有一定的关系。

绝经后骨质疏松；高通量RNA测序；免疫功能

绝经后骨质疏松属于原发性骨质疏松，源于绝经期卵巢功能的退化，成骨细胞介导的骨形成与破骨细胞介导的骨吸收失平衡，导致系统性骨骼紊乱，骨密度降低，骨组织微环境退化及骨折脆性与敏感性增加〔1〕。有研究表明，绝经后骨质疏松发生过程中，雌激素缺乏导致T细胞活性增加，活化的T细胞分泌的一些炎性因子增加，如肿瘤坏死因子(TNF)-α、核因子(NF)κB受体活化因子受体的配体(RANKL),进一步提高了破骨细胞活性〔2〕。在雌激素缺乏的情况下，活化的T细胞与TNF的来源最密切。也有报道表明T细胞缺陷的裸鼠对卵巢切除后的骨量丢失具有保护作用〔3〕。说明T淋巴细胞在绝经后骨质疏松的发生发展中占重要的一部分。

在骨质疏松的发生发展中，骨保护素(OPG)/RANKL/RANK 信号通路起重要作用〔4〕。与RANKL竞争结合RANK,调节破骨细胞的生成，而B淋巴细胞产生OPG占骨髓的64%，骨保护素负性调节破骨细胞的生成。研究表明，B淋巴细胞缺乏的小鼠表现出持续的骨质疏松及骨髓中OPG水平下降，这种情况可以被B淋巴细胞重建所补救〔5〕。

高通量测序是近年来一个热门的技术，主要特点是测序通量高，测序时间和成本显著下降〔6〕。本研究运用高通量RNA测序技术，从基因水平进一步探讨骨质疏松与免疫功能间的关系。

1 材料和方法

1.1 实验对象 选取6例绝经后妇女，符合以下条件：(1)年龄50～60岁，绝经2年以上、10年以下，自愿加入该实验研究，所有受试者均来自厦门大学附属中山医院门诊病人。(2)骨质疏松组：符合1994年世界卫生组织确立的骨质疏松诊断标准，即双能X线吸收仪测定股骨颈或腰椎骨密度T值<-2.5标准差。(3)健康对照组：双能X线吸收仪测定股骨颈或腰椎骨密度T值>-1标准差。(4)以下情况两组均予排除：近期或以前经过抗骨质疏松治疗(双磷酸盐类，特立帕肽，雷尼酸锶，雷洛昔芬)，近1个月使用过治疗骨质疏松中药者，近3个月内采用激素替代治疗和使用降钙素(密钙息和益钙宁等)，1年内应用唑来膦酸注射液者；以前或现在患影响骨或免疫系统疾病，包括早绝经(绝经年龄<40岁)、甲状旁腺功能亢进、未经治疗的甲状腺功能亢进、骨软化症、类风湿关节炎、痛风、多发性骨髓瘤、脊柱结核、肿瘤患者、人类免疫缺陷病毒(HIV)感染、乙型或丙型慢性肝炎、慢性肾功能或肝衰竭，或严重畸形、残废、丧失劳动力；以前或现在采用影响骨与免疫系统的疗法：长期类固醇治疗，免疫抑制治疗(化学疗法，甲氨蝶呤，环孢素，TNF-α阻断剂等)；合并新鲜骨折及严重多处骨折患者。

1.2 主要试剂、仪器 Trizol(美国Ambion 公司)；人淋巴细胞分离液CBD(中国斯坦姆生物公司)；RNA建库TruSeq® RNA LT Sample Prep Kit v2(美国Illumina公司)，测序试剂 TruSeq PE Cluster Kit v3-cBot-HS，TruSeq SBS Kit v3-HS(200-cycles)(美国 illumina公司)；捕获磁珠AmpureBeads(美国Beckman公司)；Qubit定量 Quant-iTTMPicoGreen® dsDNA Assay Kit (美国life公司)Hisq2000高通量测序仪；CBOT簇生成仪(美国Illumina公司)；Nanodrop(美国Thermo公司)；Qubit(美国Invitrogen公司)；普通PCR仪(美国Life公司)；电泳仪(中国天能公司)；凝胶成像系统(中国天能公司)

1.3 实验步骤

1.3.1 提取RNA 空腹抽取6例研究对象肘静脉血适量，加入等体积的人淋巴细胞分离液CBD，高速离心(2 000 r/min)分离出外周血单个核细胞，生理盐水洗涤离心1次，弃上清，加1 ml TRIzol试剂混匀，再经纯化提取外周血单个核细胞RNA，测定浓度。

1.3.2 RNA质检 用1 μl Nanodrop仪器(美国Thermo公司)对提取的RNA进行质检。

1.3.3 构建测序文库 将RNA分子随机打断并使用随机引物将RNA片段反转录成cDNA分子，然后两端加上特定的接头。

1.3.4 锚定桥接 测序过程都在Illumina公司提供的流动槽中进行。一个标准的流动槽含有8个通道,每个通道的内表面有无数被固定的单链接头。上述步骤得到的带接头的DNA 片段变性成单链后与测序通道上的接头引物结合形成桥状结构，用于接下来的桥式PCR扩增。加入每个通道的是之前构建好的DNA片段文库。

1.3.5 预扩增 添加未标记的dNTP和普通Taq酶进行固相桥式PCR扩增，单链桥型待测片段被扩增成为双链桥型片段。通过变性释放出互补的单链，锚定到附近的固相表面。通过不断循环，将会在流动槽的固相表面上获得上百万条成簇分布的双链待测片段。

1.3.6 单碱基延伸测序 在测序的流动槽中加入四种荧光标记的dNTP、DNA聚合酶以及接头引物进行扩增，在每一个测序簇延伸互补链时，每加入一个被荧光标记的dNTP就能释放出相对应的荧光，测序仪通过捕获荧光信号，并通过计算机软件将光信号转化为测序峰，从而获得待测片段的序列信息。计算机软件将同一位置坐标上检测到的碱基依次连接，形成读段。

1.4 统计学方法 应用 SPSS13.0 软件行t检验。

2 结 果

2.1 采用Tophat软件对预处理后样本序列与参考基因组序列进行全基因组序列比对，计算序列的全基因 组覆盖情况和计算RNA 的表达量，使用基因组为hg19，比对参数使用默认。结果显示各个样本的定位百分比均在90%以上，即样本RNA合格，可以用于下一步的测序步骤。见表1。

表1 六个样本的质检比对概述

项目骨质疏松组185R186R188R非骨质疏松组190R192R334D原始读段数目20,541,05419,148,91119,861,15519,060,99718,987,93316,223,061定位百分比94.61%94.16%94.27%94.14%93.82%92.59%丰度5.42%4.25%5.06%6.81%4.57%2.92%未定位百分比5.39%5.84%5.73%5.86%6.18%7.41%

原始读段数目：各个样本原始读段数目；定位百分比：定位上的读段数占原始读段总数的百分比；丰度：定位在高丰度基因(例如核糖体rRNA)上的读段数占原始读段总数的百分比；未定位百分比：未定位上的读段数占原始读段的总数百分比

2.2 基因体富集分析(GO) 根据P<0.05筛选出两组的差异表达基因，将获得的差异基因列表或者转录本列表通过DAVID (http://david.abcc.ncifcrf.gov/)进行GO(Gene Ontology)富集分析，研究差异基因在GO中的分布状况将阐明实验中样本差异在基因功能上的体现。GO分为三类：生物学过程(biological process)，分子功能(molecular function)和细胞组成(cellular componen)。结果显示这些差异表达基因在生物功能方面，大多与免疫功能相关，如防御反应、免疫反应、炎症反应、白细胞活化等。见表2。

2.3 全基因组及代谢途径数据库(KEGG)通路(Pathnay)富集分析 在生物体内，不同基因相互协调行使其生物学功能，通过Pathway 显著性富集分析能够确定差异基因参与的最主要的生化代谢途径和信号转导途径。KEGG (Kyoto Encyclopedia of genes and Genomes，www.genome.jp/kegg/)是有关Pathway的主要公共数据库。将获得的差异基因列表或者转录本列表通过DAVID (http://david.abcc.ncifcrf.gov/)进行KEGG通路富集分析。结果显示差异基因富集在与免疫相关的信号通路上。 见表3。

表2 通过P<0.05筛选得到绝经后骨质疏松组与非骨质疏松组差异表达基因，对差异表达基因进行GO功能富集

目录名称P值基因生物功能防御反应1.60E-07S100A8,IL-15,IL1RAP,LTF,TNFAIP6,CD40LG,TNFRSF4,CD24,IL2RA,IL-8,RNASE3,S100A12,CD19,P2RX1,CD14…免疫反应4.46E-07CD8A,CD8B,IL-15,LTF,DPP8,PRG2,CD40LG,TNFRSF4,CD24,TLR10,IL-2RA,IL-8,TNFRSF13C,CD79B,CD14…炎症反应2.16E-05CCL3,S100A8,IL-15,TNFRSF4,FOS,MEFV,CD24,TLR10,IL-2RA,IL-8,S100A12,TNFAIP6,CD14…白细胞活化5.83E-05IL-8,CD8A,CD8B,IL-15,TNFRSF4,IFNAR1,AZU1,CD40LG,LAX1,CD2,CD24,KIR3DL1,DPP4…淋巴细胞活化2.38E-04CD3G,CD8A,CD8B,THEMIS,KIR2DS4,IL-15,TNFRSF4,IFNAR1,CD40LG,LAX1,MS4A1,CD2,CD24,KIR3DL1,DPP4…分子功能氧气运输活性4.21E-06HBM,HBQ1,HBG1,HBA2,HBG2,HBA1,HBB,HBD细胞因子结合1.70E-05IL-2RA,CMKLR1,TNFRSF25,TGFBR1,ELANE,CXCR1,CXCR2,TNFRSF4,IL-11RA,IFNAR1,CCR6,CCR5,IL1RAP,TGFBR3,IL5RA,THBS1…氧气结合2.94E-05HBM,HBQ1,HBG1,HBA2,CYP4F3,CYP2E1,HBG2,HBA1,HBB,SOD2HBD…碳水化合物结合0.0011838ENPP3,FCER2,CLEC17A,CD69,CD22,CLEC4C,CD24,PRG2,CLC,TNFAIP6,SI-GLEC1,KLRG1,CLEC12A,IGF2R,CLEC12B,TGFBR3,CD14…血红素结合0.002341HBM,CYP2D6,FADS3,HBA2,HBA1,CYP2E1,CYP27A1,HBQ1,MPO,CYP4F3,HBG1,HBG2,HBB,HBD…细胞组成质膜部分1.42E-10SEC31B,SLC22A18,RALB,SLC4A1,KLRD1,KIR2DL1,SYT11,NOTCH2,SLC4A10,SYTL2,TP53I13,SORL1,SYNGR1,KCNIP2,PILRB,RAP2A…细胞表面3.48E-04CD69,CD2,RC3H2,CD22,MS4A2,FCER1A,IL-2RA,ELANE,TNFRSF13C,NOTCH2,CD19,CCR5,TGFBR3,AMOT,CD226,TP53I13,CTSG…胞外区域0.004373UTS2,MMP8,PGLYRP1,CBLN3,RETN,LYNX1,COLQ,CD40LG,CCL3FCER2,PI3,TGFBR3,CD8A,CD8B,MMP28,MMP25,LNPEP,CAMP,CA2…

表3 根据P<0.05,筛选得到绝经后骨质疏松组与非骨质疏松组差异表达基因，对其进行Pathway显著性分析

KEGG通路基因数P值基因抗原呈递通路124.58E-04HLA-DQB1,CD8A,CD8B,HLA-DRB5,KIR2DS4,KIR2DL1,HS-PA1A,KIR2DL3,KLRD1,HLA-DOB,KIR3DL1,KIR2DL4自然杀伤细胞介导的细胞毒性通路158.83E-04PIK3CG,KIR2DS4,IFNAR1,TNFRSF10C,PIK3CA,KIR2DL1,PPP3CA,KIR2DL3,NFATC2,PIK3R3,SH2D1B,FCGR3B,KLRD1,KIR3DL1,KIR2DL4趋化因子信号通路160.008356457PIK3CG,CCL3,GNAI3,ROCK1,IL-8,ROCK2,CXCR1,CXCR2,FOXO3,CCR6,CCR5,ADCY9,PIK3CA,PIK3R3,PLCB1,GNG7B细胞受体信号通路90.010366428PIK3CG,FOS,CD19,CD22,PIK3CA,CD79B,PPP3CA,PIK3R3,NFATC2T细胞受体信号通路100.030121964PIK3CG,FOS,CD3G,CD8A,CD8B,CD40LG,PIK3CA,PPP3CA,PIK3R3,NFATC2

基因数：富集到该KEGG通路的基因数

3 讨 论

骨骼与免疫系统在细胞起源、空间分布、功能及调控等方面存在许多相似性，免疫与骨骼系统有一些共同的调节因子，如细胞因子、转录因子及受体〔7〕。Arron等〔8〕于2000年第一次使用骨免疫学描述免疫细胞与骨骼系统之间的相互作用。这两者间的密切联系使得在生理及病理条件下均相互作用，尤其是一个系统的病态激活影响另一个系统，如类风湿关节炎中，免疫系统的异常活化影响骨重建导致病理性骨吸收〔9〕。而且，在慢性炎症状态，骨形成与骨吸收之间的平衡偏向于破骨介导的骨吸收，尤其是近年来对OPG/RANKL/RANK〔10〕系统的深入研究，使人们对骨与免疫系统之间的相互作用机制有更深的认识。RANK属于TNF超家族，很多的组织及细胞均表达RANK，如骨骼肌、肝脏、乳腺，单核巨噬细胞系、B、T淋巴细胞、树突细胞等。RANK的活化可刺激破骨前体细胞向成熟的破骨细胞分化及成熟的破骨细胞活化。RANKL 同样属于TNF超家族。骨组织及骨外组织细胞均可表达RANKL，如成骨细胞、骨髓基质细、骨细胞、乳腺等〔11〕。而且，免疫系统的很多细胞也可表达RANKL，如单核细胞、中性粒细胞、树突细胞、B、T淋巴细胞。已知这些免疫细胞可以产生促炎性细胞因子，通过加强RANK/RANKL信号通路导致骨破坏〔10〕。

高通量RNA测序，在测序过程中，由于RNA随机片段化和采用随机引物进行反转录等，所得cDNA片段较均匀地取自各转录本。因此，RNA测序得到的读段可以视为对转录组随机打断后的随机采样，这是人们使用统计模型来对转录组和基因表达作参数估计分析的理论基础〔12〕。

有研究表明，骨质疏松性骨折的发生与很多慢性炎症性疾病相关，如类风湿关节炎、强直性脊柱炎、克罗恩病、溃疡性结肠炎、慢性阻塞性肺疾病及一些病毒感染〔13〕。Breuil等〔14〕观察26例伴有骨质疏松性骨折的绝经后骨质疏松女性及24例年龄和雌激素水平匹配的健康对照组的免疫细胞的表型和功能特性。他们首次观察到B淋巴细胞数目减少，由此推测B细胞的这种改变可能与体内骨髓微环境改变有关，并且独立于年龄和雌激素状态。本实验GO富集分析显示两组差异表达基因在生物功能方面大多与免疫功能相关，如防御反应、免疫反应、炎症反应、白细胞活化等。KEGG通路分析结果表明，绝经后骨质疏松的发生与免疫功能改变有一定的关系。由此推测绝经后，患者自身免疫功能与骨微环境相互作用，引起骨形成与骨吸收之间的平衡偏向于破骨介导的骨吸收。

目前两种主要治疗骨质疏松的药物是抑制骨吸收和刺激新骨形成。抑制骨吸收药物包括双膦酸盐类，如依替膦酸钠、阿仑膦酸钠、利塞膦酸钠和唑来膦酸;雌激素与选择性雌激素受体调节剂——雷洛昔芬;鲑鱼降钙素和狄诺塞麦。目前唯一的刺激新骨形成药物是特立帕肽〔15〕。本实验从基因表达水平研究骨质疏松与免疫功能的关系，发现两组绝经后病人的差异表达基因大部分集中在与免疫功能相关的信号通路上。

1 Willson T,Nelson SD,Newbold J,etal.The clinical epidemiology of male osteoporosis:a review of the recent literature〔J〕.Clin Epidemiol,2015;7:65-76.

2 Hotokezaka H,Sakai E,Ohara N,etal.Molecular analysis of RANKL-independent cell fusion of osteoclast-like cells induced by TNF-alpha,lipopolysaccharide,or peptidoglycan〔J〕.J Cell Biochem,2007;101(1):122-34.

3 D'Amelio P,Grimaldi A,Di Bella S,etal.Estrogen deficiency increases osteoclastogenesis up-regulating T cells activity:a key mechanism in osteoporosis〔J〕.Bone,2008;43(1):92-100.

4 Lacey DL,Boyle WJ,Simonet WS,etal.Bench to bedside:elucidation of the OPG-RANK-RANKL pathway and the development of denosumab〔J〕.Nat Rev Drug Discov,2012;11(5):401-19.

5 Ma M,Chen X,Lu L,etal.Identification of crucial genes related to postmenopausal osteoporosis using gene expression profiling〔J〕.Aging Clin Exp Res,2016；28(6)：1067-74.

6 Shendure J,Ji H.Next-generation DNA sequencing〔J〕.Nat Biotechnol,2008;26(10):1135-45.

7 Walsh MC,Kim N,Kadono Y,etal.Osteoimmunology:interplay between the immune system and bone metabolism〔J〕.Annu Rev Immunol,2006;24(1):33-63.

8 Arron JR,Choi Y.Bone versus immune system〔J〕.Nature,2000;408(6812):535-6.

9 Takayanagi H.Mechanistic insight into osteoclast differentiation in osteoimmunology〔J〕.J Mol Med (Berl),2005;83(3):170-9.

10 Herman S,Kronke G,Schett G.Molecular mechanisms of inflammatory bone damage:emerging targets for therapy〔J〕.Trends Mol Med,2008;14(6):245-53.

11 Narducci P,Bareggi R,Nicolin V.Interaction of human recombinant tumor necrosis factor-related apoptosis-inducing ligand and osteoprotegerin could contribute to enhancement of the erosive processes induced by human synovial cells〔J〕.J Rheumatol,2009;36(8):1837-9.

12 Tripathy S,Jiang RH.Massively parallel sequencing technology in pathogenic microbes〔J〕.Methods Mol Biol,2012;835(2):271-94.

13 Pietschmann P,Mechtcheriakova D,etal.Immunology of osteoporosis:a mini-review〔J〕.Gerontology,2016;62(2):128-37.

14 Breuil V,Ticchioni M,Testa J,etal.Immune changes in post-menopausal osteoporosis:the Immunos study〔J〕.Osteoporos Int,2010;21(5):805-14.

15 Johnson GL.Denosumab (prolia)for treatment of postmenopausal osteoporosis〔J〕.Am Fam Physician,2012;85(4):334-6.

〔2016-01-19修回〕

(编辑 袁左鸣)

Relationship between postmenopausal osteoporosis and immune function

TAO You-Di, WANG Xiao, CHEN Jian.

College of Medicine，Xiamen University，Xiamen 361000，Fujian，China

Objective To detect the relationship between postmenopausal osteoporosis and immune function.Methods The study involved 6 women, including 3 postmenopausal osteoporosis and 3 postmenopausal non-osteoporotic women. BMD measurement was performed by dual-energy X-ray absorptiometry. After an overnight fast, morning blood was drawn from the candidates. Total RNA was extracted from peripheral blood mononuclear cells. RNA quality and purity were confirmed with a Nanodrop spectrophotometer and electrophoresis. RNA library was constructed. With Illumina Hiseq 2000 sequencing platform, dates and searched for differentially expressed genes(DEGs)were collected. Afterwards, functional and pathway analyses were conducted.Results The DEGs were mainly enriched in GO term of defense response (e.g.,IL-15,IL-1 RAP,TNFAIP6,TNFRSF4,IL-2 RA,IL-8),immune response(e.g., CD8A,CD8B,IL-15,TNFRSF13C,CD79B,CD14),inflammatory response(e.g., IL-15,CD24,TLR10,TNFAIP6),leukocyte activation(e.g., IL-8,CD8A,CD8B),lymphocyte activation(e.g., CD3G,CD8A,CD8B,THEMIS,KIR2DS4,KIR3DL1).Conclusions The postmenopausal osteoporosis is associated with changes of immune function.

Postmenopausal osteoporosis; High-throughput RNA sequencing; Immune function

国家自然科学基金资助项目(No.81272168);福建省医学创新课题资助项目(No.2012-CXB-32)

1 厦门大学附属中山医院康复医学科

陈 健(1963-)，男，博士，主任医师，硕士生导师，主要从事骨质疏松研究。

陶有娣(1989-)，女，在读硕士，主要从事骨质疏松研究。

R3

A

1005-9202(2016)23-5812-04；

10.3969/j.issn.1005-9202.2016.23.017