ChIP-seq技术分析突变IHH-GLI1信号传导通路对下游转录调控的影响

王 挺，沈 陆，魏慕筠，周 伟，李 沫，贺 林，秦胜营

(上海交通大学生命科学技术学院Bio-X研究院，中国 上海 200030)

A1型短指症(Brachydactyly type A1，BDA1)是一种罕见的遗传畸形，主要表现为手指和脚趾的中骨短小或缺失。研究发现，印度刺猬(Indian hedgehog，IHH)基因与BDA1密切相关，该基因上的杂合错义突变会导致软骨发育不全(从而导致BDA1发病)[1]。携带IHH基因E95K突变的BDA1模型小鼠具有与人类患者一样的个体矮小和典型的A1型短指/趾症表型，表现为第二至第四指中指节严重缩短，第五指中指节缺失[2]。既往研究认为E95K突变导致了IHH突变蛋白信号能力下降，使信号作用距离显著增加，削弱了其与受体Patched 1(PTCH1)和Huntingtin interacting protein 1(HIP1)的相互作用，从而影响软骨细胞的增殖和分化并最终导致骨骼生长减少[3]。IHH主要由软骨细胞表达，该蛋白与其受体PTCH1或其他因子相互作用后通过平滑蛋白(SMO)转导信号，最终作用于含锌指结构域的转录调节因子GLI(GLI 1～3)调节下游基因的表达。转录调节因子GLI1是Hedgehog(HH)信号通路的中央激活因子[4,5]。近年来，HH通路中GLI激活剂能够直接结合某些下游靶基因的重大发现对HH/GLI下游调节机制的理解提供了新见解[6]。最近的研究表明，IHH/GLI信号传导到GLI的梯度活动可以协调发育的模式，GLI转录因子接收IHH信号的强弱则决定了不同的转录输出[7]。基于基因芯片技术的研究表明BDA1相关的突变IHH信号通路会改变GLI介导的转录调控模式，进而导致了BDA1的发生[8]。然而目前尚未有研究从全基因组范围内精细分析突变IHH-GLI1信号下游转录调控，而这对于了解BDA1病理机制和靶向治疗是必不可少的。

近年来，染色质免疫共沉淀测序技术(Chromatin immunoprecipitation followed by sequencing， ChIP-seq)在转录因子调控的靶基因鉴定中的应用越来越多[9]。相对于基因芯片技术，ChIP-seq技术具有高分辨率和高精准度的优势，能够捕捉到未知的序列，发现和寻找到新的转录因子结合位点。研究人员可将获得的与目的蛋白结合的数百万序列标签准确地定位到全基因组上，从而全面、准确地鉴定与蛋白质结合的基因片段信息。该技术现已被广泛应用于疾病的诊断以及疾病治疗靶点的探究中[10]。目前对IHH信号介导的BDA1发病机制和有效治疗靶点仍不清楚，因此，本研究基于ChIP-seq技术从基因组范围内全面地对突变IHH信号下GLI1介导的基因转录调控进行分析。本研究向小鼠胚胎来源的间充质干细胞系C3H10T1/2中加入人IHH-N野生型重组蛋白(WT)和突变型蛋白(MT, p.E95K)以激活IHH信号通路。IHH蛋白作用2天后，收集并裂解细胞，使用GLI1抗体进行染色质免疫沉淀(ChIP)实验并将DNA产物用于ChIP-seq分析。通过生物信息学分析对GLI1蛋白结合的靶基因进行注释和定位，进一步挖掘两组间不同的GLI1转录因子结合基因，并对这些基因进行功能注释和KEGG通路分析；同时将ChIP-seq分析结果在基因芯片数据中进行验证，以构建IHH/GLI1下游通路网络。

1 材料和方法

1.1 材料

1.1.1 菌株和细胞 重组人IHH蛋白的大肠杆菌BL21(DE3)表达菌株，包括表达野生型蛋白(氨基酸残基28-202)和突变蛋白(E95K)两种菌株，构建过程见文献[3]，由本研究院保存。小鼠胚胎来源的间充质干细胞株(C3H10T1/2)购自ATCC。

