间充质干细胞成骨分化早期染色质开放区域的动态变化*

刘建云， 张杰， 马百成， 吴萍， 熊建军

间充质干细胞成骨分化早期染色质开放区域的动态变化*

刘建云， 张杰， 马百成， 吴萍， 熊建军△

（九江学院基础医学院，江西 九江 332000）

在全基因组水平观察骨髓间充质干细胞（MSCs）成骨分化早期染色质开放区域的动态变化。原代培养人骨髓MSCs至第6代，成骨诱导剂分别刺激0、3、5和7 d。收集各组细胞，运用ATAC-seq技术（基于高通量测序的转座酶可接近染色质分析）检测染色质开放区域，并对各组间染色质开放区域数量、开放区域的DNA基序以及邻近基因功能进行生物信息学分析。在MSCs成骨定向分化早期的不同时点，细胞染色质开放区域数量发生显著动态变化；各组功能性区域peak分布主要富集在转录起始点附近； motif分析显示各组细胞间参与转录调控的活化DNA发生细微改变；新增染色质开放区域的生物学过程主要集中在各种GTP酶活性的调节、细胞外基质组成、Wnt信号通路等；而新增的基因显著性富集信号通路主要集中在Rap1信号通路、黏着斑、黏着连接等。在骨髓MSCs成骨定向分化早期，各时段的染色质开放区域发生了显著的特征性改变。

间充质干细胞；成骨分化；ATAC-seq；生物信息学

骨质疏松症是全球性公共健康问题，以骨量减少、微观结构破坏、骨脆性增加和易发生骨折为显著的临床特征。经典观点认为，骨质疏松的发生与破骨细胞产生增多、活性增强，致使骨量丢失与骨质量下降有关［1］。新近研究证实，来源于骨髓的间充质干细胞（mesenchymal stem cells， MSCs）定向分化为成骨细胞（osteoblast， OB）和脂肪细胞（adipocyte， AD）的比例失衡（成骨细胞减少、脂肪细胞增多），是骨质疏松发生发展的另一关键因素［2-4］。成骨诱导剂（地塞米松联合β-甘油磷酸钠、维生素C）是体外诱导MSCs向成骨细胞分化最经典、最常用的方法。经过21~28 d诱导， MSCs由骨祖细胞经前成骨细胞分化为成骨细胞，表达特异性成熟标志物，如碱性磷酸酶（alkaline phosphatase， ALP）、I型胶原蛋白及骨钙素（osteocalcin）等［5］。利用这一诱导模型，一些调控成骨分化的重要信号分子和通路得以发现和阐明。在MSCs的分化过程中伴随着众多基因的转录激活或抑制，表观遗传在转录水平动态地调控相关基因的表达，对细胞功能发挥了至关重要的作用［6］。基于高通量测序的转座酶可接近染色质分析（assay for transposase-accessible chromatin with high-throughput sequencing， ATAC-seq）是近年来兴起的用于研究染色质开放性的表观遗传学工具，通过获得染色质上开放区域的位置和活跃的调控序列，有助于在全基因组范围内鉴定激活的转录因子及其作用规律［7-8］。本项目以MSCs定向分化为成骨细胞的早期过程为研究对象，运用ATAC-seq技术分析其中的染色质开放性变化，为研究成骨定向分化的调控机制及鉴定有效靶基因提供了新的切入点。

材料和方法

1　细胞培养和处理

细胞来源于本实验室保存的人源性原代培养骨髓MSCs。培养环境： 37℃， 5% CO2。每隔48 h传代１次。选择第6代MSCs进行成骨诱导，诱导剂成分为：DMEM培养液中含10%胎牛血清、100 mmol/L地塞米松、10 mmol/L甘油磷酸钠和50 µg/L维生素C。成骨诱导剂分别刺激0 d （MSC-0d）、3 d （OB-3d）、5 d （OB-5d）和7 d （OB-7d）后，收集细胞，在4℃环境下以500×离心5分钟去除上清液，随后以50 µL冰冷PBS洗涤1次，去除上清液，再以50 µL 冰冷lysis buffer悬浮细胞，立即离心10 min；去除上清液，进行后续转座反应。

