Transcriptome correlation analysis identifies two unique craniosynostosis subtypes associated with IRS1 activation

B. D. Stamper, B. Mecham, S. S. Park, H. Wilkerson, F. M. Farin, R. P. Beyer, T. K. Bammler, L. M. Mangravite, M. L. Cunningham


The discovery of causal mechanisms associated with nonsyndromic craniosynostosis has proven to be a difficult task due to the complex nature of the disease. In this study, differential transcriptome correlation analysis was used to identify two molecularly distinct subtypes of nonsyndromic craniosynostosis, termed subtype A and subtype B. In addition to unique correlation structure, subtype A was also associated with high IGF pathway expression, whereas subtype B was associated with high integrin expression. To identify a pathologic link between altered gene correlation/expression and the disease state, phosphorylation assays were performed on primary osteoblast cell lines derived from cases within subtype A or subtype B, as well as on primary osteoblast cell lines with novel IGF1R variants previously reported by our lab (Cunningham ML, Horst JA, Rieder MJ, Hing AV, Stanaway IB, Park SS, Samudrala R, Speltz ML. Am J Med Genet A 155A: 91–97, 2011). Elevated IRS1 (pan-tyr) and GSK3β (ser-9) phosphorylation were observed in two novel IGF1R variants with receptor L domain mutations. In subtype A, a hypomineralization phenotype coupled with decreased phosphorylation of IRS1 (ser-312), p38 (thr-180/tyr-182), and p70S6K (thr-412) was observed. In subtype B, decreased phosphorylation of IRS1 (ser-312) as well as increased phosphorylation of Akt (ser-473), GSK3β (ser-9), IGF1R (tyr-1135/tyr-1136), JNK (thr-183/tyr-187), p70S6K (thr-412), and pRPS6 (ser-235/ser-236) was observed, thus implicating the activation of IRS1-mediated Akt signaling in potentiating craniosynostosis in this subtype. Taken together, these results suggest that despite the stimulation of different pathways, activating phosphorylation patterns for IRS1 were consistent in cell lines from both subtypes and the IGF1R variants, thus implicating a key role for IRS1 in the pathogenesis of nonsyndromic craniosynostosis.

  • craniosynostosis
  • subtyping
  • IRS1
  • IGF1
  • metagenomics

craniosynostosis occurs in ∼1/2,500 live births and is a pathologic condition in which fusion of the calvaria occurs prematurely (8). Craniosynostosis can be subdivided into two general categories, syndromic and nonsyndromic forms. Numerous genetic mutations with known modes of inheritance have been identified in association with syndromic craniosynostosis (30); however, the pathologic mechanisms associated with nonsyndromic craniosynostosis are poorly understood. Genetic risk factors associated with nonsyndromic craniosynostosis have been identified covering a wide range of molecular functions, suggesting several independent pathologic mechanisms are at work (32).

One mechanism proposed to be involved in nonsyndromic craniosynostosis is hyperstimulation of IGF1R signaling (5). IGF1R is expressed in both bone and developing cranial sutures where it regulates cell proliferation and differentiation, most likely by mediating the biologic activity of IRS1 though IGF1 and IGF2 binding (3, 34). Both of these ligands bind directly to IGF1R instigating tyrosine kinase activity propagating signal transduction cascades that are generally considered antiapoptotic in nature (18). One of the major intracellular substrate proteins directly downstream of IGF1R is IRS1, which undergoes activating tyrosine phosphorylation upon binding to activated IGF1R (42). Signaling cascades involving IRS-1 are highly complex as IRS1 has numerous sites at which it can be phosphorylated (43). In humans, IRS1 tyrosine phosphorylation is generally activating (e.g., tyr-612, tyr-632, tyr-941), whereas serine phosphorylation is generally deactivating (e.g., ser-312, ser-337, ser-616) (45). Activating phosphorylation patterns of IRS1 then lead to the recruitment of SH2-domain-containing signal transducers that propagate IGF1R-mediated processes (45).

