Unique sex-based approach identifies transcriptomic biomarkers associated with non-syndromic craniosynostosis.

BACKGROUND: The premature fusion of one cranial suture, also referred to as non-syndromic craniosynostosis, most commonly involves premature fusion of the sagittal, coronal, or metopic sutures, in that order. Population-based epidemiological studies have found that the birth prevalence of single-suture craniosynostosis is both suture- and sex-dependent.
METHODS: Transcriptomic data from 199 individuals with isolated sagittal (n = 100), unilateral coronal (n = 50), and metopic (n = 49) synostosis were compared against a control population (n = 50) to identify transcripts accounting for the different sex-based frequencies observed in this disease.
RESULTS: Differential sex-based gene expression was classified as either gained (divergent) or lost (convergent) in affected individuals to identify transcripts related to disease predilection. Divergent expression was dependent on synostosis sub-type, and was extensive in metopic craniosynostosis specifically. Convergent microarray-based expression was independent of synostosis sub-type, with convergent expression of FBN2, IGF2BP3, PDE1C and TINAGL1 being the most robust across all synostosis sub-types.
CONCLUSIONS: Analysis of sex-based gene expression followed by validation by qRT-PCR identified that concurrent upregulation of FBN2 and IGF2BP3, and downregulation of TINAGL1 in craniosynostosis cases were all associated with increased RUNX2 expression and may represent a transcriptomic signature that can be used to characterize a subset of single-suture craniosynostosis cases.

Introduction

Craniosynostosis is the pathologic fusion of calvarial bones and occurs in approximately 1/2500 live births. A genetic component to the disease is likely given a 7%–10% recurrence risk.1 Craniosynostosis cases can be classified into two categories, syndromic and non-syndromic forms. Over one hundred forms of syndromic craniosynostosis have been identified, including Apert, Crouzon, Muenke, Pfeiffer, and Saethre-Chotzen syndromes, which are associated with various phenotypic manifestations in addition to premature fusion of calvarial sutures.2 In contrast, individuals with non-syndromic craniosynostosis lack non-sutural phenotypes.

Mutations in genes such as FGFR1-3, TWIST1, EFNB1, FBN1, MSX2, RAB23, RECQL4, and TGFBR1-2 have been associated with syndromic craniosynostosis,3 whereas markers for the pathogenesis of non-syndromic forms of the disease have proven more difficult to identify.3–5 Improving our understanding of what causes non-syndromic craniosynostosis remains an important endeavor considering non-syndromic forms of this disease account for approximately 85% of all cases. In fact, recent reports suggest the incidence of all forms of non-syndromic craniosynostosis, especially metopic cases, is on the rise.6,7 Generally, half of non-syndromic cases involve premature fusion of the sagittal suture, whereas premature coronal and metopic suture closure occurs in 22% and 15% of cases, respectively.8

The differences between syndromic and non-syndromic forms of craniosynostosis including clinical features, classified mutations, and incidence suggest that these two forms of craniosynostosis are unique. In this paper, another characteristic difference between syndromic and non-syndromic craniosynostosis is investigated, namely sex predilection. In non-syndromic craniosynostosis, males have an overall higher frequency of developing craniosynostosis than females. In isolated sagittal and metopic craniosynostosis there is a near four-fold increased incidence among males, whereas in coronal cases the male/female ratio is nearly even.8 The observation that maleness has been shown to predispose individuals to craniosynostosis suggests that if transcriptomic changes are driving this predilection, then females with male-patterned expression of select disease- related genes may be predisposed to developing non-syndromic craniosynostosis.