2　转座反应与纯化

细胞置于冰上，使用转座反应体系（25 μL 2× reaction buffer from Nextera kit、2.5 μL Nextera Tn5 transposase from Nextera kit和22.5 μL nuclease-free H2O）悬浮细胞核。37℃孵育30 min；使用Qiagen MinElute PCR Purification Kit纯化DNA； 10 µL elution buffer洗脱产物；取转座后DNA片段配制成PCR反应体系（10 µL DNA、10 µL nuclease-free H2O、2.5 µL PCR Primer 1、2.5 µL Barcoded PCR Primer 2和25 µl NEBNext High-Fidelity 2× PCR Master Mix）。反应条件：72℃延伸5 min； 98℃变性30 s （1循环）；98℃变性10 s， 63℃退火30 s， 72℃延伸1 min （共5循环）； 72℃延伸5 min； 4℃冷却。最后使用Qiagen MinElute PCR Purification Kit纯化PCR产物［9］。库检合格后进行Illumina HiSeq测序。

3　数据分析

测序数据下机之后，获得原始测序序列（raw data），在有相关物种参考序列或参考基因组的情况下，通过标准流程进行生物信息分析。

首先使用FastQC软件来查看测序数据的质量情况，包括测序错误率分布和碱基含量分布。下机的原始序列进行去接头处理，使用BWA软件将得到的有效数据（clean data）比对到参考基因组hg38_genecode上，该软件采用BWT算法建立索引的短序列快速比对，常用于将短序列比对到大型基因组上［10］。随后比对分析后的bam文件作为输入文件，使用MACS2软件进行Call Peak（使用统计学方法计算出参考基因组中比对上的Reads显著富集的区域），筛选阈值为q value<0.05［11］。对每个peak区域从上、下游两个方向分别拓宽200 bp提取出DNA序列，使用HOMER软件的findMotifsGenome.pl工具对DNA序列预测motif，然后将预测的motif与各大数据库（HOMER和JASPAR）中已有的motif数据进行匹配，找到已发表的ChIP-seq实验中的motif和相应的转录因子［12］。基因附近信号分布图的分析使用deeptools软件，对所有基因根据基因长度进行归一化处理，计算在转录起始位点（transcription start site， TSS）上游3 kb到转录终止位点（transcription termination site， TSS）下游3 kb范围内每个位点上所有基因的平均信号值，并绘制成曲线图［13］。最后利用DAVID （Database for Annotation， Visualization and Integrated Discovery）对染色质开放区域所关联的基因进行基因本体论（Gene Ontology， GO）富集分析［14］；基于KEGG数据库对peak邻近基因进行信号通路富集分析［15］。

结果

1　细胞染色质开放区域的鉴定

在成骨诱导剂作用下，体外培养的MSCs逐渐出现形态变化，至14 d，分化为能分泌骨基质的成骨细胞（图1）。我们使用ATAC-seq来检测早期不同分化时点（0、3、5和7 d）染色质开放区域的动态变化。测序下机的原始数据经筛选后，采用BWA软件比对，各组Reads比对率均高于90%； Reads在基因区域的信号分布均富集在TSS附近（图2A）。使用MACS软件对各组Reads进行Call Peak，在MSC-0d组筛选出110 369个Reads显著富集的区域（peak），在OB-3d组筛选出113 785个peak，在OB-5d组筛选出110 267个peak，在OB-7d组筛选出68 216个peak，表明成骨分化的第7天，细胞转录活性区域显著减少（图2B）。功能性区域分布显示，位于启动子与TSS之间的peak比例大约在10%~15%之间变动，但是各组之间也存在细微的区别。值得注意的是，在成骨诱导的第7天，位于启动子与TSS的染色质开放区域的比例在各组中最高（15.98%），尽管这一时期整个peaks数最少（图2C）。