Further support regarding the role of IGF1R and its downstream targets in craniosynostosis is shown by the fact that individuals with extra copies of IGF1R have been shown to have the disease (28, 48). Additionally, a recent report from our laboratory identified three novel and two rare IGF1R variants in cases with craniosynostosis (5). Two of the novel variants had arginine to histidine substitutions mapping to the receptor L domain of IGF1R at either amino acid position 406 or 595. Based on structural modeling experiments these substitutions were believed to be gain-of-function, suggesting that mutation-induced hyperstimulation of IGF1R signaling is causal in a subset of craniosynostosis cases. IGF1-mediated signaling also plays a critical role in extracellular matrix-mediated focal adhesion and cell migration (1, 21), alterations in which have been proposed as a pathologic route to the disease state (41).

The goal of the present study was to identify molecularly distinct subtypes of craniosynostosis cases and to characterize the functional aberrations within these craniosynostosis subtypes. Subtypes were identified by genome-wide correlation analysis within gene expression profiles from osteoblast cell lines derived from craniosynostosis cases and normal controls. Phosphorylation and mineralization assays were then utilized to examine whether these two subtypes possessed distinct pathologic mechanisms. Cell lines exhibiting the signature associated with subtype A, termed core-craniosynostosis signature A (CC-A), demonstrated hypomineralization and high IGF-related transcript levels (IGF1, IGF2, IGFBP1, IGFBP3, IGFBP5), whereas cell lines exhibiting the signature associated with subtype B, termed core-craniosynostosis signature B (CC-B) demonstrated hyperstimulation of IGF1R-mediated Akt signaling and high integrin transcript expression (ITGA5, ITGAV, ITGB1, ITGB3). Cell lines from both subtypes, as well as those harboring IGF1R variants (R406H and R595H), all demonstrated activating phosphorylation patterns for IRS1.


Ethics statement and participant enrollment.

Written informed consent was obtained from all participants with nonsyndromic craniosynostosis, whereas a waiver of consent was obtained from the Seattle Children's Hospital institutional review board (IRB) for the anonymous control samples used in this study. This study is Health Insurance Portability and Accountability Act compliant, and independent prospective IRB approval was obtained from each participating center. Participants were enrolled as previously described in a prospective, four-center investigation of neurodevelopment among children with nonsyndromic craniosynostosis (38). Additional information on the enrollment process can be found in an earlier publication (41).

Osteoblast expansion and culture.

Calvaria samples from craniosynostosis cases (n = 199) were obtained from discarded tissues during surgical reconstructive procedures, whereas control calvaria samples (n = 50) were obtained from discarded tissues from anonymous surgical or autopsy specimens (Table 1). Harvested calvaria were washed, cleaned, sliced, and grown in 12-well plates containing 2 ml of Waymouth's medium (Sigma, W1625) supplemented with 2× antibiotic (100× Pen/Strep/Fungizone, Hyclone SV30079.01, lot JUA33955) and 10% FBS (Hyclone SH30070.03, lot ATK33398). Upon reaching confluence, cells were passaged, grown to confluence, and stored in liquid nitrogen. Upon use, thawed osteoblast lines were grown in complete Waymouth's medium. Once confluent, cells were passaged at a cell density of 175,000 cells per 25 cm2.

View this table:
Table 1.

Age information describing case and control populations used in the microarray experiment

RNA extraction and analysis.

For a more detailed description of this process, please see the following reference (41). In brief, following the plating of 175,000 cells per 25 cm2, each osteoblast cell line was grown to 75% confluence, washed, trypsinized, and collected by centrifugation. RNA was extracted from cell pellets with the Roche High Pure miRNA Isolation Kit according to the manufacturer's protocol (Roche 050080576001). Integrity of RNA samples was assessed with an Agilent 2100 Bioanalyzer (Agilent Technologies, Santa Clara, CA). We judged RNA integrity by observing distinct and sharp 18S and 28S ribosomal RNA peaks that were baseline separated; all samples had an RNA integrity number ≥8. We determined RNA quantity by measuring OD260 with a Thermo Scientific NanoDrop-1000 Spectrophotometer (Thermo Fisher Scientific, Wilmington, DE). The NanoDrop instrument was also used to determine purity of RNA samples by measuring optical density (OD) 260/280 and OD260/230 ratios. Only samples passing this stringent quality control were analyzed for transcriptomic changes using Affymetrix Human Gene 1.0 ST arrays. Expression data were normalized by Supervised Normalization of Microarrays (SNM package) (24), utilizing a model that included covariates known to influence gene expression (craniosynostosis suture type, sex, age at sample collection, hybridization batch, and cell growth rate as represented by time to reach 75% confluence). Analyses were restricted to the N genes exhibiting the greatest variance in expression across samples. All microarray data are MIAME compliant, and the raw dataset has been deposited in the MIAME-compliant Gene Expression Omnibus (GEO) database under accession number GSE27976 (