To this end, transcriptomic arrays from a large cohort of individuals with non-syndromic craniosynostosis were analyzed in order to identify sex-related changes in gene expression that predispose individuals to developing this disease. First, differential gene expression between control males and females was compared against differential gene expression between males and females with craniosynostosis in order to create divergent and convergent gene sets. Sex related expression was defined as divergent in craniosynostosis cases when sex differences were not observed in controls and convergent when sex differences present in controls were absent in cases (Fig. 1). Next, the divergent and convergent gene expression of affected male and female cases was compared directly against sex-matched controls to confirm that the gain or loss of transcript expression was significantly different (Fig. 1). Results from these comparisons identified sex-dependant disturbances to genes involved in Ca2+-mediated phosphatidylinositol 3-kinase/protein kinase B (PI3K/Akt) signaling as key targets that predispose females to craniosynostosis, whereas TINAGL1 expression plays a more important role in predisposing males to the disease by potentially disrupting TGF-β activity. More importantly, differential expression of genes related to these pathways can be directly linked to increased RUNX2 activity, a mechanism known to cause premature fusion of calvarial sutures.

Figure 1

Comparative approach used to identify transcripts that may predispose individuals to non-syndromic craniosynostosis based on sex.

Notes: Sex-based differences associated with affected individuals and not seen in controls are contained within the divergent gene set, whereas sex-based differences unique to the control population and not seen in affected individuals are contained within the convergent gene set. Those transcripts contained within the shared gene set are sex-based differences associated with both controls and affected craniosynostosis cases. Once identified, gene expression considered divergent or convergent was reanalyzed by comparing affected individuals to their respective sex-matched controls.

Methods

Informed consent

Written informed consent was obtained from all participants with single-suture 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 HIPAA compliant, and we obtained independent prospective IRB approval from each participating center, including Seattle Children’s Hospital, Northwestern University in Chicago, Children’s Heath Care of Atlanta, and St. Louis Children’s Hospital.

Cell culture

All primary osteoblast cell lines were derived from participants in a previously described craniosynostosis study.9 Subsets of the 249 cell lines included 50 controls and 100 sagittal, 50 coronal, and 49 metopic craniosynostosis cases. Craniosynostotic calvaria were obtained from discarded tissues during surgical reconstructive procedures, whereas control calvaria were obtained from discarded tissues from anonymous surgical or autopsy specimens. Harvested calvaria samples were then washed and expanded in Waymouth’s media (Sigma: St. Louis, MO) supplemented with 2X Penicillin/Streptomycin/Fungizone (Hyclone: Logan, UT) and 10% heat-inactivated fetal bovine serum (FBS) (Hyclone: Logan, UT). Osteoblasts were grown at 37 °C, 5% CO2, and 99% humidity, trypsinized using 0.05% Trypsin (Hyclone: Logan, UT) upon reaching 75% confluence, counted and passaged at a cell density of 175,000 cells per 25 cm2.

Cell harvest and RNA isolation

Once the re-plated cells reached 75% confluence, they were photographed for quality control purposes, washed twice with 1X PBS, and trypsinized. An equal volume of media containing FBS was added after trypsin exposure, and cells were centrifuged at 200 × g for 10 minutes at 4 °C in nuclease free 15 mL conical tubes (Corning: Lowell, MA). Following a washing step, cells were centrifuged again at 200 × g for 10 minutes at 4 °C. RNA extraction from cell pellets was performed using the Roche High Pure miRNA Isolation Kit in accordance to the manufacturer’s protocol (Roche: Indianapolis, IN). RNA was stored immediately in −80 °C and either analyzed by quantitative reverse transcriptase PCR (qRT-PCR) or submitted for microarray processing on dry ice.

Preparation of samples for quantitative real-time PCR