Figure 1.　Alizarin red staining of MSCs at different time points of differentiation (scale bar=20 μm). A: 0 d; B： 7 d; C： 14 d.

Figure 2.　Identification of the chromatin accessibility. A: signal distribution map near the gene （the horizontal axis represents the normalized gene range coordinates, and the vertical axis represents the average signal value of the site）; B: the number of peaks in each group; C: peak distribution in functional areas between MSC-0d and OB-7d groups.

2　peak序列motif分析

HOMER软件对4组细胞染色质开放区域的DNA序列进行motif分析，显示各组细胞中结合最多的motif均是亮氨酸拉链（basic leucine zipper， bZIP）转录因子家族的成员，如、、、、、-等，但是每组之间的排列顺序又略有差异（表1）。

表1　各组细胞染色质开放区内motif排序表

随后我们对各分化组与MSC-0d组之间的差异motif进行分析，显示在OB-3d组与MSC-0d组之间上调数量变化最显著的motif是基因家族、基因家族和等；在OB-5d组与MSC-0d组之间新增数量变化最为显著的motif是基因家族、基因家族和等；在OB-7d组与MSC-0d组之间， motif则发生了更为显著的变化，数量增加最为明显的是、、、、等（图3A）。RUNX2作为调控成骨早期分化的关键转录因子，我们重点追踪了其结合序列开放性的动态变化。结果显示，在成骨诱导的第3天和第5天，结合motif的开放数量在整个细胞中的排序中前移，但在第7天显著下降（图3B），表明的转录活性在第5天后逐渐弱化。

Figure 3.　Motif analysis in each group. A: the main different motifs in each group; B: the dynamic ranks of RUNX2 motif during osteogenic differentiation.

3　peak邻近基因功能GO富集分析

GO富集分析显示， MSC-0d、OB-3d、OB-5d和OB-7d组细胞染色质开放区域所关联基因的主要生物学过程（biological process， BP）并无显著差异，大多是蛋白磷酸化和代谢途径等，但各组间细胞差异染色质开放区域所关联的基因BP有所不同。在OB-3d组与MSC-0d组之间，上调最为显著的BP有正性调节GTP酶活性（positive regulation of GTPase activity）、正性调控Rho GTP酶活性（positive regulation of Rho GTPase activity）、正性调控Rac GTP酶活性（positive regulation of Rac GTPase activity）、细胞外基质组成（extracellular matrix organization）、正性调控Wnt信号通路（positive regulation of Wnt signaling pathway）等；在OB-5d组与MSC-0d组之间，上调最为显著的BP与在OB-3d组与MSC-0d组之间相似；但在OB-7d组与MSC-0d组之间，差异BP发生显著变化，上调最显著的BP为细胞黏附（cell adhesion）、细胞外基质组成、细胞-基质黏附（cell-matrix adhesion）等（图4）。

Figure 4.　Differential GO enrichment analysis of peak adjacent genes in each group.

4　peak邻近基因功能信号通路富集分析

信号通路富集分析显示， MSC-0d、OB-3d、OB-5d和OB-7d组细胞染色质开放区域所关联的主要信号通路都包括磷酸化、转录调控和代谢过程等。但各组细胞之间差异信号通路有所不同。在OB-3d组与MSC-0d组之间，上调较为显著的信号通路有Rap1信号通路、黏着斑通路、炎症介质对TRP通道的调节等。在OB-5d组与MSC-0d组之间，上调较为显著的信号通路有Rap1信号通路、黏着连接和黏着斑通路等。在OB-7d组与MSC-0d组之间，上调较为显著信号通路有PI3K-Akt信号通路、Rap1信号通路、黏着斑通路等（图5）。

Figure 5.　Differential pathway enrichment analysis of peak adjacent genes in each group.

讨论