Protein harvest for phosphorylation assays.

Following the plating of 175,000 cells per 25 cm2, eight control osteoblast cell lines as well as 16 case osteoblast cell lines (eight enriched for CC-A, and eight enriched for CC-B) (Table 2) were once again grown to 75% confluence, washed, and passaged at 100,000 cells per 9.5 cm2. Upon reaching 75% confluence, growth medium was aspirated off and washed once with 2 ml of ice-cold 1× Dulbecco's phosphate-buffered saline (DPBS, Gibco, #14190). Cells were then incubated on a rotator at 4°C in 200 μl of Beadlyte Cell Signaling Universal Lysis buffer (Millipore, #43-040) for 15 min, at which point we removed cells by scraping and collected them into microcentrifuge tubes. Cells then underwent centrifugation (16,000 g) for 20 min at 4°C. Supernatant was then transferred into two tubes: one for protein quantification and the other for the phosphorylation experiments. Both aliquots were stored at −80°C until ready to use. The aliquot for protein quantitation was thawed and assayed with the Pierce BCA protein Assay Kit (#23225) according to the manufacturer's protocol.

View this table:
Table 2.

Age information describing case and control populations used in the biologic validation experiments (phosphorylation and mineralization assays)

Phosphorylation assays.

Cell culture samples prepared for phosphoprotein and total protein assays were analyzed using the custom Non-Magnetic MILLIPLEX MAP Human phosphoprotein and total protein assays for Akt/PKB, GSK3β, IGF1R, IRS1; the premixed Magnetic MILLIPLEX MAP 11-Plex Akt/mTOR Phosphoprotein Panel; the Magnetic MILLIPLEX MAP 11-Plex Akt/mTOR total protein Panel; and the Non-Magnetic MILLIPLEX MAP Phosphoprotein 8-Plex Multi-Pathway Signaling Kit (Millipore, Billerica, MA) according to the manufacturer's instructions. Activating tyrosine phosphorylation of IRS1 was measured in the IGF1R variants, whereas deactivating IRS1 serine phosphorylation was measured in subtypes A and B due to antibody availability. In brief, a threefold dilution of each sample was carried out using the manufacturer's Assay Buffer 2. Then, 25 μl of each standard from a series of fourfold dilutions and the diluted samples were incubated with the target capturing beads on a 96-well plate for 20 h at 4°C, followed by a 1 h incubation with the appropriate 1× MILLIPLEX MAP Detection Antibody at room temperature. After the washing procedure, 1× MILLIPLEX MAP Streptavidin-PE was added to each well and incubated for 30 min. Without discarding the streptavidin-PE, we added amplification buffer and incubated it for 15 min at room temperature. Samples were washed with the assay buffer, and this wash fluid was discarded using the Magnetic Separation Block for magnetic beads assays and a vacuum manifold for the nonmagnetic beads assays. Samples were resuspended in assay buffer and the fluorescence was measured using a Luminex MAGPIX (Austin, TX) for magnetic assays and a Luminex 100 for the nonmagnetic assays. Each sample was analyzed in triplicate.

Differential gene-gene correlation analysis.