cDNA was synthesized from total RNA using the High RevertAid First Strand cDNA Synthesis Kit (Fermentas: Glen Burnie, MD) according to the manufacturers protocol. Reactions were set up in duplicate with the following components in each well of an ABI microAMP Fast Optical 96-well Reaction Plate (Life Technologies: Carlsbad, CA), 10 μL 2X SensiMix SYBR low-ROX master mix (Bioline: London, UK), 2 μL primers (1–5 μM) (Sigma: St. Louis, MO), and 8 μL sample cDNA (1.25 ng/μL). All primers were designed using PrimerBank (http://pga.mgh.harvard.edu/primerbank/) to amplify the following targets: CD24 (NM_013230), FBN2 (NM_001999), IGF2BP3 (NM_006547), PDE1C (NM_005020), RUNX2 (NM_001024630), and TINAGL1 (NM_022164). Reaction plates were then sealed and centrifuged for 1 minute at 1,000 × g. qRT-PCR was performed on an Applied Biosystems ABI7500 Fast Real-Time PCR System (Life Technologies: Carlsbad, CA). Ct values were normalized to 18s expression (NR_003286) and relative differences in mRNA expression were compared between craniosynostosis cases and controls.

Validation of divergent CD24 expression by qRT-PCR

Quantitative real-time PCR was performed on 32 randomly selected primary cell lines (24 craniosynostosis cases and 8 controls). Eight coronal, eight metopic and eight sagittal cases comprised the twenty-four total craniosynostosis cases assayed. Within each group of eight, four females and four males were assayed in order to determine whether divergent expression was consistent among all craniosynostosis sub-types. Similarly, of the eight controls assayed, four were female and four were male.

qRT-PCR on cell lines expressing high FBN2, high IGF2BP3, and low PDE1C

Quantitative real-time PCR was performed on primary cell lines in which microarray-based expression of both FBN2 and IGF2BP3 were in the top third of the entire cohort (n = 249), while maintaining PDE1C expression within the bottom third (n = 21; 17 craniosynostosis cases and 4 controls).

Microarray analysis

Raw microarray data were pre-processed and normalized with Affymetrix® Expression Console™ Software using RMA normalization (http://affymetrix.com). Microarray quality control metrics include the manufacturer’s recommended guidelines: visual inspection of probe array images, proper ranking of hybridization and Poly-A controls, and area under the curve values for a receiver operating characteristic plot comparing the positive control and negative control signal values. Other microarray quality control metrics were generated with the Bioconductor R package named “aroma.affymetrix” and included the relative log expression (RLE) values (assessing the potential spread or shifting of expression values), and the normalized unscaled standard errors (NUSE) (assessing variability of genes across arrays) (R: A Language and Environment for Statistical Computing, R Foundation for Statistical Computing, http://www.R-project.org).10 All microarray data have been deposited in the Gene Expression Omnibus Database under accession number GSE27976 (http://www.ncbi.nlm.nih.gov/geo/).

Statistical analysis

From the normalized data, genes with significant evidence for differential expression were identified using the limma package11 in Bioconductor.12 P-values were calculated with a modified t-test in conjunction with an empirical Bayes method to moderate the standard errors of the estimated log-fold changes. P-values were adjusted for multiplicity using Bioconductor’s implementation of the Benjamini-Hochberg13 P-value adjustment method, which allows for selecting statistically significant genes while controlling the estimated false discovery rate.

Results

Identification of divergent, convergent, and shared transcripts

249 primary osteoblast cell lines were developed from individuals with non-syndromic craniosynostosis and controls. Transcriptomic analysis of these cell lines identified sex-based differential gene expression that persisted during in vitro culture. Divergent genes unique to craniosynostosis were highly dependent on the synostosis sub-type, in that the number of sex-related, craniosynostosis-related transcripts was variable among coronal, metopic, and sagittal cases (Table 1). In fact, only CD24 downregulation was considered significant overall three sutures (Supplemental Table 1). Divergent CD24 expression was verified by qRT-PCR (Supplemental Table 2), confirming female cases had lower expression of CD24 than their male counterparts (fold change = −2.8, P = 0.06) with no significant differences seen between female and male controls (P = 0.61). Expression of COL11A1, F2R, FAM38B, PLA2G16, and VCAM1 was consistent among coronal and metopic cases, whereas significant expression of all other divergent transcripts was unique to one specific synostosis sub-type (Supplemental Table 1). These data suggest that divergent sex-related gene expression in non- syndromic craniosynostosis is generally suture-dependent, and that divergent expression is most prominent in metopic cases since 99 transcripts were identified in this synostosis sub-type alone (Table 1).

Table 1

Number of divergent, shared, and convergent transcripts identified based on craniosynostosis sub-type (gene symbols for the transcripts represented here are listed in Table 2).

	Transcripts in divergent gene set	Transcripts in shared gene set	Transcripts in convergent gene set
Coronal	6	21	31
Metopic	99	21	31
Sagittal	1	20	32

Unlike the divergent gene set, transcript expression within the convergent gene set was not dependent on the synostosis sub-type. This was evidenced by the fact that sex-related expression in 30 of the 32 convergent genes was consistent among all three synostosis sub-types (Table 2). These results suggest that losing sex-based differences in transcriptional activity is consistent among affected individuals and is therefore independent of synostosis sub-type. Like the convergent gene set, transcript expression within the shared gene set was also not dependent on the synostosis sub-type. This was evidenced by the fact that in 20 of the 22 genes were shared among controls and all three forms of craniosynostosis (Table 2). Furthermore, of the 20 genes shared among controls and all three craniosynostosis types, all 20 are sex-linked (Supplemental Table 3). GAGE12C and a pseudogene of AARSD1 are the only X-linked genes in the shared set and are upregulated in females, whereas the remaining 18 Y-linked genes are upregulated in males (Supplemental Table 3). By definition, the shared gene set includes genes in which expression was found to be differentially regulated between control males and females as well as case males and females. While this gene set may not provide insight in differential sex-related gene expression between cases and controls, it does validate our methodology in that differential expression of X- and Y-linked genes is consistent with respect to the sex of the individual.

Table 2

Divergent, shared and convergent transcripts that underwent large (|log2FC| >1.5) and significant (P < 0.05) changes when comparing either male and female craniosynostosis cases directly or male and female control populations directly.

Divergent gene set	Shared gene set	Convergent gene set
Coronal cases (6):CDH6, FAM38B, NEFM, OGN, S100A4, VCAM1Metopic cases (99):ABCA8, ADAMTS5, ADAMTSL3, AK5, ANLN, ARHGAP11 A, ARHGDIB, ASPM, AURKA, BRIP1, BUB1, BUB1B, CA12, CASC5, CCNA2, CCNB1, CCNB2, CCNE2, CD24, CDC2, CDC20, CDCA2, CDKN3, CENPF, CENPI, CEP55, CLDN1, CLDN11, COL11A1, COL14A1, CRISPLD1, CTSC, CYTL1, DEPDC1, DLGAP5, DTL, ESCO2, EXO1, EYA4, F2R, F2RL2, FAM111B, FAM38B, FLG, FNDC1, HAPLN1, HELLS, HEY2, HIST1H1B, HIST1H2BM, HIST1H3B, HJURP, HMCN1, HMMR, IFI30, KIAA0101, KIF11, KIF14, KIF15, KIF20A, KIF20B, KIF2C, LMNB1, LOC100289612, MELK, MKI67, MOXD1, MXRA5, MYCT1, NCAPG, NEIL3, NUF2, NUSAP1, PAMR1, PBK, PLA2G16, PLK1, PLK4, PLXDC2, PRC1, PRR11, PTGER2, PTTG1, RAD51AP1, RRM2, SCRG1, SGCG, SGOL1, SHCBP1, SKA1, SKA3, SLITRK6, SPC25, SULF2, TOP2 A, TRIP13, TTK, VCAM1, ZIC1Sagittal cases (1):CD24	All cases (20):BPY2, CYorf15A, CYorf15B, DDX3Y, EIF1AY, GAGE12C, KDM5D, JPX, TTTY14, NLGN4Y, PRKY, RBMY1A1, RBMY1B, RBMY2EP, RPS4Y1, RPS4Y2, TSPY1, USP9Y, UTY, ZFY	All cases (30):ACAN, ACTG2, CNTNAP3, CHI3L1a, COL4A1, COL4A2, DSG2b, ENPP2, FBN2, FGL2, FLT1, GCNT4, GFRA1, GREM2, ID4, IGF2BP3, LAMC2, LOXL4, LPPR4, LUZP2, MCAM, MYOCD, OXTR, PDE1C, PENK, PLA2G5, SEMA3C, SEMA3D, SFRP4, SLC7A2, SULT1E1, TINAGL1

Notes:

Convergent CHI3L1 expression was observed in metopic and sagittal cases only;

Convergent DSG2 expression was observed in coronal and sagittal cases only.

Divergent genes differentially expressed from controls

The identification of divergent synostotic expression may provide insight into genes that shape the observed sex frequencies found in the disease state, however, it is also of interest to investigate whether divergent synostotic expression differs significantly from control populations. To this end, affected individuals were separated based on sex and compared against their respective sex-matched controls to determine whether divergent synostotic gene expression differed from controls to a significantly large extent (|fold change| >1.5, P < 0.05) (Supplemental Table 4). VCAM1 was the only divergent transcript differentially expressed between sex-matched cases and controls to a large and significant extent in males with all forms of non-syndromic craniosynostosis (Table 3). It is also worth noting that increased VCAM1 expression was observed in female coronal and sagittal cases compared to female controls as well, suggesting that VCAM1 upregulation may not be a male-specific craniosynostotic effect (Table 3).

Table 3

Genes differentially expressed between all affected males or all affected females compared to their sex-matched controls (|% change| >50%, P < 0.05).

Transcript	Gene set	Female_male (% change)				Female_female (% change)			Male_male (% change)
		CON_CON	COR_COR	MET_MET	SAG_SAG	COR_CON	MET_CON	SAG_CON	COR_CON	MET_CON	SAG_CON
FBN2	Convergent	−76	–	–	–	95	107	57	–	–	–
IGF2BP3	Convergent	−88	–	–	–	88	157	88	–	–	–
PDE1C	Convergent	74	–	–	–	−108	−96	−54	–	–	–
TINAGL1	Convergent	−71	–	–	–	–	–	–	−113	−108	−52
CD24	Divergent	–	−49a	−56	−75	–	−76	−95	–	–	–
VCAM1	Divergent	–	−64	−55	–	78	–	141	103	61	85

Note:

Significant difference, but below % change threshold value.

Convergent genes differentially expressed from controls

Analogous to the analysis performed on the divergent gene set, affected individuals were separated based on sex and compared against sex-matched controls to determine whether convergent gene expression in craniosynostosis differed from controls to a significantly large extent. It was expected that transcript expression would be more consistent between sex-matched cases and controls based on the fact that convergent expression was independent of synostosis sub-type. Four convergent genes (IGF2BP3, FBN2, PDE1C and TINAGL1) lost significantly large differences in expression between males and females in all cases, while still maintaining large (|fold change| > 1.5) and significant (P < 0.05) differences in expression from sex-matched controls (Fig. 2). Compared to sex-matched controls, female cases had lower expression of PDE1C and higher expression of IGF2BP3 and FBN2, while male cases had decreased expression of TINAGL1. The identification of these four transcripts where sex-specific expression was lost, yet differential expression from sex-matched controls was maintained, may indicate critical pathogenic roles for these genes. In other words, male and female cases are becoming more similar to each other, while simultaneously distancing themselves from control-like expression patterns where sex differences once existed. Finally, all four of these transcripts are capable of modulating RUNX2 expression/activity. Because of this, the dataset was mined to determine whether sex-based differences in RUNX2 expression were present. The search revealed RUNX2 to be a significantly convergent transcript that did not meet our threshold cutoff (|% change| < 50%), as female controls expressed 32% more RUNX2 than their male counterparts (P < 0.05). No significant sex-based differences were observed when comparing males and females in each of the craniosynostosis sub-types.

Figure 2

Consistently convergent gene expression.

Notes: Convergent transcripts were defined as those genes that were differentially expressed between male and female controls (|fold change| >1.5, P < 0.05), but lost these sex-related differences when male and female cases were compared. Four genes fit these criteria and were consistent in all three synostosis sub-types (coronal, metopic and sagittal) to a large and significant extent.

Quantitative real-time PCR on cases with high FBN2, high IGF2BP3, and low PDE1C expression

qRT-PCR was performed on 21 primary cell lines (17 craniosynostosis cases and 4 controls) that were identified in the microarray study to have high FBN2, high IGF2BP3, and low PDE1C levels (ie, female cases becoming male-like) (Fig. 2 and Supplemental Table 5). These samples were specifically chosen to investigate whether consistently convergent gene expression in craniosynostosis cases leads to dysregulation of RUNX2 expression/activity, due to the fact that the three genes are known to modulate its activity. First, convergent expression of FBN2, IGF2BP3, PDE1C, RUNX2, and TINAGL1 was confirmed in the selected cases by qRT-PCR, based on the fact that female and male cases showed no significant differences in the expression of these transcripts (Supplemental Table 5). When compared to controls, qRT-PCR was unable to detect significant downregulation of PDE1C in cases suggesting concurrent upregulation of FBN2 and IGF2BP3 is more critical to the disease state. This was found to be true in that three-fold FBN2 upregulation and four-fold IGF2BP3 upregulation, in addition to nearly five-fold TINAGL1 downregulation were observed in cases compared to controls (Table 4). Furthermore, qRTPCR not only confirmed that these same case lines expressed RUNX2 to a greater extent than matched controls, but also found this upregulation to be more than six-fold (Table 4). These results suggests that the ↑FBN2, ↑IGF2BP3, ↓TINAGL1 transcriptomic signature represents a polygenic signature (and potential cause) of craniosynostosis that acts through cooperative RUNX2 upregulation.

Table 4

qRT-PCR validation of differential expression between craniosynostosis cases and controls both concurrently expressing the high FBN2, high IGF2BP3 and low PDE1C (P < 0.05).

Transcript	Case_control (fold change)	Primers
FBN2	3.2	fwd: CTGAAGGCGGGTTTCTAGCGrvs: CAAATCGGGACAATGCACTGG
IGF2BP3	4.5	fwd: TATATCGGAAACCTCAGCGAGArvs: GGACCGAGTGCTCAACTTCT
RUNX2	6.5	fwd: TGGTTACTGTCATGGCGGGTArvs: TCTCAGATCGTTGAACCTTGCTA
TINAGL1	−4.8	fwd: ATGGGACCCACTCAGTCAAGArvs: GTTGGCCGCAGTCCAGTATTT

Discussion

Unlike the majority of syndromic forms, nonsyndromic craniosynostosis is much more common, is not associated with predictable non-calvarial phenotypic manifestations, and has unique sex-based frequencies depending on the suture that is prematurely fusing.2,8 One notable exception is the P250R fibroblast growth factor receptor 3 (FGFR3) mutation in Muenke syndrome, which is associated with a more severe phenotype in females compared to males.14 To investigate sex-based differences in the frequency of non-syndromic craniosynostosis, a unique comparative transcriptomic approach was developed to identify transcriptional changes driving sex predilection (Fig. 1). This model was necessary because there are sex specific effects differentially expressed among controls that are not related to craniosynostotic mechanisms (ie, genes on the sex chromosomes and those genes regulated by genes on the sex chromosomes). These normal differences need to be “removed” in order to understand sex related effects that are specific to the disease state. As proof of principle, our approach was internally validated by the fact that only sex-linked genes were identified in the shared region of the Venn diagrams (Supplemental Table 3).

Sex-related transcript expression that “diverged” in craniosynostotic individuals was dependent on the synostosis sub-type. In fact, extensive divergent expression was exclusive to metopic cases (Table 1). It was expected that metopic and sagittal cases would share similar divergent gene expression profiles considering metopic and sagittal craniosynostosis not only have defined male to female ratios greater than 3:1, but also share many other epidemiological similarities.8 It is worth noting that expression of one transcript (CD24) was found to be divergent in both sagittal and metopic cases, suggesting that reduced levels of CD24 expression may be associated with developing craniosynostosis in females (Table 3). In fact, recent evidence suggests that CD24 may inhibit cell invasiveness,15 suggesting that loss of CD24 expression may promote cell invasion, perhaps at the osteogenic front, and act as a mechanism for premature suture closure.

In males, divergent upregulation of VCAM1 was observed in coronal and metopic cases, but not in sagittal cases. However, both male and female sagittal cases had increased VCAM1 expression compared to their sex-matched controls (Table 3). Even though VCAM1 expression was divergent in coronal and metopic cases and upregulated in all male cases compared to sex-matched controls, the increase seen in affected females with coronal and sagittal synostosis makes it unlikely that upregulation of this transcript is a male-specific effect. While male cases generally have higher VCAM1 levels than female cases, female coronal and sagittal cases still have higher VCAM1 expression than their sex-matched controls to a significantly large extent (Table 3). These results suggest that while increased VCAM1 expression is related to craniosynostosis, its upregulation may be sex-independent. This observation substantiates an earlier report from our lab identifying VCAM1 upregulation as a novel biomarker in single-suture craniosynostosis.16

Unlike divergent expression, convergent expression in craniosynostotic individuals was independent of the synostosis sub-type (Table 2). Convergent transcripts were then mined to identify gene expression that was also significantly different between the cases and their sex-matched controls (Supplemental Table 4). Only three of the convergent transcripts were differentially expressed in all three synostosis sub-types and confirmed by qRT-PCR; IGF2BP3 and FBN2 in females, and TINAGL1 in males (Table 3). Convergent expression of IGF2BP3 and FBN2 is interesting on two levels. First, gain of IGF2BP3 and FBN2 expression in female cases results in “male-like” expression of these transcripts (Fig. 2). The fact that males have a higher incidence for developing craniosynostosis would suggest that females with “male-like” expression of craniosynostosis-associated genes may have a higher risk of developing the disease. Second, increased IGF2BP3 and FBN2 expression in case females would have a the same effect on RUNX2 activity, a transcription factor that is expressed in fusing cranial sutures.17 Mutations of the RUNX2 modulator, TWIST1, are associated with a syndromic form of craniosynostosis known as Saethre-Chotzen syndrome.18 TWIST1 is a negative regulator of RUNX2, therefore loss of function mutations in TWIST1 lead to increased RUNX2 activity accounting for the presumed mechanism of action by which these mutations lead to non-syndromic craniosynostosis.19 This mechanism is consistent with the finding that duplication of the RUNX2 gene is associated with metopic craniosynostosis.20

In osteoblasts, RUNX2 activation can also occur through Ca2+-mediated PI3K/Akt activation and subsequent FOXO1/4 phosphorylation.21,22 The male-like expression of IGF2BP3 and FBN2 in female craniosynostosis cases can impact this pathway at various points to promote RUNX2 activity (Fig. 3). Increased IGF2BP3 expression promotes IGF2 translation leading to PI3K/Akt-mediated RUNX2 activity.23 Upregulation of FBN2 has been shown to increase bone mass by restricting RUNX2 microRNA expression and processing, thereby attenuating the translational repression of RUNX2.24 Another mechanism by which increased FBN2 expression leads to increased RUNX2 activity is through the negative regulation of TGF-β signaling, a cascade known to repress RUNX2 activity.25,26 In fact, type 1 Loeys-Dietz syndrome, which is caused by mutations in TGF-β receptors, is associated with craniosynostotic phenotypes.27 Finally, high FBN2 expression may increase osteoblast migration, as alterations to intracellular Ca2+ levels have been shown to promote cell migration along fibrillin paths.28 Mutations in FBN2 cause congenital contractual arachnodactyly,29 and while there are currently no known FBN2 mutations that cause craniosynostosis, mutations in FBN1 cause Marfan Syndrome30 and possibly Shprintzen-Goldberg Syndrome,31,32 a syndromic form of craniosynostosis with Marfanoid features that is associated with arachnodactyly.29,33 The fact that alterations to FBN1 and FBN2 gene function lead to similar phenotypes, coupled with the fact that FBN1 and FBN2 share overlapping functionality,24,34 strongly implicates the role fibrillin dysregulation in the pathogenesis of craniosynostosis. Taken together, craniosynostosis cases with high FBN2 and IGF2BP3 expression compared to controls with high expression of these same transcripts demonstrated a six-fold increase in RUNX2 expression (Table 3). This suggests that significantly high FBN2 and IGF2BP3 expression coupled with RUNX2 upregulation might act as a pathogenic mechanism for craniosynostosis (Fig. 3).

Figure 3

Proposed RUNX2 activation model for differentially expressed convergent transcripts.

With respect to convergent expression that was differentially expressed between male cases and male controls, only TINAGL1 downregulation surfaced as a potential marker associated with all male forms of nonsyndromic craniosynostosis. TINAGL1 is a matricellular protein that has been shown to promote cell adhesion and invasion through its interactions with cell surface receptors and cell matrix proteins such as laminin and integrin α5-β1.35,36 Our previous findings implicating extracellular matrix-mediated focal adhesion as the strongest signal associated with non-syndromic craniosynostosis cases strongly supports the identification of TINAGL1 as a transcript associated with the disease state.16 Furthermore, TINAGL1 overexpression has been shown to couple TGF-β to SMAD2/3 phosphorylation in endothelial cells, a mechanism by which TGF-β attenuates RUNX2 function in osteoblasts.37,38 Presumably, lower levels of TINAGL1 observed in craniosynostotic males would lead to decreased TGF-β signaling through SMAD2/3, and increased RUNX2 function (Fig. 3).37,38 This mechanism is further supported by the fact that cases with high FBN2 and IGF2BP3 expression have low TINAGL1 levels in addition to high RUNX2 expression (Table 4).

In conclusion, non-syndromic craniosynostosis is a complex disease with multiple causes, of which only a handful have been discovered. By utilizing this unique sex-based transcriptomic approach, our lab was able to generate novel candidate genes that predispose individuals to developing the disease state. Based on the likelihood that non-syndromic craniosynostosis is a polygenic disease, sub-populations within each synostosis sub-type are likely driving the significant changes in gene expression that surfaced using this approach. Finding commonality among these diverse targets is critical to understanding the mechanism by which premature fusion of calvarial sutures occurs. Our approach was able to identify convergent sex-based gene expression that cooperatively upregulates RUNX2 expression (Table 4), implicating the ↑FBN2, ↑IGF2BP3, ↓TINAGL1 transcriptomic signature in the development of single-suture craniosynostosis. A potential limitation of this study is that while care was taken to age-match the case and control samples, the primary control lines in the entire cohort were generally from older children due to the limited availability of calvaria samples from controls less than one year of age.16 Future work will focus on large-scale sequencing of the candidates presented here, as well as any transcriptional regulators (promoter regions, transcription factors, nuclear receptors, etc.) in individuals that are known to have divergent/convergent transcriptomic profiles to identify causal mutations. The identification of transcripts associated with osteogenic differentiation and proliferation in this study suggests that these candidates are of interest, and that this novel approach has potential utility in the identification of gene candidates and transcriptomic signatures in other diseases with disparate sex-based frequencies.

Supplementary Table 1

(divergent): Expression of the 103 divergent transcripts differentially expressed to a large and significant extent between the sexes in craniosynostosis cases and not in controls.

Supplementary Table 2

Individual qRT-PCR data for CD24 expression in randomly selected cases (n = 24) and controls (n = 8).

Supplementary Table 3

(shared): Expression of 21 shared transcripts differentially expressed to a large and significant extent between the sexes in both controls and craniosynostosis cases.

Supplementary Table 4

(convergent): Expression of 32 convergent transcripts differentially expressed to a large and significant extent between the sexes in controls and not in craniosynostosis cases.

Supplementary Table 5

Individual qRT-PCR data for cases and controls in which FBN2 and IGF2BP3 expression were both in the top third of the entire cohort (n = 249), while maintaining PDE1C expression within the bottom third.