- Research article
- Open Access
p.E95K mutation in Indian hedgehog causing brachydactyly type A1 impairs IHH/Gli1 downstream transcriptional regulation
BMC Genetics volume 20, Article number: 10 (2019)
Brachydactyly type A1 (BDA1, OMIM 112500) is a rare inherited malformation characterized primarily by shortness or absence of middle bones of fingers and toes. It is the first recorded disorder of the autosomal dominant Mendelian trait. Indian hedgehog (IHH) gene is closely associated with BDA1, which was firstly mapped and identified in Chinese families in 2000. Previous studies have demonstrated that BDA1-related mutant IHH proteins affected interactions with its receptors and impaired IHH signaling. However, how the altered signaling pathway affects downstream transcriptional regulation remains unclear.
Based on the mouse C3H10T1/2 cell model for IHH signaling activation, two recombinant human IHH-N proteins, including a wild type protein (WT, amino acid residues 28–202) and a mutant protein (MT, p.E95k), were analyzed. We identified 347, 47 and 4 Gli1 binding sites in the corresponding WT, MT and control group by chromatin immunoprecipitation and the overlapping of these three sets was poor. The putative cis regulated genes in WT group were enriched in sensory perception and G-protein coupled receptor-signaling pathway. On the other hand, putative cis regulated genes were enriched in Runx2-related pathways in MT group. Differentially expressed genes in WT and MT groups indicated that the alteration of mutant IHH signaling involved cell-cell signaling and cellular migration. Cellular assay of migration and proliferation validated that the mutant IHH signaling impaired these two cellular functions.
In this study, we performed integrated genome-wide analyses to characterize differences of IHH/Gli1 downstream regulation between wild type IHH signaling and the E95K mutant signaling. Based on the cell model, our results demonstrated that the E95K mutant signaling altered Gli1-DNA binding pattern, impaired downstream gene expressions, and leaded to weakened cellular proliferation and migration. This study may help to deepen the understanding of pathogenesis of BDA1 and the role of IHH signaling in chondrogenesis.
The Hedgehog (HH) signaling pathway is a key signaling pathway that regulates a wide range of developmental processes in vertebrate and invertebrate organisms [1,2,3]. HH proteins in mammals consist of three homologues: Indian hedgehog (IHH), Sonic hedgehog (SHH) and Desert hedgehog (DHH). In mammal developing cartilage, IHH is primarily expressed by chondrocytes. IHH protein interacts with its receptor Patched 1 (PTCH1) or other partners, followed by the transduction of the signaling through Smoothened (SMO) protein. This process ultimately regulates expression of downstream genes by Gli zinc finger domain-containing transcription regulators (Gli1 to − 3) [4, 5]. Transcription regulator Gli1 is a central activator of the HH pathway, while Gli2 and Gli3 serve as activators or as repressors depending on HH signaling level [6,7,8]. HH/Gli signaling pathway regulates multiple aspects of endochondral bone formation and controls the growth of skeleton by coordinating chondrocyte proliferation and differentiation . Studies of chondrogenesis suggested that there existed a complicated network of various signal interactions  involving IHH signaling. Gli-mediated regulation of IHH signaling in vertebrates controls downstream gene expressions during the cartilage development . The significant discoveries of the direct targets of Gli-activator in HH signaling made in last few years [12,13,14,15,16] have provided new insights to our understanding of general mechanism of HH/Gli downstream regulation.
Heterozygous missense mutations in IHH result in brachydactyly type A1 (BDA1; OMIM 112500), a condition characterized by the shortening of digits due to hypoplasia/aplasia of the middle phalanx [17, 18]. In 2009 and 2011, two studies have revealed that the BDA1-related mutant IHH impaired its interaction with PTCH1 and HIP1 and lowered the responses close to the source but increased those at a distance. The final enhanced IHH signaling range might cause reduced bone growth [5, 19]. Recent work has indicated that HH signaling transduced into a gradient of Gli activity that orchestrates patterning of development, while relatively small differences in the ration of Gli forms (Gli activators and repressors) dictated distinct transcriptional outputs . However, how the BDA1-related mutant IHH signaling changes Gli-mediated transcriptional regulation remains unclear.
In this study, we aimed to initiate an effort to characterize the alteration of Gli1-mediated transcriptional regulation caused by E95K mutant signaling. First, C3H10T1/2 cells were induced with wild type (WT) or E95K mutant (MT) IHH protein to activate IHH signaling. After 2 days’ incubation, chromatin immunoprecipitation (ChIP) was conducted and then the DNA products were used to interrogate ChIP-chip analysis for identifying and mapping gli1-binding sites. Then we performed microarray-based gene expression analysis on RNA derived from the induced cells. All the data was integrated by bioinformatics analyses to construct IHH/Gli1 downstream pathway networks. In sum, our results demonstrated significant differences of cartilage-specific Gli1-mediated transcriptional regulation between WT group and MT group. This discovery will help us to unveil the underlying intracellular mechanisms that modulate pathogenesis with BDA1 and to promote the cartilage defect therapy.
IHH protein expression and purification
The pGEX-2 T-based recombinant human IHH-N proteins, including wild type protein (amino acid residues 28–202) and a mutant protein (E95k), were expressed by the BL21 (DE3) strain of Escherichia coli (E. coli). We obtained the two transformed strains from Bio-X Institutes, Shanghai Jiao Tong University, which had been published earlier . After coated plates, monoclonal colony was selected and incubated in 5 ml LB liquid medium with ampicillin at 37 °C overnight. Then the medium was transferred to 400 ml 2 × YT liquid medium with ampicillin and incubated at 37 °C for 4 h. OD600 reached 0.6 when the strains were induced with 0.2 mM IPTG at 32 °C for 5 h.. Pellets were collected by centrifugation at 8000×g for 10 min at 4 °C and re-suspended in 200 mM NaCl and 20 mM Tris-HCl, pH 7.5 with cocktail protein inhibitors (Sigma-Aldrich), followed by sonication. After centrifugation at 25000×g for 30 min at 4 °C, the supernatant was collected. Then the proteins were purified according to the Bulk and RediPack GST Purification Modules (GE-Healthcare), followed by cleavage of the GST tag using thrombin. The purified IhhN proteins were verified by western blotting with IHH Antibody (SC-1782, Santa Cruz) and filtered by 0.22 μm filter (Millipore).
IHH activation of C3H10T1/2 cell
The induction was performed following the method that had been described in previous study . C3H10T1/2 cells (obtained from ATCC) was cultured in 12-well plates in a growth medium containing MEM/EBSS, 10% fetal bovine serum (FBS), and 2% penicillin/streptomycin. WT or MT IhhN (E95K) was added into the growth mediums (1x = 750 nM IhhN protein finally) when cells grow to 5 × 104 per well. After 48 h of further incubation at 37 °C (5% CO2), the cells were harvested for further analysis. Expression of Gli1, Ptch1 and Ihh, as markers of IHH signaling pathway, were detected by real time PCR. All the induction assays were performed in triplicate.
We used EZ-ChIP Kit (Millipore) to perform chromatin immunoprecipitation. After isolating sheared crosslinked chromatin, each sample (1 × 105cells) was incubated overnight with 5 μg polyclonal Gli1 antibody (AB3444, Millipore) or matched IgG controls. Then 60 μL Protein G Agarose was added to each sample and incubated for 1 h at 4 °C with rotation. After incubation, Protein G Agarose was pelleted by brief centrifugation and washed according to the instruction manual. The pellets were incubated with 100 μL elution buffer (10 μL 20% SDS, 20 μL 1 M NaHCO3 and 170 μL sterile, distilled water) at room temperature for 15 min twice and the supernatant was collected. To each sample collection (totally 200 μL), 8 μL 5 M NaCl was added and incubated at 65 °C overnight to reverse the DNA-protein crosslink. On the next day, the samples were sequentially treated with RNAse A and Proteinase K. DNAs were purified using PCR Purification Kit (Qiagen) and amplified using GenomePlex® Whole Genome Amplification Kit (Sigma-Aldrich).
For each ChIP, 2 μg DNA was labeled with Cy5-dUTP for amplified ChIP sample and Cy3-dUTP for input sample and hybridized on MM9 ChIP 3x720K RefSeq Promoter microarray (Roche Nimblegen) according to the manufacturer’s recommendation. Three independent biological sample groups were employed in the analysis and one of the three groups was chosen to conduct further analyses based on its IHH activation level. Arrays were scanned using a Nimblegen MS 200 Microarrary Scanner at 2 μm resolution. Array signals were extracted using NimbleScan Software version 2.6 (Roche Nimblegen) according to the manufacturer’s instruction.
Microarray-based gene expression analysis
We used TRIZOL™ solution (Invitrogen) to extract the total RNA from cell samples according to the manufacturer’s instruction. RNAs were examined by Agilent 2100 Bioanalyzer to assure its RIN to be more than 7. MicroRNA microarray analyses were performed using 200 ng total RNA on Mouse GE 4x44K v2 Microarray (Agilent Technologies). The amplification, labeling, hybridization and scanning of samples were accomplished by Agilent Technologies (Shanghai, China) using standard protocols. The data were extracted by Feature Extraction version 12.5 (Agilent Technologies) and analyzed by Genespring GX version 12.6 (Agilent Technologies).
Validation of Results from Microarray-based Analyses by Quantitative Real-time PCR
All qPCR data generated were compliant with current Minimum Information for Publication of Quantitative Real-Time PCR Experiment (MIQE) standards according to Bustin et al. .
RNA was reverse transcribed to cDNA using the transcriptor high-fidelity cDNA synthesis kit (Roche), with random hexamers according to the manufacturer’s instructions. We used qPCR with SYBR Green to validate results of microarrays. Primers for ChIP-chip and gene expression microarray were designed by Beacon Designer version 7.9 (PREMIER Biosoft). The PCR reaction was performed using the FastStart Universal SYBR Green Master (Rox) Kit (Roche Diagnostics) on ViiA™ 7 Real-Time PCR System (Applied Biosystems). Each sample in triplicates was assessed in denaturation at 94 °C for 10 min, following by 40 cycles of denaturation at 94 °C for 15 s and annealing at 60 °C for 40 s. Expression levels were assessed in triplicate, normalized to Gapdh and Hprt control levels. All the sequences of the validation primers are shown in Additional file 1.
Microarray data process
All the microarray data generated were compliant with current Minimum Information About a Microarray Experiment (MIME) standards according to Brazma et al. .
ChIP-chip data process
Raw data of microarrays were generated in pair format using NimbleScan software version 2.6 (Roche Nimblegen). Peak calling was accomplished by MA2C peak calling algorithm (version 1.0.0, P = 1e-5, Robust method, C = 2) . The raw data had been submitted to NCBI Gene Expression Omnibus (GSE74022) . Resulting peak files in bed and wig formats were used as input for Cis-regulatory Element Annotation System  to annotate peaks with genome features within Cistrome  galaxy environment.
To discover de novo motifs of Gli1, MEME-CHIP  analysis were conducted on around 500 bp sequences on surrounding peak centers for peak sets. Cluster analysis of peaks were performed following functional enrichment analysis with Database for Annotation, Visualization and Integrated Discovery (DAVID) [28, 29] and Gene Ontology Analysis [30,31,32]. We used Reactome platform to conduct pathway analysis of the predicted genes .
To identify correlation between Gli1-binding peaks and potential target genes, we annotated each peak to any gene within 25 kb from peaks to transcription start site (TSS) region according to the UCSC browser within the Galaxy environment.
Gene expression microarray data process
The differentially expressed genes (DEGs) among three groups were identified and clustered using Genespring GX version 12.6 (Agilent Technologies). DEGs were filtered by a criteria of fold change> 2.0. The raw data could be obtained from GEO (GSE74021).
Integrated analysis of microarray data
We used Ingenuity Pathway Analysis (IPA, Qiagen Corp.) to integrate transcriptional DNA-binding sites and mRNA expression for prediction of interactions between signaling pathways. We used Galaxy (https://usegalaxy.org/) for sequence manipulating, correlation computing and genomic interval operation [34,35,36]. We used Binding and Expression Target Analysis to combine ChIP-chip with differential gene expression data to infer direct target genes . We used VennDiagram in R to draw venn diagrams .
Cell migration and proliferation test
For cell migration test, C3H10T1/2 cells were seed into 6-well tissue culture plates and grew in MEM/EBSS medium with 10% FBS until they reach ~ 70–80% confluence as a monolayer. We scratched the monolayer with a new 10 μl pipette tip across the well. After scratching, the cells were washed twice with medium and replenished with fresh medium. After 48 h of incubation, the cells were washed twice with 1x PBS and fixed with 3.7% paraformaldehyde. The gap distances were evaluated quantitatively by ImageJ software. Experiment in each group was run in triplicate.
For cell proliferation test, we used MTT Cell Proliferation and Cytotoxicity Assay Kit (Beyotime, China). C3H10T1/2 cells were seed into 96-well tissue culture plate and each well contained ~ 1000 cells within 100 μl MEM/EBSS medium with 10% FBS. After 48 h of culturing with or without IHH-induction, each well was added with 10 μl MTT and the incubation continued for 4 h. Then 100 μl formanzan solution was added into each well. After 4 h of incubation, the absorbance was calculated at 570 nm. Experiment in each group was repeated in triplicate.
Activation of IHH signaling pathway in C3H10T1/2 cells induced by IhhN proteins
According to previous BDA1 cell model , after incubation of C3H10T1/2 cells with 750-nM IhhN protein for 48 h, targets of activated IHH signaling, such as Gli1, Ptch1 and Ihh, were examined using quantitative real time PCR in two groups (Fig. 1). Ihh and Gli1 were up-regulated after induction. The expression of Gli1 was higher in WT group than in MT group (3.63 vs 2.07, p = 0.051). These results suggested that E95K mutant IHH weakened the activation of IHH signaling. Interestingly, we observed that Ptch1, a 12-transmembrane domain receptor belonging to the Hedgehog (HH) signaling, was up-regulated in WT group but down-regulated in MT group (2.27 vs 0.67, p = 0.093). The binding between IHH and PTCH1 contributes to the blocking of inhibitory effects of PTHC1 on the activity of SMO, thus activating the IHH signaling downstream targets .
Genome-wide characterization of Gli1-binding sites reveals the differences of downstream regulation in response to IHH signaling between two groups
Based on data from ChIP-chip experiments, we conducted the peak calling using MA2C algorithm  and identified 347 peaks in WT group, 47 peaks in MT group and 4 peaks in control group (peak selection threshold = p-value <10e-5) (Additional file 2). Among WT, MT and control group there was few overlapping of Gli1 binding sites (Fig. 2a), which indicated that the alteration of IHH signaling had a great impact on Gli1-mediated downstream regulation. We performed de novo motifs searching on the sequences derived from ChIP-chip peaks. However, no significant Gli1 binding consensus sequences were found in any dataset. To verify the ChIP-chip results, a subset of Gli1-binding sequences was validated by real time PCR (Additional file 1).
Characterization of cis-regulation targets mediated by Gli1 in response to IHH signaling
Based on the NCBI RefSeq genes database, we annotated sequences within 25 k bp of Gli1 binding sites to characterize potential Gli1-regulated targets and identified 795 genes in WT group, 109 genes in MT group and 11 genes in control group using the methods described in earlier study  (Additional file 3). Gli1 appeared to regulate more targets in WT group than in MT group and there existed a very small number of overlapping candidate genes in two groups (Fig. 2b). The result of gene enrichment analysis showed that there were significant enrichments of cell-to-cell signaling and interaction related biological processes among candidate genes in WT group, while there were no significant enrichments found in MT and control group (Table 1).
Then we performed microarray-based analysis for mRNA expression profile. Differential genes expression (DEGs, fold change > 2 compared to control group) was analyzed in WT and MT group (Additional file 4). We identified genes with the largest differences in gene functions between WT group and MT group (Table 2). Expression of seven genes with this signature was validated by real-time PCR (Additional file 5).
To characterize targets of Gli1 cis-regulation, we matched up-regulated genes among DEGs with putative Gli1-targeting genes from ChIP-chip. The result showed that Cfd was the only one gene shared in the two datasets while 32 matched genes in WT group and 4 matched genes in MT group were found (Additional file 6).
E95K mutant signaling impaired cellular migration and proliferation
According to our previous results, Gli1 targets and DEGs were enriched in cell-to-cell signaling and interaction as well as cellular activities (Tables 1 and 2). Wnt5a, an enhancer of cell migration and invasion, was down regulated in MT group. Cell migration and proliferation assay was performed to analyze the impact of altered IHH signaling on cellular activity. We measured the migration speed with high (10%) or low (0.25%) FBS concentration. We calculated the absorbance of MTT at 570 nm after 48 h of cell growth based on a gradient of IhhN protein incubation. Cells in WT group had a higher migration speed (Fig. 3) and proliferation speed (Fig. 4) than those in MT group. The ability of cellular proliferation increased when IhhN concentrations increased (Fig. 4).
Pathway analysis reveals potential pathogenesis of E95K-induced BDA1
We used Reactome platform to conduct pathway analysis and showed top 25 significant pathways (Additional file 7 and Additional file 8) . Predicted Gli1-target genes were enriched in Runx2-related pathways in MT group (6 of top 25 pathways, p < 0.05, Additional file 8) but not in WT group (Additional file 7).
The Gli1 binding position in MT group (chr17: 44,872,717-44,873,366) locates at exon1 or intron1/2 of some Runx2 transcript variants. Runx2 is a transcription factors and encodes a nuclear protein with a Runt DNA-binding domain. This protein is essential for osteoblastic differentiation and skeletal morphogenesis and acts as a scaffold for nucleic acids and regulatory factors involved in skeletal gene expression. Both loss-of-function and gain-of-function of Runx2 would result in brachydactyly. It was reported that intragenic duplication encompassing exons 3–5 or 3–6 of Runx2 gene was associated with brachydactyly phenotype [41, 42]. Real-time PCR results showed that Runx2 had a lower expression level in both WT and WT groups compared to the control group. However, Runx2 expression in MT appeared to be higher than in WT group (Additional file 9). And gene expression profile showed that Mmp13, the downstream target of Runx2, had a lower expression in WT than in MT group (Additional file 4). These results indicated that the E95K mutant IHH signaling had an impaired inhibition effect on Runx2 expression and its downstream regulation, which was likely to cause brachydactyly.
In this study, we reported that the E95K mutant IHH signaling led to a significant alteration of Gli1-DNA binding pattern and weakened cellular function of migration and proliferation. Our results suggested that E95K mutation impaired the potency of IHH signaling at transcriptional regulation. This discovery was consistent with previous report that had shown that E95K mutant IHH impaired interactions with its partners in mouse model with digit abnormalities . Our findings would be helpful to explain the mechanism of IHH/Gli1 downstream regulation altered by IHH p.E95K mutation.
In earlier studies, three mutant IHH protein, including E95K, D100E and E131K, were found in Chinese BDA1 families . Further analysis of the processing and intracellular maintenance of BDA1 mutant IhhN proteins suggested that E95K mutant IHH protein had a similar stability as wild type IHH protein . In this study, we minimized the possibility of differences caused by degradation of IHH proteins during induction. In addition, we focused on Gli1, a central positive transcriptional regulator without repressor forms among Gli genes , in order to simplify the analyses in our study. We also selected 750 nM concentration of IHH protein for signaling induction since it had been reported that activation of the signaling was most sensitive at this concentration .
Approaches for genome-scale transcription factor binding screening have been well developed in these years. The approaches had been successfully applied in researches on yeast , embryonic stem cells  and Drosophila embryos . After that, integrative analysis by combining DNA-binding sites and gene expression data advanced the precise identification of transcription regulator target genes [12, 40, 48]. In this study, we managed to combine ChIP-chip data and gene expression data to propose a Gli1-mediated model of BDA1. There were some existing bioinformatics tools for integrative analysis such as BETA (Binding and Expression Target Analysis) for ChIP and mRNA expression data . However, these methods required large datasets, such as experimental replications, to obtain credible results. Therefore, they were not applicable for a small number of specific samples in this study. It had been reported that integrative analyses like combining ChIP data with MCA analysis would greatly improve the optimal predictive rate of regulator target genes .
In this study, we failed to find Gli1 consensus binding sequence GACCACCCA from predicted peaks of ChIP-chip. This may indicated that Gli1 combined with co-transcriptors functioning as regulator in IHH/Gli1 signaling. This assumption was supported by a few reports [12, 49]. Our results suggested that transduction of E95K mutant IHH signaling altered Gli1’s interaction with transcriptional co-regulators and makes Gli1-DNA binding abnormal. Several studies had reported that graded SHH signaling resulted in distinct transcriptional regulation [50, 51]. However, there were few studies reported on how the alteration of IHH signaling affected its downstream transcriptional outputs. In this study, we demonstrated that mutant IHH signaling would lead to a distinct transcriptional pattern compared to wild type IHH signaling (Fig. 2).
We found that Runx2, a gene associated with brachydactyly, was differentially expressed between WT and MT group. Based on the evaluation of functional annotations of top 25 pathways from pathway analysis, Runx2-related pathway was shown to be one of the most relevant candidates to explain the variation in E95K-induced BDA1 and IHH downstream regulation. More functional analyses are to be conducted to confirm our discovery. Our future research plan will involve two experimental Runx2-target drugs, 3-Sulfinoalanine and O-Benzylsulfonyl-Serine, in order to rescue experiment of BDA1 cell model..
Together, by integrative analyses of Gli1-binding sites scanning, mRNA expression profiling and cellular assays, this study revealed that the E95K mutant IHH signaling significantly changed the potency and effect of its Gli1-mediated downstream transcriptional regulation and identified potential novel targets of IHH/Gli1 pathway related to chondrogenesis. We proposed explanations of how in vivo the E95K mutation alters Gli1-mediated gene regulation and finally impairs IHH signaling capacity in BDA1 model. Our result would provide helpful information to understand the role of IHH signaling in chondrogenesis and pathogenesis in BDA1.
Brachydactyly type A1
- E. coli :
Indian Hedghog N-terminus protein
Model-based Analysis of 2-Color Arrays
Transcription start site
Ingham PW, McMahon AP. Hedgehog signaling in animal development: paradigms and principles. Genes Dev. 2001;15(23):3059–87. https://doi.org/10.1101/gad.938601 PubMed PMID: 11731473.
McMahon AP, Ingham PW, Tabin CJ. Developmental roles and clinical significance of hedgehog signaling. Curr Top Dev Biol. 2003;53:1–114 PubMed PMID: 12509125.
Huangfu D, Anderson KV. Signaling from Smo to ci/Gli: conservation and divergence of hedgehog pathways from Drosophila to vertebrates. Development. 2006;133(1):3–14. https://doi.org/10.1242/dev.02169 PubMed PMID: 16339192.
Hooper JE, Scott MP. Communicating with hedgehogs. Nat Rev Mol Cell Biol. 2005;6(4):306–17. https://doi.org/10.1038/nrm1612 PubMed PMID: 15803137.
Gao B, Hu J, Stricker S, Cheung M, Ma G, Law KF, et al. A mutation in Ihh that causes digit abnormalities alters its signalling capacity and range. Nature. 2009;458(7242):1196–200. https://doi.org/10.1038/nature07862 PubMed PMID: 19252479.
Bai CB, Joyner AL. Gli1 can rescue the in vivo function of Gli2. Development. 2001;128(24):5161–72 PubMed PMID: 11748151.
Pan Y, Bai CB, Joyner AL, Wang B. Sonic hedgehog signaling regulates Gli2 transcriptional activity by suppressing its processing and degradation. Mol Cell Biol. 2006;26(9):3365–77. https://doi.org/10.1128/MCB.26.9.3365-3377.2006 PubMed PMID: 16611981; PubMed Central PMCID: PMC1447407.
Wuelling M, Vortkamp A. Transcriptional networks controlling chondrocyte proliferation and differentiation during endochondral ossification. Pediatr Nephrol. 2010;25(4):625–31. https://doi.org/10.1007/s00467-009-1368-6 PubMed PMID: 19949815.
Long F, Zhang XM, Karp S, Yang Y, McMahon AP. Genetic manipulation of hedgehog signaling in the endochondral skeleton reveals a direct role in the regulation of chondrocyte proliferation. Development. 2001;128(24):5099–108 PubMed PMID: 11748145.
Goldring MB, Tsuchimochi K, Ijiri K. The control of chondrogenesis. J Cell Biochem. 2006;97(1):33–44. https://doi.org/10.1002/jcb.20652 PubMed PMID: 16215986.
Milla LA, Cortes CR, Hodar C, Onate MG, Cambiazo V, Burgess SM, et al. Yeast-based assay identifies novel Shh/Gli target genes in vertebrate development. BMC Genomics. 2012;13:2. https://doi.org/10.1186/1471-2164-13-2 PubMed PMID: 22214306; PubMed Central PMCID: PMC3285088.
Vokes SA, Ji H, McCuine S, Tenzen T, Giles S, Zhong S, et al. Genomic characterization of Gli-activator targets in sonic hedgehog-mediated neural patterning. Development. 2007;134(10):1977–89. https://doi.org/10.1242/dev.001966 PubMed PMID: 17442700.
Wang X, Zhao Z, Muller J, Iyu A, Khng AJ, Guccione E, et al. Targeted inactivation and identification of targets of the Gli2a transcription factor in the zebrafish. Biology open. 2013;2(11):1203–13. https://doi.org/10.1242/bio.20136262 PubMed PMID: 24244857; PubMed Central PMCID: PMC3828767.
Lee EY, Ji H, Ouyang Z, Zhou B, Ma W, Vokes SA, et al. Hedgehog pathway-regulated gene networks in cerebellum development and tumorigenesis. Proc Natl Acad Sci U S A. 2010;107(21):9736–41. https://doi.org/10.1073/pnas.1004602107 PubMed PMID: 20460306; PubMed Central PMCID: PMC2906878.
Bergeron SA, Milla LA, Villegas R, Shen MC, Burgess SM, Allende ML, et al. Expression profiling identifies novel Hh/Gli-regulated genes in developing zebrafish embryos. Genomics. 2008;91(2):165–77. https://doi.org/10.1016/j.ygeno.2007.09.001 PubMed PMID: 18055165; PubMed Central PMCID: PMC2683590.
Vokes SA, Ji H, Wong WH, McMahon AP. A genome-scale analysis of the cis-regulatory circuitry underlying sonic hedgehog-mediated patterning of the mammalian limb. Genes Dev. 2008;22(19):2651–63. https://doi.org/10.1101/gad.1693008 PubMed PMID: 18832070; PubMed Central PMCID: PMC2559910.
Gao B, Guo J, She C, Shu A, Yang M, Tan Z, et al. Mutations in IHH, encoding Indian hedgehog, cause brachydactyly type A-1. Nat Genet. 2001;28(4):386–8. https://doi.org/10.1038/ng577 PubMed PMID: 11455389.
Liu M, Wang X, Cai Z, Tang Z, Cao K, Liang B, et al. A novel heterozygous mutation in the Indian hedgehog gene (IHH) is associated with brachydactyly type A1 in a Chinese family. J Hum Genet. 2006;51(8):727–31. https://doi.org/10.1007/s10038-006-0012-6 PubMed PMID: 16871364.
Ma G, Yu J, Xiao Y, Chan D, Gao B, Hu J, et al. Indian hedgehog mutations causing brachydactyly type A1 impair hedgehog signal transduction at multiple levels. Cell Res. 2011;21(9):1343–57. https://doi.org/10.1038/cr.2011.76 PubMed PMID: 21537345; PubMed Central PMCID: PMC3193471.
Stamataki D, Ulloa F, Tsoni SV, Mynett A, Briscoe J. A gradient of Gli activity mediates graded sonic hedgehog signaling in the neural tube. Genes Dev. 2005;19(5):626–41. https://doi.org/10.1101/gad.325905 PubMed PMID: 15741323; PubMed Central PMCID: PMC551582.
Bustin SA, Benes V, Garson JA, Hellemans J, Huggett J, Kubista M, et al. The MIQE guidelines: minimum information for publication of quantitative real-time PCR experiments. Clin Chem. 2009;55(4):611–22. https://doi.org/10.1373/clinchem.2008.112797 PubMed PMID: 19246619.
Brazma A, Hingamp P, Quackenbush J, Sherlock G, Spellman P, Stoeckert C, et al. Minimum information about a microarray experiment (MIAME)-toward standards for microarray data. Nat Genet. 2001;29(4):365–71. https://doi.org/10.1038/ng1201-365 PubMed PMID: 11726920.
Song JS, Johnson WE, Zhu X, Zhang X, Li W, Manrai AK, et al. Model-based analysis of two-color arrays (MA2C). Genome Biol. 2007;8(8):R178. https://doi.org/10.1186/gb-2007-8-8-r178 PubMed PMID: 17727723; PubMed Central PMCID: PMC2375008.
Edgar R, Domrachev M, Lash AE. Gene expression omnibus: NCBI gene expression and hybridization array data repository. Nucleic Acids Res. 2002;30(1):207–10 PubMed PMID: 11752295; PubMed Central PMCID: PMC99122.
Ji X, Li W, Song J, Wei L, Liu XS. CEAS: cis-regulatory element annotation system. Nucleic acids research. 2006;34(Web Server issue):W551–4. https://doi.org/10.1093/nar/gkl322 PubMed PMID: 16845068; PubMed Central PMCID: PMC1538818.
Liu T, Ortiz JA, Taing L, Meyer CA, Lee B, Zhang Y, et al. Cistrome: an integrative platform for transcriptional regulation studies. Genome Biol. 2011;12(8):R83. https://doi.org/10.1186/gb-2011-12-8-r83 PubMed PMID: 21859476; PubMed Central PMCID: PMC3245621.
Machanick P, Bailey TL. MEME-ChIP: motif analysis of large DNA datasets. Bioinformatics. 2011;27(12):1696–7. https://doi.org/10.1093/bioinformatics/btr189 PubMed PMID: 21486936; PubMed Central PMCID: PMC3106185.
Huang da W, Sherman BT, Lempicki RA. Systematic and integrative analysis of large gene lists using DAVID bioinformatics resources. Nat Protoc. 2009;4(1):44–57. https://doi.org/10.1038/nprot.2008.211 PubMed PMID: 19131956.
Huang da W, Sherman BT, Lempicki RA. Bioinformatics enrichment tools: paths toward the comprehensive functional analysis of large gene lists. Nucleic Acids Res. 2009;37(1):1–13. https://doi.org/10.1093/nar/gkn923 PubMed PMID: 19033363; PubMed Central PMCID: PMC2615629.
Ashburner M, Ball CA, Blake JA, Botstein D, Butler H, Cherry JM, et al. Gene ontology: tool for the unification of biology. The Gene Ontology Consortium. Nat Genet. 2000;25(1):25–9. https://doi.org/10.1038/75556 Epub 2000/05/10. PubMed PMID: 10802651; PubMed Central PMCID: PMCPMC3037419.
Mi H, Huang X, Muruganujan A, Tang H, Mills C, Kang D, et al. PANTHER version 11: expanded annotation data from Gene Ontology and Reactome pathways, and data analysis tool enhancements. Nucleic Acids Res. 2017;45(D1):D183–D9. https://doi.org/10.1093/nar/gkw1138 Epub 2016/12/03. PubMed PMID: 27899595; PubMed Central PMCID: PMCPMC5210595.
The Gene Ontology C. Expansion of the Gene Ontology knowledgebase and resources. Nucleic Acids Res. 2017;45(D1):D331–D8. https://doi.org/10.1093/nar/gkw1108 Epub 2016/12/03. PubMed PMID: 27899567; PubMed Central PMCID: PMCPMC5210579.
D'Eustachio P. Reactome knowledgebase of human biological pathways and processes. Methods Mol Biol. 2011;694:49–61. https://doi.org/10.1007/978-1-60761-977-2_4 Epub 2010/11/18. PubMed PMID: 21082427.
Giardine B, Riemer C, Hardison RC, Burhans R, Elnitski L, Shah P, et al. Galaxy: a platform for interactive large-scale genome analysis. Genome Res. 2005;15(10):1451–5. https://doi.org/10.1101/gr.4086505 PubMed PMID: 16169926; PubMed Central PMCID: PMC1240089.
Goecks J, Nekrutenko A, Taylor J, Galaxy T. Galaxy: a comprehensive approach for supporting accessible, reproducible, and transparent computational research in the life sciences. Genome Biol. 2010;11(8):R86. https://doi.org/10.1186/gb-2010-11-8-r86 PubMed PMID: 20738864; PubMed Central PMCID: PMC2945788.
Blankenberg D, Von Kuster G, Coraor N, Ananda G, Lazarus R, Mangan M, et al. Galaxy: a web-based genome analysis tool for experimentalists. Curr Protocols Mol Biol. 2010;Chapter 19:Unit 19 0 1–21. https://doi.org/10.1002/0471142727.mb1910s89. PubMed PMID: 20069535; PubMed Central PMCID: PMC4264107.
Wang S, Sun H, Ma J, Zang C, Wang C, Wang J, et al. Target analysis by integration of transcriptome and ChIP-seq data with BETA. Nat Protoc. 2013;8(12):2502–15. https://doi.org/10.1038/nprot.2013.150 PubMed PMID: 24263090; PubMed Central PMCID: PMC4135175.
Chen H, Boutros PC. VennDiagram: a package for the generation of highly-customizable Venn and Euler diagrams in R. BMC Bioinform. 2011;12:35. https://doi.org/10.1186/1471-2105-12-35 PubMed PMID: 21269502; PubMed Central PMCID: PMC3041657.
Agren M, Kogerman P, Kleman MI, Wessling M, Toftgard R. Expression of the PTCH1 tumor suppressor gene is regulated by alternative promoters and a single functional Gli-binding site. Gene. 2004;330:101–14. https://doi.org/10.1016/j.gene.2004.01.010 PubMed PMID: 15087129.
McDade SS, Patel D, Moran M, Campbell J, Fenwick K, Kozarewa I, et al. Genome-wide characterization reveals complex interplay between TP53 and TP63 in response to genotoxic stress. Nucleic Acids Res. 2014. https://doi.org/10.1093/nar/gku299 PubMed PMID: 24823795.
Moffatt P, Ben Amor M, Glorieux FH, Roschger P, Klaushofer K, Schwartzentruber JA, et al. Metaphyseal dysplasia with maxillary hypoplasia and brachydactyly is caused by a duplication in RUNX2. Am J Hum Genet. 2013;92(2):252–8. https://doi.org/10.1016/j.ajhg.2012.12.001 Epub 2013/01/08. PubMed PMID: 23290074; PubMed Central PMCID: PMCPMC3567275.
Al-Yassin A, Calder AD, Harrison M, Lester T, Lord H, Oldridge M, et al. A three-generation family with metaphyseal dysplasia, maxillary hypoplasia and brachydactyly (MDMHB) due to intragenic RUNX2 duplication. Eur J Hum Genet. 2018;26(9):1288–93. https://doi.org/10.1038/s41431-018-0166-7 Epub 2018/06/13. PubMed PMID: 29891876; PubMed Central PMCID: PMCPMC6117264.
Regl G, Neill GW, Eichberger T, Kasper M, Ikram MS, Koller J, et al. Human GLI2 and GLI1 are part of a positive feedback mechanism in basal cell carcinoma. Oncogene. 2002;21(36):5529–39. https://doi.org/10.1038/sj.onc.1205748 PubMed PMID: 12165851.
Guo S, Zhou J, Gao B, Hu J, Wang H, Meng J, et al. Missense mutations in IHH impair Indian hedgehog signaling in C3H10T1/2 cells: implications for brachydactyly type A1, and new targets for hedgehog signaling. Cell Mol Biol Lett. 2010;15(1):153–76. https://doi.org/10.2478/s11658-009-0040-2 Epub 2009/12/22. PubMed PMID: 20024692.
Banerjee N, Zhang MQ. Identifying cooperativity among transcription factors controlling the cell cycle in yeast. Nucleic Acids Res. 2003;31(23):7024–31 PubMed PMID: 14627835; PubMed Central PMCID: PMC290262.
Lee TI, Jenner RG, Boyer LA, Guenther MG, Levine SS, Kumar RM, et al. Control of developmental regulators by Polycomb in human embryonic stem cells. Cell. 2006;125(2):301–13. https://doi.org/10.1016/j.cell.2006.02.043 PubMed PMID: 16630818; PubMed Central PMCID: PMC3773330.
Alekseyenko AA, Larschan E, Lai WR, Park PJ, Kuroda MI. High-resolution ChIP-chip analysis reveals that the Drosophila MSL complex selectively identifies active genes on the male X chromosome. Genes Dev. 2006;20(7):848–57. https://doi.org/10.1101/gad.1400206 PubMed PMID: 16547173; PubMed Central PMCID: PMC1472287.
McDade SS, Henry AE, Pivato GP, Kozarewa I, Mitsopoulos C, Fenwick K, et al. Genome-wide analysis of p63 binding sites identifies AP-2 factors as co-regulators of epidermal differentiation. Nucleic Acids Res. 2012;40(15):7190–206. https://doi.org/10.1093/nar/gks389 PubMed PMID: 22573176; PubMed Central PMCID: PMC3424553.
Winklmayr M, Schmid C, Laner-Plamberger S, Kaser A, Aberger F, Eichberger T, et al. Non-consensus GLI binding sites in hedgehog target gene regulation. BMC Mol Biol. 2010;11:2. https://doi.org/10.1186/1471-2199-11-2 PubMed PMID: 20070907; PubMed Central PMCID: PMC2830928.
Ericson J, Rashbass P, Schedl A, Brenner-Morton S, Kawakami A, van Heyningen V, et al. Pax6 controls progenitor cell identity and neuronal fate in response to graded Shh signaling. Cell. 1997;90(1):169–80 Epub 1997/07/11. PubMed PMID: 9230312.
Gritli-Linde A, Lewis P, McMahon AP, Linde A. The whereabouts of a morphogen: direct evidence for short- and graded long-range activity of hedgehog signaling peptides. Dev Biol. 2001;236(2):364–86. https://doi.org/10.1006/dbio.2001.0336 Epub 2001/07/31. PubMed PMID: 11476578.
The authors gratefully acknowledge the support and advice by Guangyi Zhao, Ran Ran, Xiaoxu Wang, Jiajun Shi, Lianyi Zhao, Hao Wu, Yuqi Wei, Luan Chen at Shanghai Jiao Tong University in Shanghai.
This work was supported by grants from the 863 Program (2012AA02A515, 2012AA021802), the National Nature Science Foundation of China (81773818, 81273596, 30900799, 81671326), National key research and development program (2016YFC0905000, 2016YFC0905002, 2016YFC1200200, 2016YFC0906400) for design of the study and data collection and analysis, The 4th Three-year Action Plan for Public Health of Shanghai (The Project No.: 15GWZK0101), Shanghai Pujiang Program (17PJD020) for interpretation of data, Shanghai Key Laboratory of Psychotic Disorders (13dz2260500) for manuscript writing.
Availability of data and materials
All data generated or analysed during this study are included in this published article and its supplementary information files or uploaded to public database. The raw data of ChIP-chip could be obtained from NCBI Gene Expression Omnibus (GSE74022). The raw data of gene expression array could be obtained from GEO (GSE74021).
Ethics approval and consent to participate
Consent for publication
The authors declare that they have no competing interests.
Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.
ChIP Primers. Primer sequences used for ChIP validation. (XLSX 11 kb)
Peak Calling of Gli1-binding sites. ChIP-chip peak calling. (XLS 66 kb)
Predicted Gli1 target genes. ChIP-chip predicted genes. (XLS 569 kb)
Gene Expression Profile. Differentially expressed genes between three groups. (XLS 2660 kb)
Primers for Real Time PCR. Primers for gene expression validation. (XLSX 11 kb)
Matched Genes. matched genes between ChIP and DEGs. (XLSX 9 kb)
Pathway Analysis of Predicted Gli1 target genes in WT. Top 25 pathway analysis of predicted Gli1 target genes in WT. (PDF 1618 kb)
Pathway Analysis of Predicted Gli1 target genes in MT. Top 25 pathway analysis of predicted Gli1 target genes in MT. (PDF 1633 kb)
Relative expression of Runx2 with IHH concentration increasing. Relative expression of Runx2 with increased IHH concentrations in WT and MT group compared to control group. (TIFF 489 kb)