In this section, we describe the methods used to identify pairs of genes whose correlation structure is significantly different between cases and controls. To accomplish this, we analyzed differences between the within condition correlations for pairs of genes. To define the values, we carried out the following procedure: for any gene pair i and j, s_ij was defined as their correlation in cases and r_ij their correlation in controls, and the difference between these two values as d_ij. Matrix D was created with element i,j equal to d_ij. Next, a bootstrap procedure was used to test the null hypothesis that the joint distribution of the absolute value of the observed correlation differences is not larger than the absolute value of correlation differences obtained after permuting the case vs. control class label. Specifically, for any randomly selected gene pair i and j, the case vs. control label was permuted, and the difference in correlations with respect to this permuted label was then calculated. This permuted statistic was called p_ij, and matrix P was similarly defined. A one-sided Kolmogorov-Smirnov (KS) test was then used to determine whether or not the joint distribution of |D| is >|P|. The resulting P value from this test (P < 1e-6) suggests the null hypothesis can be refuted, indicating that craniosynostosis has a significant effect on the correlation structure, or connectedness, between genes. Given this difference, we identified pairs of genes that underwent significant correlation changes between cases and controls. Statistically this amounts to a test of the hypothesis that d_ij = 0. To accomplish this, cutoffs for refuting this hypothesis were defined as the 5th and 95th percentiles of the joint distribution of all p_ij. This resulted in lower and upper cutoffs of −0.26 and 0.26, respectfully. Note that if d_ij is >0.26, it implies genes i and j are more correlated in cases than controls, if it is <−0.26 it implies these genes are more correlated in cases than controls; otherwise it implies there is no difference. For each gene, the number of genes with which it undergoes a significant gain in correlation in the cases was counted. Any gene with at least two such significant gains was considered significant (837 total).

Expression of subtype signatures in additional studies.