骨髓MSCs是在骨髓中发现的非造血干细胞群，其显著特征在于具有强大的自我更新能力与多向分化潜能［16］。在分化早期，骨髓MSCs定向分化为特异性祖细胞是决定谱系定向的关键环节，因此，我们将观察的染色质开放区的时点设定为成骨诱导后的第3天、第5天和第7天。

MSCs向成骨细胞分化过程中伴随着众多基因的转录激活或抑制。在调控基因转录的启动子区和增强子区，“开放”和“封闭”染色质分别代表单个基因的转录激活和抑制状态。“开放”的染色质区域显示较低的核小体密度甚至无核小体，并与修饰的核心组蛋白的特定组合；相反，转录抑制的基因通常在“封闭”染色质结构域内［17-18］。ATAC-seq技术检测的内容为基因组中的“开放”区域，即与转录因子结合调控基因表达的序列，对于了解组织分化和分化的分子机制提供了新的手段。本研究显示，来源于骨髓间充质的第六代干细胞中有超过十万个染色质开放区域，并且主要富集在转录起始点附近，表明MSCs在体外培养条件下就具有较高的基因转录活性。但是由于体外培养环境与体内生长环境的差异，我们尚不能以此来推测MSCs在体内生长时的染色质开放状态。在成骨诱导剂作用第3天和第5天，染色质开放区域的数量趋于稳定，而到了第7天，显著下降，表明具有转录活性的基因明显减少。我们推测成骨细胞分化的第7天可能是细胞功能的分水岭：在成骨诱导剂的作用的前7 d，MSCs以定向分化为主，转录因子活跃，参与的基因数量众多；而7 d之后，细胞趋于成熟，功能上以分泌细胞外基质和矿化为主，所以激活基因数量显著下降。这一推测需要相应的RNA-seq进行验证。

除了染色质开放数量，各组间功能性区域peak分布也有所不同。在成骨诱导的第3天和第5天，功能性区域peak分布在启动子与转录起始点之间的染色质开放区域的比例微小变动，但在第7天，这一比例明显升高，提示细胞功能发生了明显变化。值得注意的是，虽然基因间检测到的peak较多，但是由于基因间区在基因组上占比极大，所以基因间区上微弱的信号一般不是真的调控因子结合位点。

与染色质开放数量趋势相似的是各组间motif变化。有意思的是，各组细胞中占主导地位的转录因子结合序列均是bZIP基因家族成员（如、、和-等）。bZIP是碱性亮氨酸拉链（basic leucine zipper）蛋白大家族，普遍存在于动植物及微生物中，通过二聚化亮氨酸拉链区对基因转录发挥广泛的调控作用［19］。显著的motif变化发生在成骨诱导的第7天，排序靠前的motif中出现转录因子和，这两种转录因子如何参与细胞功能的转变尚不清楚。

比各组间motif比较更有实际意义的是组间差异motif。测序结果显示，与未诱导的MSC相比，成骨诱导第3天和第5天显著变动的motif主要有、和bZIP家族中的等，提示其可能是调控成骨分化早期的重要转录因子。其中是为人们熟知的促进成骨分化的关键转录因子。基因属于TEA/ATTS结构域家族的成员，也被证实与成骨分化有关［20］。在成骨诱导第7天，组间差异motif变化与前5 d相比更为显著。和的motif数量下降，取而代之的则是基因家族和基因家族，表明主要参与成骨早期的细胞分化，同时也表明MSCs的成骨分化到了第7天后可能进入新的转录调控阶段。

GO富集分析的目的是对基因进行注释和分类。在MSC-0d、OB-3d、OB-5d和OB-7d各组细胞，GO富集分析的主要BP并无显著差异，大多是蛋白磷酸化、生物代谢等，说明以上这些功能可能是各类细胞基础生理活动所共有。而各组间差异peak邻近基因GO的富集分析显示不同的基因功能。在成骨诱导剂作用的第3天和第5天，上调的基因多与GTP酶活性有关。事实上，已有研究认为GTP酶活性是MSCs成骨定向分化的重要参与者［21］。此外，细胞外基质的形成、Wnt信号通路功能和骨骼发育相关功能也被激活，这与成骨方向的分化是密切相关的［22-24］。第7天，新增的基因功能转变为细胞黏附、细胞外基质组织、细胞与基质黏附等活动，表明细胞功能出现明显变化。此外，我们在peak邻近基因的信号通路富集分析和差异peak邻近基因信号通路富集分析也看到类似的变化。