1.1.2 主要试剂和仪器 MEM/EBSS培养基、 胎牛血清(FBS)、蛋白酶抑制剂和青霉素/链霉素双抗购自Gibco公司； GenomePlex©全基因组扩增试剂盒购自Sigma-Aldrich公司； EZ-ChIP 试剂盒、0.22 μm 滤膜、GLI1多克隆抗体均购自Millipore公司；IHH 抗体(Santa Cruz)； GST 标签纯化装置(GE-Healthcare)；蛋白 G 琼脂糖(Pierce)、RNA酶和蛋白K酶、TRIZOLTM试剂均购自Invitrogen；PCR 纯化试剂盒 (Qiagen)；NEBNext DNA 超快速文库制备试剂盒 (NEB)； 转录子高保真cDNA合成试剂盒、SYBR Green (Rox) 试剂盒购自Roche公司； Illumina测序试剂和用品均购自诺唯赞生物科技有限公司； ViiATM7实时荧光PCR 仪(Applied Biosystems)；安捷伦2000生物分析仪以及芯片分析所有试剂、仪器均来自安捷伦公司。

1.2 方法

1.2.1 IHH蛋白表达和纯化 将重组人IHH-N蛋白的大肠杆菌BL21(DE3)表达菌株(包括野生型蛋白氨基酸残基28-202；突变蛋白E95K)涂布平板后，选择单克隆菌落，并在5 mL含氨苄青霉素的LB液体培养基中于37 ℃孵育过夜。然后将培养基转移到含有氨苄青霉素的400 mL 2×YT液体培养基中，并在37 ℃下孵育4 h。当诱导菌株的OD值达到0.6时，添加0.2 mmol·L-1IPTG在32 ℃条件下诱导蛋白表达5 h。在4 ℃下以8 000×g离心15 min，加入含有蛋白酶抑制剂的200 mmol·L-1氯化钠和20 mmol·L-1Tris-盐酸(pH 7.5)溶液重悬细菌超声破碎提取蛋白。在4 ℃下以25 000×g离心30 min后，收集上清液。 然后根据GST纯化装置说明书纯化蛋白质，使用凝血酶裂解GST标签。纯化的IHH-N蛋白通过Western blot进行验证，并用0.22 μm滤膜进行过滤。

1.2.2 C3H10T1/2细胞中IHH信号通路的激活 C3H10T1/2细胞的IHH激活按照文献[5]的方法进行诱导。将C3H10T1/2细胞培养在添加MEM/EBSS完全培养基(包含MEM/EBSS基础培养基、10%胎牛血清(FBS)和2%青霉素/链霉素)的12孔板中，在37 ℃(5%CO2)的恒温培养箱中培养。当细胞生长到每孔5×104时，将WT或MT IHH-N(E95K)添加到生长培养基中(最终IHH-N蛋白为1×=750 nmol·L-1)。 在37 ℃(5%CO2)下进一步孵育48 h后，收获细胞用于后续的进一步分析。 通过实时定量PCR检测IHH信号通路标志物分子Gli1，Ptch1和Ihh的表达。所有诱导测定均一式三份进行。

1.2.3 染色质免疫共沉淀(ChIP)样品的准备 使用EZ-ChIP试剂盒进行染色质免疫沉淀。分离剪切的交联染色质后，将每个样品(1×105个细胞)与5 μg多克隆GLI1抗体(或匹配的IgG对照)孵育过夜。然后将60 μL Protein G磁珠添加到每个样品中，并在4 ℃旋转孵育1 h。孵育后，通过短暂离心将Protein G磁珠沉淀，并根据说明书进行洗涤。将沉淀物与100 μL稀释液(10 μL 20%SDS，20 μL 1 mol·L-1NaHC03和70 μL无菌蒸馏水)在室温下孵育15 min，收集上清液。向每个样品收集液(总计200 μL)中加入8 μL的5 mol·L-1NaCl，并在65 ℃下温育过夜，以逆转DNA-蛋白质交联。第二天，所有样本依次用RNA酶A和蛋白酶K处理。酶处理后的样本中加入酚氯仿抽提DNA，使用乙醇沉淀法纯化 DNA, 添加TE 试剂溶解DNA 后进行测序文库的构建。