To gain insight into the potential biologic consequences of these subtypes, additional datasets enriched for expression of genes within CC-A or CC-B were identified from the 1,700 unique human microarray studies deposited within the GEO ( that used the Affymetrix HG_U133plus2 platforms. For each study, unsupervised normalization was implemented with the SNM package to normalize expression data for intensity-dependent array effects (24). The first singular value for each probe sets data was calculated and defined as the feature information content (FIC). For each study, a KS test was used to compare the FIC values for the CC-A or CC-B genes with the FIC values for all other genes. A set of 42 studies was identified in which the genes comprising CC-A and CC-B had significantly higher FIC scores than expected by chance (Table 3). Among the most significant studies was GSE12264 (9), which contains expression data derived from three osteoblast cell lines cultured in medium that promotes mineralization at four time points. Comparing the gene-wise effects of mineralization between the two groups suggested that genes in CC-A tend to have reduced expression during mineralization, whereas genes in CC-B tend demonstrate enhanced expression.

View this table:
Table 3.

GEO studies from human samples in which genes comprising CC-A were more active than would be expected by chance

Testing mineralization signature.

This section describes the approach used to measure the expression of the CC-A and CC-B genes in the craniosynostosis samples. Variations in expression of each signature were determined across samples by first mean centering each gene and then calculating two one-sided KS-test statistics. The first is a test of the hypothesis that the expression levels of CC-A genes are larger than those of the CC-B; the second tests the opposite hypothesis, namely the expression of CC-B genes is larger than the expression of CC-A. Using a hypergeometric test, we observed that samples with large KS-test statistics tended to be cases more often than expected by chance: 28 of the top 30 for CC-A > CC-B (P < 0.01) and 30 of the top 30 for CC-B > CC-A (P < 0.01). We selected eight cell lines representing each subtype for biologic analysis by randomly selecting samples exhibiting high subtype-specific KS statistics. In addition, eight cell lines were selected at random from control samples regardless of averaged subtype-specific KS statistic.

Mineralization assays.

Following the plating of 175,000 cells per 25 cm2, each osteoblast cell line was grown to 75% confluence, washed, and passaged at 20,000 cells per 1.9 cm2 in supplemented Waymouth's media. Cells were then cultured for 21 days to allow for mineralization to occur, at which point the media was carefully removed, and cells were washed with 500 μl of 1× DPBS. Cells were then fixed in 300 μl of 10% paraformaldehyde at room temperature for 15 min. The fixative was removed by gently washing cells in 1 ml of distilled water three times (5 min per wash). Following the wash steps, cells were exposed to 1 ml of Alizarin red stain for 20 min, after which excess dye was removed by three washes with deionized water. Quantitative mineralization assays were then performed using the Osteogenesis Assay Kit (Millipore, #ECM815) according to the manufacturer's protocol.


Increased IRS1 and GSK3β phosphorylation in IGF1R variants.

Phosphorylation assays were carried out to determine whether increased activation of IGF1R-mediated signaling (p-IGF1R, p-IRS1, p-Akt, and p-GSK3β) occurred in cell lines harboring purported gain-of-function (or activating) IGF1R mutations. Levels of p-IGF1R and p-Akt were unable to be determined as the fluorescent signal for these targets was indistinguishable from background. However, increased IRS-1 tyrosine phosphorylation (activating) was observed in both variants, and increased GSK3β phosphorylation at ser-9 (deactivating) was observed in Case 1 (Table 4), suggesting that a differential phosphorylation pattern was present in cell lines harboring these variants compared with control cell lines.

View this table:
Table 4.

Increased IRS1 and GSK3β phosphorylation in IGF1R variants compared with controls

Differential mRNA correlation analysis of craniosynostosis data.

We identified molecular alterations associated with disease state by assessing changes in gene-gene correlation structure within expression data derived from osteoblast cell lines derived from 199 craniosynostosis cases and 50 controls. As described in materials and methods, the observed differences between cases and controls were significantly larger than expected by chance, suggesting that transcriptional regulation is influenced by disease state. One specific example of a pair of genes with significantly larger correlations in cases than controls is illustrated in Fig. 1, in which correlation between IGF1 and IGF1R expression was observed in cell lines derived from craniosynostosis cases, but not from controls. These case-specific correlation changes were further characterized using a subset of 837 genes that demonstrated significant case-specific increases in gene-gene correlations with at least two other genes. Interestingly, this gene subset was significantly enriched for a set of 41 genes known to be involved in biologic processes altered by craniosynostosis (n = 17 of 41, P < 0.03) (Fig. 2). To assess the nature of the craniosynostosis-induced changes in correlation structure across this gene subset, we performed classification using singular value decomposition, and two distinct clusters of genes that differ in expression properties within craniosynostosis cases were identified (Fig. 3). These gene sets are referred to as the core-craniosynostosis signatures A and B (CC-A and CC-B, respectively). The gene lists comprising CC-A and CC-B can be found in Table 5. Taken together, these findings clearly illustrate that nonsyndromic craniosynostosis has a significant effect on the correlation structure, or connectedness between a subset of genes, including genes active in biological processes known to be involved in craniosynostosis, and that the affected genes cluster into two distinct groups.

Fig. 1.

Differential gene correlation for IGF1 and IGF1R, stratified by cases and controls. Illustration of the criteria that was developed to identify genes (e.g., IGF1R) that undergo significant changes in their correlation structure with another gene (e.g., IGF1) in cases compared with controls, and whose gene-gene correlation is also significant. Each circle in the figure represents the expression level for IGF1 and IGF1R in a single sample.

Fig. 2.

Heat map comparing gene expression levels between the 2 subtypes. Transcript levels of 17 genes previously shown to be altered by craniosynostosis were found to be differentially expressed between representative cases for subgroup A (n = 8) and subgroup B (n = 8) (P < 0.01); black represents upregulation, and white represents downregulation.

Fig. 3.

Identification of unique subtypes with case enrichment. Individual samples were assigned a score based on the gene activity of their signatures. The direct comparison of core-craniosynostosis signature A (CC-A) with core-craniosynostosis signature B (CC-B) illustrates the fact that these signatures not only represent 2 distinct subtypes of patients but also include a disproportionately high number of cases.

View this table:
Table 5.

Gene lists comprising each subtype signature

Differential mineralization between craniosynostosis subtypes.

To identify biologic processes that may be influenced by the alterations in gene expression associated with subtypes A and B, publicly accessible gene expression data from >1,700 studies were mined to identify additional experimental conditions under which the signature genes represented by CC-A or CC-B were highly active. These genes were highly active in 42 human studies including GSE12264 (9), a study in which expression was assessed in three osteoblast lines exposed to a mineralization-promoting medium. To assess the effect of mineralization-promoting medium exposure on the expression of CC-A and CC-B, genome-wide differential expression analysis was performed within this dataset. Expression of genes within CC-A tended to be reduced following mineralization, whereas expression of genes within CC-B tended to be increased following mineralization (P < 0.01) (Fig. 4). This result utilized a data-driven approach to establish the hypothesis that transcriptional changes in response to craniosynostosis may be related to alterations in bone mineralization. To test this hypothesis, subsets of osteoblast cell lines derived from cases were selected based on enrichment for either CC-A or CC-B (Fig. 3). Quantitative mineralization assays using Alizarin red stain were performed on osteoblast cell lines representing each subtype and controls (Table 2) (10). Results from these experiments validated the hypothesis, in that the average level of Alizarin red staining increased in the following order: subtype A < controls < subtype B (Fig. 5). However, the only significant differences in mineralization were observed between the two subtypes (P < 0.05), and subtype A and controls (P < 0.1).

Fig. 4.

Differential expression of CC-A and CC-B genes following mineralization. Expression levels were compared between 0 and 24 h time points for the CC-A and CC-B signatures. The resulting fold changes for each signature are presented as density plots. The x-axis describes the relationship between the sign of the fold change and the condition with increased expression. Expression of genes within CC-A tended to be reduced following mineralization, whereas expression of genes within CC-B tended to be increased following mineralization (P < 0.01).

Fig. 5.

Differential mineralization observed in craniosynostosis subtypes. The metaGenomics analysis identified GSE12264 in which CC-A was highly active. Data-driven statistics describing the effect of mineralization on gene expression suggested that gene expression in CC-A tends to be reduced during mineralization, whereas gene expression in CC-B tends to increase, which led to the hypothesis that subtype B would result in a hypermineralization phenotype, whereas subtype A would result in a hypomineralization phenotype. Quantitative mineralization was performed on representative subtype A, subtype B, and control lines, and measured using Alizarin red stain. Values are means ± SE of eight cell lines per group; *P < 0.05, †P < 0.1.

Phosphorylation of signaling cascades is related to craniosynostosis subtypes.

Based on the fact that IGF1R variants were associated with increased phosphorylation of IRS1-mediated Akt/GSK3β signaling, we investigated the influence of subtype on phosphorylation within this signaling cascade using the same representative osteoblast cell lines assayed in the mineralization studies (Table 2). We did not observe any significant differences in total protein concentrations between the subtypes and controls for any target except for IRS1 when comparing subtype B and controls (Table 6). Therefore, phosphoprotein levels were used as a measure of activation rather than a phosphorylated-to-total protein ratio. Subtype B demonstrated activation of Akt-mediated signaling through significant increases in activating phosphorylation of Akt, IGF1R, JNK, p70S6K, and RPS6, as well as significant decreases in phosphorylation of inactivating IRS1 (ser-312) and GSK3β (ser-9) residues (P < 0.05) (Table 6). Significant regulation with common directionality of seven targets associated with IRS1-mediated Akt activation highly implicates activation of this pathway in subtype B. Similar to subtype B, subtype A demonstrated a significant decrease in the inactivating ser-312 phosphorylation of IRS1 relative to control (P < 0.05) (Table 6). However, the widespread activation of Akt-mediated signaling observed in subtype B was not seen in subtype A, as decreased phosphorylation of p38 and p70S6K was the only significant change compared with control that was observed (P < 0.05). Not only was p70S6K phosphorylation significantly different between both subtypes and controls, but also the directionality of p-p70S6K phosphorylation was opposite for the two subtypes.

View this table:
Table 6.

Protein targets found to be significantly different between craniosynostosis subtypes and controls


Previous work from our laboratory identified two novel IGF1R variants resulting in arginine to histidine substitutions in the receptor L domain from a population of nonsyndromic craniosynostosis cases (5). Both variants resulted in mutations on the extracellular surface of IGF1R, suggesting these changes may affect ligand binding. During IGF1R signaling, autophosphorylation of tyrosine residues, namely 1131, 1135, and 1136, occurs following a conformational change induced by ligand binding (14). The phosphorylation of these residues not only enhances the tyrosine kinase activity of IGF1R but also provides a docking site for SH2 domain containing proteins such as IRS1 (6). While no significant changes in pan-tyrosine phosphorylation were observed in the IGF1R mutants, higher levels of activating IRS1 phosphorylation (pan-tyrosine) were observed. Active IRS1 is capable of binding PI3K leading to the stimulation of Akt-mediated signaling and antiapoptotic effects (17), including the deactivating GSK3β phosphorylation at ser-9 (44), which was also observed in the R406H IGF1R variant (Table 4).

Based on the fact that two downstream targets of IGF1R-mediated signaling were significantly affected in the IGF1R mutants (Table 4), it was hypothesized that a population of nonsyndromic craniosynostosis cases would share a similar pathologic mechanism. To this end, a transcriptomic dataset was mined for unique expression patterns within subpopulations of nonsyndromic craniosynostosis cases. A prior analysis of this dataset was able to identify relatively small, but highly significant differences between cases and controls in transcripts associated with extracellular-matrix-mediated focal adhesion, including significant upregulation of IGF1 expression in coronal cases (41). However, this analysis was unable to detect significantly large gene expression changes in the majority of transcripts related to IGF1R-mediated signaling. This result was not entirely unexpected based on the fact that transcript expression is not necessarily predictive of protein activation (e.g., phosphorylation) (39). To identify unique subpopulations of cases, correlation analysis was used as an alternative strategy to detect distinct patterns of transcriptomic regulation that were strongly associated with craniosynostosis cases. Following this process, two distinct patterns emerged, termed subtype A and subtype B, which were not only unique to craniosynostosis cases, but also unique from each other. This was evidenced by the fact that those cases enriched for CC-B (e.g., yellow triangles) had exceptionally low CC-A enrichment, whereas those cases enriched for CC-A (red triangles) had exceptionally low CC-B enrichment (Fig. 3).

Based on preliminary data from phosphorylation assays on novel IGF1R variants (Table 4), similar assays were performed on subtype A and subtype B (Table 6). Analysis of phosphoprotein targets within the IRS1-mediated Akt signaling cascade revealed little dysregulation of Akt-mediated signaling in subtype A (Fig. 6). Significant differences in phosphoprotein levels were only seen in two targets within this pathway, p-IRS1 (activating) and p-p70S6K (inactivating). In contrast, subtype B cases had widespread hyperstimulation of this pathway (Fig. 6). These cases showed significant increases in IGF1R tyrosine phosphorylation as well as decreased ser-312 phosphorylation of IRS1, both of which implicate activation of this pathway. Furthermore, levels of phosphorylated JNK were observed in the subtype B, a kinase known to phosphorylate IRS1 specifically at ser-312 (Table 6) (35).

Fig. 6.

Phosphorylation cascades significantly altered in craniosynostosis subtypes.

Active IRS1 binds PI3K (11, 26), which in turn produces phosphoinositol 3,4,5-triphosphate, which binds to Akt, mediating its activity (7, 47). In order for full activation, Akt must also be phosphorylated at ser-473 (13, 37), at which point Akt is capable of phosphorylating numerous target proteins (22), including GSK3β directly at serine-9 and p70S6K at thr-389 (4, 31). Phosphorylation of GSK3β at ser-9 is deactivating, thus instigating numerous downstream consequences including disinhibition of RUNX2 and β-catenin (16, 49). Increased RUNX2 and Wnt signaling has the potential to profoundly impact osteoblast development (33). In fact, disruption of either factor has been implicated as a pathologic mechanism in craniosynostosis (2, 19, 20, 25, 40, 46). Finally, phosphorylated p70S6K phosphorylates RPS6 at multiple threonine residues, upon which activated RPS6 has the potential to drastically alter the biosynthesis of translational products (36).

The fact that subtype B possessed significant increases in activating phosphorylation as well as significant decreases in deactivating phosphorylation throughout the IGF1R/IRS1/Akt/GSK3β/p70S6K/RPS6 axis strongly supports the idea that hyperstimulation of this pathway leads to the development of craniosynostosis in this subtype (Fig. 6). It is also worth noting that active and phosphorylated Akt also stimulates the activity of additional targets not assayed here, including FOXO proteins and BMP2 (15, 27). Future investigations into the role of these two downstream Akt targets is of specific interest as both proteins are known to modulate RUNX2 activity as well (15, 23, 50).

The discovery that the pathologic mechanism of craniosynostosis is likely related to Akt-mediated signaling in subtype B suggests these cases share a common cause and should be analyzed as a group. While extensive dysregulation of Akt-mediated signaling was not observed in subtype A (Fig. 6), decreased levels of p-IRS1 (ser-312) were seen in both subtypes, suggesting IRS1 may be a critical trigger for the disease state. On the other hand, subtype A had significant downregulation of the activating thr-412 phosphorylation on p70S6K, the exact opposite effect seen in the subtype B (Table 6). That cases in both subtypes are highly correlated with the disease suggests that p70S6K dysregulation may have less to do with the pathogenesis of the disease and perhaps more to do with the different correlation structure observed between the two subtypes since p70S6K plays an important role in regulating both translation and transcription (12).

In an attempt to identify a pathologic mechanism associated with the unique correlation structure of subtype A, the gene set comprising CC-A was submitted for metaGenomics analysis. In one particular study (9), CC-A was exceptionally robust, suggesting that its correlation structure would be associated with a hypomineralization phenotype. Quantitative mineralization assays on randomly selected cases within subtype A confirmed this prediction (Fig. 5). At first, a hypomineralization phenotype may seem counterintuitive as a mechanism leading to craniosynostosis. However, making the important distinction between bone mineralization, the process of mineralizing the extracellular matrix, and bone formation, which is an entirely different process, is critical (29). The fact that subtype A demonstrates decreased mineralization suggests that specific dysregulation of matrix-mediated processes may be occurring in these cases, a pathologic mechanism previously suggested by our laboratory (41).

In conclusion, the two most robust (and unique) subtypes identified through transcriptome correlation analysis were associated with those cases expressing exceptionally high levels of either IGF-related transcripts (subtype A) or integrin transcripts (subtype B). While differential expression is of definite interest, it is important to recognize that subtype A and subtype B were identified based on gene correlation structure and not transcript levels. Based on the hypothesis that several independent mechanisms lead to the occurrence of several different forms of craniosynostosis, correlated transcriptomic data were used to classify cases into relatively exclusive subtypes that would provide insight into the pathogenesis the disease. The unique correlation structures of CC-A and CC-B revealed two unique mechanisms related to the disease state in a subset of craniosynostosis cases. Hyperstimulation of IRS1-mediated Akt signaling was found to be significantly associated with subtype B, whereas hypomineralization was associated with subtype A. Despite different phosphorylation levels, it is worth reiterating that both subtypes demonstrated a decrease in the inhibitory phosphorylation of IRS1 at ser-312, which would implicate IRS1 activation as a feature in the majority of craniosynostosis cases. Future work will focus on the precise regulatory mechanisms associated with the dysregulation in gene correlation that was observed and investigate its impact on developing craniosynostosis.


This work was supported by National Institutes of Health Grants NIH/NIDCR R01 DE-018227, NIH/NIEHS P30ES-07033, and NIH/NICHD P30HD-02274; the Jean Renny Endowment for Craniofacial Research; and a grant from the Washington State Life Sciences Discovery Fund.


No conflicts of interest, financial or otherwise, are declared by the author(s).


Author contributions: B.D.S., B.M., S.S.P., F.M.F., R.P.B., T.K.B., L.M.M., and M.L.C. conception and design of research; B.D.S., B.M., S.S.P., and H.W. performed experiments; B.D.S., B.M., S.S.P., H.W., R.P.B., T.K.B., L.M.M., and M.L.C. analyzed data; B.D.S., B.M., R.P.B., T.K.B., L.M.M., and M.L.C. interpreted results of experiments; B.D.S., B.M., L.M.M., and M.L.C. prepared figures; B.D.S. drafted manuscript; B.D.S., B.M., S.S.P., L.M.M., and M.L.C. edited and revised manuscript; B.D.S., B.M., S.S.P., L.M.M., and M.L.C. approved final version of manuscript.


  1. 1.
  2. 2.
  3. 3.
  4. 4.
  5. 5.
  6. 6.
  7. 7.
  8. 8.
  9. 9.
  10. 10.
  11. 11.
  12. 12.
  13. 13.
  14. 14.
  15. 15.
  16. 16.
  17. 17.
  18. 18.
  19. 19.
  20. 20.
  21. 21.
  22. 22.
  23. 23.
  24. 24.
  25. 25.
  26. 26.
  27. 27.
  28. 28.
  29. 29.
  30. 30.
  31. 31.
  32. 32.
  33. 33.
  34. 34.
  35. 35.
  36. 36.
  37. 37.
  38. 38.
  39. 39.
  40. 40.
  41. 41.
  42. 42.
  43. 43.
  44. 44.
  45. 45.
  46. 46.
  47. 47.
  48. 48.
  49. 49.
  50. 50.
View Abstract