综合分析，本项目通过追踪骨髓MSCs成骨分化早期染色质开放区域的动态变化，为研究成骨定向分化的调控机制及鉴定有效靶基因提供了新的切入点，后期工作将联合RNA-seq数据和ChIP-seq数据进行深入挖掘。

［1］ Johnell O， Kanis JA. An estimate of the worldwide prevalence and disability associated with osteoporotic fractures ［J］. Osteoporos Int， 2006， 17 （12）：1726-1733.

［2］ Muruganandan S， Roman AA， Sinal CJ. Adipocyte differentiation of bone marrow-derived mesenchymal stem cells： cross talk with the osteoblastogenic program［J］. Cell Mol Life Sci， 2009， 66（2）：236-253.

［3］ Botolin S， McCabe LR. Bone loss and increased bone adiposity in spontaneous and pharmacologically induced diabetic mice ［J］. Endocrinology， 2007， 148（1）：198-205.

［4］杲丽，胡成虎，金岩. 骨髓间充质干细胞在小鼠早衰性骨质疏松中的细胞生物学功能研究［J］. 中国病理生理杂志， 2013， 29（4）：707-712.

Gao L， Hu CH， Jin Y. Age-related changes of bone marrow mesenchymal stem cells in senile osteoporosis［J］.Chin J Pathophysiol， 2013， 29（4）：707-712.

［5］ Fiorentini E， Granchi D， Leonardi E， et al. Effects of osteogenic differentiation inducers onexpanded adult mesenchymal stromal cells ［J］. Int J Artif Organs， 2011， 34（10）：998-1011.

［6］ Arnsdorf EJ， Tummala P， Castillo AB， et al. The epigenetic mechanism of mechanically induced osteogenic differentiation［J］. J Biomech， 2010， 43（15）：2881-2886.

［7］ Scott-Browne JP， Lopez-Moyado IF， Trifari S， et al. Dynamic changes in chromatin accessibility occur in CD8+T cells responding to viral infection［J］. Immunity， 2016， 45（6）：1327-1340.

［8］ Zhao Y， Zheng D， Cvekl A. Profiling of chromatin accessibility and identification of general-regulatory mechanisms that control two ocular lens differentiation pathways［J］. Epigenetics Chromatin， 2019， 12（1）：27.

［9］ Buenrostro JD， Wu B， Chang HY， et al. ATAC-seq： a method for assaying chromatin accessibility genome-wide［J］. Curr Protoc Mol Biol， 2015， 109：21.29.1-21.29.9.

［10］ Langmead B， Salzberg SL. Fast gapped-read alignment with Bowtie 2［J］. Nat Methods， 2012， 9（4）：357-359.

［11］ Zhang Y， Liu T， Meyer CA， et al. Model-based analysis of ChIP-Seq （MACS）［J］. Genome Biol， 2008， 9（9）：R137.

［12］ Heinz S， Benner C， Spann N， et al. Simple combinations of lineage-determining transcription factors prime-regulatory elements required for macrophage and B cell identities［J］. Mol Cell， 2010， 38（4）：576-589.

［13］ Ramirez F， Dundar F， Diehl S， et al. deepTools： a flexible platform for exploring deep-sequencing data［J］. Nucleic Acids Res， 2014， 42（Web Server issue）：W187-W191.

［14］ Huang DW， Sherman BT， Lempicki RA. Systematic and integrative analysis of large gene lists using DAVID bioinformatics resources［J］. Nat Protoc， 2009， 4（1）：44-57.