1.2.4 ChIP文库构建、质检及测序 取1 μg input 对照DNA和 ChIP实验DNA, 按照 NEBNext DNA 超快速文库制备试剂盒的要求规范制备测序文库。构建好的文库使用 Bioanalyzer检测文库DNA片段的完整性及插入片段大小，用定量PCR对文库DNA拷贝数进行精确定量。随后将构建合格的文库按照有效浓度及目标下机数据量的需求进行混合，使用Illumina高通量测序平台 (HiSeq) 进行双端测序。

1.2.5 实时荧光定量PCR验证 使用转录子第一链cDNA合成试剂盒将 RNA合成cDNA。依据SYBR Green Master Rox试剂盒说明书，加入特殊引物进行PCR反应。PCR反应条件： 94 ℃变性10 min，再94 ℃变性15 s，60 ℃退火40 s，总共40个循环。将反应好的PCR体系放在ViiATM7实时定量PCR系统上进行检测。每个PCR检测重复3次，使用Gapdh和Hprt作为内参基因对检测基因的表达水平进行归一化。根据2-ΔΔCT方法计算相对表达值，当2-ΔΔCT>1或<1认为检测基因上调或下调。所有验证引物的序列见表1。

表1 RT-PCR验证引物

1.2.6 ChIP-seq数据生物信息学分析 本研究使用自动化分析程序ChiLin对ChIP-seq数据进行质量控制和数据分析[11]。FastQC检查原始序列质量和GC含量，默认比对工具BWA将ChIP和对照组的FASTQ文件比对到参考基因组mm9上。完成基因组比对后，使用默认峰调用程序MACS2 (https://github.com/taoliu/MACS/)来执行窄峰(用于点源绑定)或宽峰的富集。进一步通过qvalue计算报告峰的错误发现率以及每个峰的褶皱富集程度。Cistrome平台中输入ChiLin 分析完成的bed 和wig格式的最终峰文件，通过对顺式调控元件结合序列峰在基因组上的特征进行注释，以找到下游离它最近的基因和最接近的转录起始位点。UCSC 来源的参考基因注释文件将每个峰标注到从峰值到转录起始位点(transcription start site, TSS) 25 kb范围内的任何基因上以探索GLI1结合峰与潜在靶基因之间的相关性。

1.3 GO和 KEGG 分析以及蛋白互作分析工具

采用R扩展包Cluster Profiler对GLI1结合的靶基因进行GO(gene ontology，基因本体论)注释和KEGG(Kyoto Encyclopedia of Genes and Genomes)通路富集分析，筛选标准为P<0.05。STRING(http://string-db.org/)是一个预测蛋白质间功能相关性的在线数据库[12], 本研究利用STRING提取靶基因之间的相互作用关系(score>0.4)。

1.4 数据统计

统计结果的所有统计量均以3个独立实验的平均值M±SEM表示，采用t检验比较IHH正常组与突变组之间的差异，P<0.05为显著性差异。

2 结果

2.1 IHH-N突变蛋白抑制的C3H10T1/2细胞IHH信号通路的激活

将C3H10T1/2细胞与750 nmol·L-1IHH-N蛋白孵育48 h后，使用实时定量PCR在两组中检测激活的IHH信号传导的靶标Gli1，Ptch1和Ihh的表达水平,发现诱导后Ihh和Gli1上调,WT组的Gli1表达高于MT组(3.63vs2.07，P=0.051)(图1)。这些结果表明E95K突变削弱了IHH信号的传导作用。有趣的是，笔者观察到，属于Hedgehog(HH)信号传导的12个跨膜结构域受体Ptch1在WT组中被上调，而在MT组中被下调(2.27vs0.67，P=0.093)。IHH和受体PTCH1之间的结合有利于激活IHH/GLI1介导的下游靶基因的表达(图1)。

图1 激活IHH后，Gli1，Ihh和Ptch1的相对mRNA表达水平 Fig. 1 Relative mRNA expression level of Gli1, Ihh and Ptch1 after IHH activation

2.2 IHH突变蛋白作用下GLI1蛋白-DNA结合位点的全基因组表征

笔者使用Cistrome平台对预测到的GLI1结合区域在小鼠基因组上进行了定位与注释。结果表明，在染色体位置分布上，WT组的ChIP区域在所有染色体上均匀分布，而MT组在4, 9, 11, 13染色体上的富集比例明显升高(图2A)。为了获得GLI1结合序列在基因组上的分布情况，笔者分别统计了WT, MT以及对照组的ChIP位点在基因组不同区域的比例。相对于CT组，WT和MT组的ChIP位点在基因组上启动子区的分布比例明显增加。MT组在基因组上明显富集于启动子区(≤1 000 bp)，而在基因组下游(≤1 000 bp)区域的富集率相对较小(图2B)。这表明MT组的GLI1作为转录调节因子，主要在基因近端调控区发挥作用。随后笔者对预测到的GLI1结合位点富集到基因组上的峰数目进行统计，鉴定出WT组201个峰、MT组92个峰和CT组87个峰(峰选择阈值=p值<10e-5)。在WT和MT组之间，GLI1结合位点重叠区域较少(图3A)，这表明突变的IHH信号改变了GLI1转录因子的下游调节机制。

2.3 IHH突变蛋白作用下GLI1转录因子结合靶基因的鉴定

本研究使用UCSC小鼠参考基因组注释了GLI1结合位点25 kbp之内的序列以表征潜在的GLI1调控靶基因，结果显示在GLI1结合位点鉴定到WT组332个基因、MT组151个基因和CT组101个基因(图3B)。WT组中GLI1结合的基因似乎比MT组更多，并且两组中鉴定到共同的GLI1结合基因有18个(图3B)。稳定表达GLI1的细胞基因表达谱揭示了Slx1b是GLI1结合元件的直接下游靶标[13]，这也表明ChIP-seq技术分析结果具有较高的可靠性。

进一步对在WT 和MT组中预测到的特有的GLI1结合基因进行GO注释，结果显示，WT组的GLI1结合基因主要细胞成分为细胞核和细胞内膜结合细胞器(图4A)。GLI1结合靶基因参与了多种生物学过程，包括DNA模板转录的调控、RNA生物合成的调控、细胞大分子生物合成调控等，其中大部分为调控相关基因，表明GLI1可能参与细胞分化过程中的物质合成调控过程(图4B)。而MT组中GLI1结合靶基因无明显富集的分子功能和生物学过程。进一步对两组中GLI1结合的靶基因进行KEGG通路富集分析，结果表明，与WT组相比，MT组特有的GLI1结合靶基因富集到的信号通路明显减少(P<0.05)(图5)。两组共有Rap1和ErbB信号通路,MT组主要缺失Wnt信号通路、刺猬因子(HH)信号通路、胰岛素样生长因子(IGF)信号通路、mTOR信号通路、TNF信号通路以及Th17信号通路等，其中刺猬因子(HH)信号通路、胰岛素样生长因子(IGF)信号通路和Wnt信号通路被报道是参与调控生长板软骨的生长和发育的关键信号通路。

(A)ChIP区域在染色体上的分布统计；(B) ChIP区域在全基因组位置上的分布统计，启动子区(≤1 000,≤2 000,≤3 000 bp)和下游调控区(≤1 000,≤2 000,≤3 000 bp)。

(A) GLI1结合位点的维恩图；(B)预测的GLI1结合基因的维恩图。图3 ChIP-seq分析结果Fig. 3 Results from ChIP-seq analysis

(A)细胞成分；(B)生物学进程。图中每个图形为一个GO term，其中纵坐标代表GO名称，横坐标代表富集率，气泡的大小代表富集到的基因数目，气泡颜色越深表示该GO功能的富集越明显。FDR代表矫正后的P值。

(A) WT组；(B) MT组。每条柱的长度代表每条通路富集到的基因数目，颜色越深表示该通路的富集越显著。

分别对两组中鉴定到的特有GLI1结合靶基因进行STRING分析，以提取靶基因之间的相互作用关系。结果显示，相对于WT组，MT组的GLI1结合靶基因的蛋白质互作网络节点的数量明显减少，相互作用明显减弱(图6)。进一步对MT和WT组中预测到的GLI1结合的基因进行比较，发现两组间存在465个不同的基因，对这些基因进行通路富集，列出前10个显著富集的通路及通路上的候选基因，见表2。以上结果表明，E95K突变信号极大地改变了GLI1介导的下游转录调控，从而削弱了软骨细胞增殖和分化中重要参与因子的应答。

表2 前10个明显富集的KEGG信号通路

(A)WT组特有的靶基因相互作用关系；(B)MT组特有的靶基因相互作用关系。线条粗细代表蛋白与蛋白之间的关联程度，线条越粗，代表蛋白之间关联度越大。图6 GLI1结合的蛋白质相互作用网络Fig. 6 Interaction networks of GLI1 binding proteins

2.4 ChIP-seq分析鉴定到的GLI1结合基因在基因芯片上的验证

将ChIP-seq预测的GLI1结合靶基因集与本实验室前期的ChIP-chip数据集[8]进行联合分析, 找到19个共同的候选基因。对这些候选基因进行生物学功能、通路及编码蛋白互作分析。GO分析结果显示，这些基因主要参与细胞与细胞间信号转导等生物过程；KEGG通路富集分析表明这些基因参与cAMP信号通路和刺猬因子(Hedgehog, HH)信号通路(P<0.05)；蛋白互作网络展示了功能基因编码蛋白的互作关系，结果显示刺猬信号通路下游GLI1调控的近端信号靶基因如Nkx2-5，Eya1和Foxa1之间存在较强的互作关系，可能是BDA1发生发展的关键基因(图7)。同时，笔者也筛选到了感兴趣的分子如Clec3a,Mtl5,Foxa1等，为进一步探索BDA1的发病机制和治疗靶点提供理论基础。

注：两个节点之间的线段代表蛋白质之间的相互作用。 图7 GLI1结合的差异靶基因的蛋白互作网络Fig. 7 The protein-protein interaction network built with the differentially expressed GLI1-binding target gene

3 讨论与结论

HH/GLI信号通路通过协调软骨细胞的增殖和分化来调节软骨内骨形成的多个方面并控制骨骼的生长。在软骨发育过程中，GLI介导的IHH信号在脊椎动物中调控下游基因的表达[14]。然而目前仍缺乏从全基因组范围内的分析，研究突变IHH-GLI1信号是如何影响下游转录调控模式，进而导致BDA1的发生。本研究中，笔者首先使用野生型和突变型IHH蛋白作用于小鼠间充质干细胞C3H10T1/2诱导IHH通路激活，通过实时定量PCR检测和比较两组间IHH通路相关因子转录水平变化。结果表明，E95K突变会削弱IHH信号在转录调控中的作用，这与前人关于E95K突变体IHH会在软骨发育异常的小鼠模型中损害与其结合的下游靶基因的表达的研究一致[8]。随后，使用ChIP-seq技术采用生物信息学方法鉴定和比较IHH蛋白正常组和突变组中与GLI1转录因子结合的靶基因，并采用GO注释和KEGG通路对两组中特有的靶基因进行了生物学功能表征。在ChIP-seq分析结果中，笔者发现 E95K突变削弱了IHH信号的传导作用，使富集到的GLI1结合的靶基因明显减少。缺少的信号通路主要有刺猬因子(HH)信号通路[15]、胰岛素信号通路[16]和Wnt信号通路[17]等参与调控生长板软骨的生长和发育的关键信号通路。本研究进一步解释了由于IHH p.E95K突变改变的IHH/GLI1下游调节机制，为进一步探索BDA1的发病机制和治疗靶点提供理论基础。

基于对前10个显著富集的通路的评估，笔者发现Prkaca基因参与的信号通路最多，发挥着重要的功能。Prkaca基因编码蛋白激酶A的一个催化亚基，进而参与cAMP的蛋白磷酸化，在许多细胞过程(包括分化、增殖和凋亡)中发挥着重要作用。另一个基因Mapk3也参与多个信号通路，其编码的蛋白质是细胞外信号调节激酶(MAP)家族的成员, 在信号级联反应中起作用，响应各种细胞外信号调节各种细胞过程，例如增殖、分化和细胞周期进程。这也表明E95K突变信号改变了GLI1介导的下游转录调控表达模式，其信号下游应答因子主要参与与软骨形成相关的细胞分化、增殖和凋亡等生物学过程。此外，KEGG通路分析提示GLI1结合的候选靶基因在突变IHH-GLI1信号通路的响应应答中发挥着重要作用，如参与软骨发育关键信号通路Wnt 信号通路中的Lgr6，Axin2，Map3k7，Dkk1，Prkaca，Apc及Ctbp2；刺猬因子(HH)信号通路中的Spop及Ptch2；IL-17信号通路中的Hsp90aa1及Fosb；胰岛素信号通路中的Sh2b2，Calm2，Ppp1r3a及Mapk3等。

高通量测序和基因芯片等现代生物技术的发展以及生物信息学应用为我们从分子水平揭示E95K突变改变的IHH/GLI1下游调控机制以解释BDA1的发病机理提供了很好的分析手段。这两种技术的结合提高了转录调节靶基因的精准鉴定。笔者将ChIP-seq预测的GLI1结合靶基因集与本实验室前期的基因芯片数据集进行联合分析, 找到了19个共同的候选基因。GO和KEGG分析显示这些基因主要参与细胞与细胞间信号转导等生物过程以及参与cAMP和刺猬因子(HH)等重要信号通路(P<0.05)。

蛋白互作网络展示了功能基因编码蛋白的互作关系，笔者发现刺猬因子(HH)信号通路下游GLI1调控的近端信号靶基因Nkx2-5，Eya1和Foxa1之间存在较强的互作关系。Nkx2-5基因被报道是骨髓间充质干细胞(MMSC)向骨髓成心肌细胞分化的早期形成关键基因[18]。有研究表明MMSC是成体多能干细胞，在不同诱导条件下可以分化为软骨、骨、心肌等多种组织的细胞。因此，笔者猜想本研究中E95K突变信号激活了C3H10T1/2细胞向其他细胞分化的参与因子，从而削弱了IHH信号在转录调控中的作用，导致BDA1的发生。然而，具体的作用机制还需要在未来的工作中进一步验证。另外，本研究筛选到的Clec3a基因是成骨作用和骨化作用因子，其过表达会在创伤性颞下颌关节强直模型动物出现过度成骨现象，从而导致颞下颌关节强直疾病的发生[19]；Mtl5是参与骨关节炎和骨质疏松症的候选基因[20]；Foxa1参与骨癌的发展过程[21]。

综上所述，本研究明确了IHH蛋白突变影响IHH通路蛋白转录，并通过ChIP-seq对GLI1转录因子结合的靶基因位点进行分析，结合基因芯片数据的结果，分析和鉴定了IHH野生型蛋白和E95K突变型蛋白导致的GLI1-DNA结合靶基因的差异。GO分析和KEGG分析提示这些基因具有骨化作用的生物学功能并参与了刺猬信号骨发育关键信号通路。此研究结果对体内E95K突变如何改变GLI1介导的基因调控并最终损害BDA1模型中的IHH信号传导能力做出了新的解释，这将为理解IHH信号在BDA1的软骨形成作用以及探索BDA1的发病机制和治疗靶点提供有用的信息。