［15］ Kanehisa M， Goto S， Sato Y， et al. KEGG for integration and interpretation of large-scale molecular data sets［J］. Nucleic Acids Res， 2012， 40（Database issue）：D109-D114.

［16］ Pittenger MF， Mackay AM， Beck SC， et al. Multilineage potential of adult human mesenchymal stem cells［J］. Science， 1999， 284（5411）：143-147.

［17］ Long HK， Prescott SL， Wysocka J. Ever-changing landscapes： transcriptional enhancers in development and evolution［J］. Cell， 2016， 167（5）：1170-1187.

［18］ Tsompana M， Buck MJ. Chromatin accessibility： a window into the genome［J］. Epigenetics Chromatin， 2014， 7（1）：33.

［19］ Hurst HC. Transcription factors 1： bZIP proteins［J］. Protein Profile， 1995， 2（2）：101-168.

［20］ Hakelien AM， Bryne JC， Harstad KG， et al. The regulatory landscape of osteogenic differentiation［J］. Stem Cells， 2014， 32（10）：2780-2793.

［21］ Strzelecka-Kiliszek A， Mebarek S， Roszkowska M， et al. Functions of Rho family of small GTPases and Rho-associated coiled-coil kinases in bone cells during differentiation and mineralization［J］. Biochim Biophys Acta Gen Subj， 2017， 1861（5 Pt A）：1009-1023.

［22］ Ward DF Jr， Salasznyk RM， Klees RF， et al. Mechanical strain enhances extracellular matrix-induced gene focusing and promotes osteogenic differentiation of human mesenchymal stem cells through an extracellular-related kinase-dependent pathway［J］. Stem Cells Dev， 2007， 16（3）：467-480.

［23］ Baksh D， Boland GM， Tuan RS. Cross-talk between Wnt signaling pathways in human mesenchymal stem cells leads to functional antagonism during osteogenic differentiation ［J］. J Cell Biochem， 2007， 101（5）：1109-1124.

［24］ Deng ZL， Sharff KA， Tang N， et al. Regulation of osteogenic differentiation during skeletal development［J］. Front Biosci， 2008， 13：2001-2021.

Dynamic changes of chromatin accessibility in early osteogenic differentiation of human mesenchymal stem cells

LIU Jian-yun， ZHANG Jie， MA Bai-cheng， WU Ping， XIONG Jian-jun

（，，332000，）

To observe the genome-wide dynamic profiling of chromatin accessibility in early osteogenic differentiation of human bone marrow mesenchymal stem cells （MSCs）.The 6th generation of primarily cultured human bone marrow MSCs were treated with osteogenic inducer for 0 d， 3 d， 5 d and 7 d， and then were harvested to receive the assay for transposase-accessible chromatin with high-throughput sequencing （ATAC-seq）. The numbers of open chromatin regions， the specific activated DNA motifs in open chromatin regions and the functional enrichment analysis of peak adjacent genes were analyzed by bioinformatics.At different time points in the MSCs at early stage of osteogenic differentiation， the number of open chromatin regions changed significantly， and functional peak distributions in each group were mainly enriched near the transcription start sites. Motif analysis revealed subtle alterations of activated DNA involved in transcriptional regulation between each group. Gene ontology analysis of peak's adjacent gene function showed that the biological processes of the newly opened chromatin regions were mainly concentrated on regulation of several GTPase activity， extracellular matrix organization， Wnt signaling pathway， and so on， while enrichment pathways in different genes included Rap1 signaling pathway， focal adhesion， adherens junction， and so on.During osteogenic differentiation of MSCs， the open chromatin regions change significantly at the early stage.

Mesenchymal stem cells； Osteogenic differentiation； ATAC-seq； Bioinformatics

R681； R363

A

10.3969/j.issn.1000-4718.2020.11.016

1000-4718（2020）11-2037-06

2020-03-14

2020-05-22

国家自然科学基金资助项目（No.81860165）

Tel: 0792-8577050； E-mail: jcyx_xiongjianjun@jju.edu.cn

（责任编辑：卢萍，罗森）