Gene Expression Changes Between Patent and Fused Cranial Sutures in a Nonsyndromic Craniosynostosis Population

Samples

Prior to study enrollment, all patients were examined by a medical geneticist for determination of the status of their craniosynostosis as syndromic or nonsyndromic. Children with only a single suture that was fused, without malformations in other organ systems and molecular genetic testing for FGFR 1, 2, 3 and TWIST mutations failed to reveal any mutation, were considered nonsyndromic. Written informed consent was obtained from the parents of 7 infants diagnosed with NSC prior to initial craniosynostosis surgery. The study protocol was approved by the institutional review board at Virginia Commonwealth University. All sutures were completely fused, which was verified with computed tomographic scans preoperatively and by observation intraoperatively. Bone fragments were recovered from synostosed and patent cranial sutures of each patient undergoing surgical correction. Using a 2-mm Kerrison ronqeur, samples were taken within 1 mm of the location where the ectocranial ridging of the fused suture was most prominent. Patent samples were taken with the same instrument across the patent suture to include no more than 1 mm from the open suture (Supplemental Figure 1S). The fragments were immediately stored in RNA later (Qiagen) solution and then kept at −80°C until ready for isolation. The synostosed fragments of each participant were kept separate from the patent sutural material. Therefore, 14 total samples were included in the study: 7 synostosed and 7 patent samples. The synostosed sutures were derived from 5 sagittal, 1 unilateral coronal and 1 metopic. Participants included male and female infants younger than 1 year (Table 1).

RNA isolation, quality, and integrity

Frozen bone was pulverized and homogenized in TRIzol reagent (Life Technologies). A standard TRIzol protocol was followed to isolate RNA followed by isopropanol precipitation and 70% ethanol wash (https://tools.invitrogen.com/content/sfs/manuals/trizol_reagent.pdf). Purity and concentration of the total RNA were tested using a Nanodrop spectrophotometer. The A260/A280 ratio ranged from 1.45 to 1.95 (average = 1.7). Integrity of the RNA was determined using the Agilent 2100 Bioanalyzer that indicated RNA degradation with RINs ranging from 2.0 to 6.80 (average 4.0). Thus, the Illumina WGE-DASL assay was chosen for this study as it is designed specifically for degraded RNA expression profiling. RNA was not pooled for this study.

Gene expression, preliminary data analysis, and statistical methods

Global gene expression was measured using the Illumina WGE-DASL BeadChip Assay. Reactions were performed according to the manufacturer's WGE-DASL Manual (Illumina). In brief, total RNA was used to generate cDNA using biotinylated oligo-dT18 and random primers. Pooled assay specific oligonucleotides targeting a 50 nucleotide sequence of each cDNA and provided in the WGE-DASL assay pool were annealed to the targeted cDNAs during a 16-hour incubation. This was followed by enzymatic extension and ligation steps, as described in the manufacturer's protocol. Ligated products were PCR amplified and labeled with a universal Cy3-coupled primer after which single-stranded labeled products were precipitated and then hybridized to HumanRef-8 expression BeadChips (Illumina). The bead arrays were washed and coated. Finally, the raw expression intensities, as measured by fluorescence, were captured, by scanning the arrays, on a BeadArray Reader (Illumina).

Raw image files were imported into GenomeStudio (Illumina). All BeadArray data were normalized using the Quantile method without background subtraction. Overall hybridization and array performance were assessed using GenomeStudio QA tools. Normalized annotated sample results were exported in .txt format and imported into Partek Genomics Suite (Build 6.6) for identification of differentially expressed mRNA transcripts.

Changes in mRNA transcript expression were assessed for the synostosed sutures using intensity data of RNA derived from patent suture samples as baseline. Therefore, we compared the gene expression in synostosed sutures to patent suture. Genes were considered to have increased or decreased expression in the prematurely fused suture set relative to the patent suture set. The paired t test was used to determine differentially expressed genes with significance set at ≤.05 and the Benjamini and Hochberg correction for multiple testing was conducted with adjusted P value set at .10. To refine the gene list, fold changes were calculated by using the unlogged intensities to calculate the ratio of change in synostosed suture when using patent suture expression as baseline or control. Genes significant with false discovery rate (FDR) of 0.10 or smaller and fold changes FDR and fold changes (ratios) −2.0 or less and 2.0 or more were retained for further investigation and are referred to as our set of FDR 0.10 genes. Genes were observed approaching significance with FDR of 0.10 and fold changes (ratios) −2.0 or less and 2.0 or more. Therefore, the FDR was relaxed to 0.15 to include these genes. This set, which includes the FDR 0.10 genes, is referred to as our FDR 0.15 genes. We also conducted paired t test on a subset of the samples contributed only from patients with sagittal synostosis. Both males and females were included in the study. Synostosed suture derived from the prematurely fused sagittal suture while the patient's paired patent suture control derived from lambdoidal suture. Total sample size was 10: 5 synostosed and 5 patent. Again, baseline expression was determined in patent sutures with change in expression in synostosed suture relative to this. Analysis was performed using Partek software and the paired t test subroutine. Significance was set at P ≤ .05 with FDR at 0.10 to correct for multiple testing. We also considered a qualitative method of calculating fold changes. Genes significant with FDR 0.10 were then further filtered by implementing a fold change stringency of −2.0 or less and 2.0 or more to the list. We also applied the fold change stringency to the list of significantly differentially expressed genes with relaxed FDR of 0.15. All microarray data have been submitted to the Gene Expression Omnibus database at NCBI with the accession number GSE50796.

Validation

To confirm global gene expression analysis using the Illumina WGE-DASL BeadChip platform, qRT-PCR assay was performed for 4 genes from our gene lists. Reverse transcription of total RNA using the SYBR Green qPCR reagents (Applied Biosystems Inc) was performed. Specifically, first strand cDNA was generated from 2.5 μg of total RNA with 5 ng/μL of oligo dT and random primers and 0.5 mM dNTPs. The reaction was heated to 65°C for 5 minutes and placed on ice. First strand reaction buffer (1X), 0.01 M DTT, and 40 U of RNase inhibitor (Promega) were added and the reaction was incubated at 37°C for 2 minutes. MMLV-RT was added (200 U) (Invitrogen) and the sample was incubated at 37°C for 1 hour followed by a final incubation of 95°C for 5 minutes. qRT-PCR for validation of gene expression changes was performed in triplicate using the Applied Biosystems (ABI) 7500 Fast Real Time System. Primers were designed and optimized using standard PCR. Reactions were set up using 4 μL of cDNA diluted 1:8, 3 pmol/μl primers and 10 μL of 2X SYBR Green qPCR Master Mix (ABI) in a total volume of 20 μL. A 2-step PCR protocol was used with a 10-minute denaturation step at 95°C and then 40 cycles of 95°C for 15 seconds and 60°C for 1 minute. A melt curve analysis was performed to ensure PCR specificity. For each experimental sample, 18S was also run as an internal control. Ribosomal 18S RNA is a common housekeeping gene widely used as an endogenous control including in other studies using human sutures.⁷ Output was analyzed using SDS 7500 Fast software version v1.3.1. Values were obtained using a standard curve of dilutions of cDNA derived from RNA isolated from normal bone and normalized with 18S. The relative standard curve method was used to analyze the expression results and used to calculate fold changes. Table 2 presents the genes tested, primer sequences, and expected amplicon length for the validation experiments (Supplementary Figure 2S).

Functional analysis

FDR 0.10 and 0.15 gene lists and fold changes were evaluated through the Ingenuity Pathway Analysis (IPA) (Ingenuity® Systems, www.ingenuity.com) specifically de novo network analyses and predictions of upstream activators and inhibiters were most useful in analysis of the gene sets. The de novo network analysis used the gene lists derived from the microarray data from these experiments and found known relationships between the genes that created the pathways (Figs 1 and 2). Genes not on the list may be inserted into the pathways to make connections between the genes. Venny was used to create Venn diagrams to identify overlap of genes between studies.^7-9,12 Chi-square analysis with Yates correction was used to determine the significance of overlapping genes (GraphPad).

Figure 1. Pathway analysis. De novo pathway analysis using the false discovery rate (FDR) 0.10 gene list. Pathways were built with Ingenuity Pathway Analysis taking the genes with significant gene expressions changes from this experiment and identifying known relationships between the genes. Green indicates genes downregulated and red indicates upregulated genes from this study. The white genes did not have significant gene expression changes but link the genes in the network that are significantly changed. The intensity of the color indicates the degree of the expression change. Solid lines indicate direct interaction, while dashed lines indicate indirect interaction. Lines indicate binding between molecules, while arrows indicate molecule that can act on another which may or may not also bind.

Figure 2. Pathway analysis. De novo pathway analysis using the false discovery rate (FDR) 0.15 gene list. Pathways were built with Ingenuity Pathway Analysis taking the genes with significant gene expressions changes from this experiment and identifying known relationships between the genes. Green indicates genes downregulated and red indicates upregulated genes from this study. The white genes did not have significant gene expression changes but link the genes in the network that are significantly changed. The intensity of the color indicates the degree of the expression change. Solid lines indicate direct interaction, while dashed lines indicate indirect interaction. Lines indicate binding between molecules, while arrows indicate molecule that can act on another which may or may not also bind.

Correlation between technical replicates

To test the accuracy and precision of the expression data generated from the Illumina BeadArray, we included technical replicates for 5 patent and 5 fused RNAs. RNA from a single isolation for each sample was divided into separate aliquots and processed through the complete Illumina assay, hybridized to a BeadArray, and scanned. Correlations between duplicate sample arrays were all greater than 0.95 (see Supplementary Table 1S). This suggests a high degree of precision providing an additional validation of our results.

Differentially expressed genes

We conducted global gene expression analysis, using RNA from 7 pairs of fused and patent infant cranial suture samples in an effort to identify dysregulated genes. When comparing mRNA transcript expression in prematurely fused human infant cranial sutures to patent sutures, 68 genes were significantly differentially expressed with FDR 0.10 and fold changes ≤ −2.0 and ≥ 2.0; 5 were overexpressed while 63 were underexpressed. With relaxed FDR (0.15), 136 genes were significantly differentially expressed and had fold changes −2.0 or less and 2.0 or more; 21 had increased expression, while 114 showed decreased expression (Table 3).

The top-10 downregulated genes significant at FDR 0.10 and measured by fold change are shown in Table 4. With FDR relaxed to 0.15, SFRP2, GPNMB, SCN2B, LRRC17, EDIL3, and FAP are still included, but DKK2, RORA, RGS22, and FGFBP2 are also in the top 10 significantly underexpressed with fold changes −2.0 or less. The gene lists are not identical since some of the genes added in the more relaxed criteria have greater fold changes than the genes they are displaced in the FDR 0.10 list. Identified genes with increased expression in synostosed cranial suture include DOCK2, CDCA5, CDCA7, MOSC1, and CEACAM1 (Table 5). The top 10 significantly overexpressed genes at FDR 0.15 and having fold changes 2.0 or more include EPB42, CCDC81, KEL, OLR1, CEACAM1, VPREB1, RAB3IP, HBQ1, CEACAM1 (alternative transcript), and TRIM10, a different list than FDR 0.10 except for CEACAM1.

To confirm gene expression changes, quantitative real-time PCR (qRT-PCR) was performed with 4 downregulated genes from our list, SFRP2, GPNMB, MN1, and DKK2 using samples from patients 5 and 6. The amount of RNA from the remaining samples was not sufficient for validation studies and we did not validate these in independent samples due to the difficulty in obtaining additional tissue. SFRP2, GPNMB, and DKK2 were all downregulated comparing patent to fused sutures in both samples (Supplemental Figure 1S and Supplemental Table 2S). DKK2 was chosen to validate as it was identified with the more relaxed criteria; thus, FDR 0.15 genes may still be biologically relevant in craniosynostosis. MN1 was identified in the FDR 0.01 data set but was only 2-fold downregulated according to the microarray. It was downregulated in patient 5 but slightly upregulated in patient 6. Additional samples would be necessary to determine whether MN1 is generally downregulated as sutures fuse.

Correlation of expression between patients

Even though our samples were derived from only NSC individuals, we used different patent and fused suture type combinations in this study as well as samples from both male and female individuals. To determine whether this combination of different suture type tissue is an issue, we used hierarchical clustering with our significantly differentially expressed genes across samples (Fig 3). We found that in all but 1 case, the expression of the fused sutures cluster together as well as the patent suture irrespective of suture type or gender. In general, neither suture type nor sex clustered together within a group (patent vs fused). The 1 exception was individual 7, whose fused suture clustered with the patent suture group most closely with patent suture from patient 4, the only sagittal patent suture studied.

Figure 3. Heat map and hierarchical clustering. Clustering of individual suture samples using significantly differentially expressed genes. Individual suture samples are listed on the left with patient numbers (Table 1) and P referring to patent sutures and F to fused sutures.

To further identify the molecular signature of NSC and genes most important in maintaining patency across suture type and sex, we refined our gene list to the top-10 significantly upregulated and downregulated transcript probe sets for the entire sample. Hierarchical clustering as conducted using log 2 transformed normalized signal for only these probe sets to create a heat map and dendrogram (Fig 4). Both samples and genes were clustered. Measure of dissimilarity was Pearson's dissimilarity with complete linkage. Using only the top-10 up- and downregulated genes 2 major branches separated fused and patent samples. Within suture status, sex and suture type did not influence clustering. Suture expression represented by log 2 transformed signal converted to z score with mean of 0 and SD of 1 clearly delineates patent from fused. Thus, our data show that we can identify genes that are generally changed between patent and fused sutures, irrespective of sex or suture type.

Figure 4. Heat map and hierarchical clustering. Clustering of individual suture samples and top 10 genes significantly upregulated and downregulated genes. Individual suture samples are listed on the left with patient numbers (Table 1) and P referring to patent sutures and F to fused sutures.

Differentially expressed genes in sagittal sutures only

Results of paired t test considering sagittal suture synostosis only resulted in 2 significantly differentially expressed genes when implementing FDR 0.10 with P ≥ .05 and fold change stringency −2.0 or less and 2.0 or more (Table 2S). With relaxed FDR to 0.15, 119 genes were significant and had fold changes −2.0 or less and 2.0 or more; 114 were downregulated in prematurely fused sagittal suture while 6 demonstrated increased expression (Supplemental Table 3S). SFRP2 and SPHKAP were downregulated and significant at FDR 0.10. Interestingly, SPHKAP is not identified in the original analysis, perhaps indicating a gene important in sagittal suture fusion. Considering genes significant with FDR 0.15, there were only 6 upregulated genes, which include APOBEC3F, AMMECR1, IKZF3, MYO3A, OR10H5, and IRE4. The top-10 downregulated genes were B4GALT4, C10ORF26, CCD074A, EPHB2, GPC3, LAMA2, RORA, and SFRP2.

Comparing the sagittal only gene set with the entire set reveals that SFRP2 is significantly downregulated in prematurely fused suture across the entire set at FDR 0.10 (fold changes ≤ −2.0). This suggests that SFRP2 has a significant role as a universal mechanism in maintaining suture patency. Considering genes differentially expressed with the relaxed FDR, 35 genes are found in common. These genes likely represent a set of genes important in suture patency across all sutures.

Gene functions

To determine whether the gene expression data have biological meaning, the significant genes in both FDR 0.10 and 0.15 sets were analyzed through the Spring 2013 release of Ingenuity Pathway Analysis including Upstream Analysis, which identifies upstream regulators that can explain the observed gene expression changes (Ingenuity^® Systems, www.ingenuity.com). To test if it is valid to include genes with less stringent criteria (FDR 0.15), we see that the significant canonical pathways overlap in the 2 groups of genes (Table 6). Both gene lists (FDR 0.10 and FDR 0.15) identified the same canonical pathway (osteoblasts, osteoclasts, and chondrocytes in rheumatoid arthritis). In fact, it is the only significant canonical pathway at FDR 0.15, which justifies including this list of genes. The identification of these specific pathways would be expected in genes important in suture biology.

The top de novo networks generated with the 2 different data sets connect to interesting cellular signaling genes (Figs 1 and 2). Genes in both pathways include ERK1/2 and PDGFBB. The pathway created with FDR 0.10 gene list also include p38 MAPK, Akt, and IRS1, while the top network with FDR 0.15 genes also includes TGFβ, all genes important in cellular morphology, cell assembly and organization, connective tissue, and known to be important in suture biology.

Finally, we looked at putative upstream regulators of these differentially expressed genes (Table 7). Interestingly, the most significant molecule and the ability to predict how the molecule might be regulated (z scores) in both gene sets is the inhibition of WNT3A. The only other change found in both lists is the inhibition of TGFβ1. It seems that in NCS Wnt and TGFβ signaling might be involved in suture patency and FGFs may not be the activator upstream of suture closure in NS cases.

Genes across studies

To determine whether our gene list corroborates previous studies, we compared the gene list from this study with genes identified in 3 previous human studies (Fig 5). First, we compared our gene list with study results of differentially expressed genes between fused sutures and suture tissue that were dedifferentiated in culture.⁷ There was no overlap with our most relaxed gene list (FDR 0.15) with the top-50 genes from that study. A large number of differentially expressed genes were identified in that study (∼2000) so the genes identified from this study could have overlapped with the genes in the larger list. The second list was of single sutures with tissue also grown in tissue culture and compared to a set of unrelated controls.⁹ Again, only the top genes from that study were analyzed (47 genes) and 2 genes overlap with our FDR 0.15 gene list, GEM, MAB21L2 (P = .007) (Fig 5). Interestingly, comparing the gene list from this study with a much smaller list of significantly differentially expressed genes (28 genes) using primary suture tissue, 6 genes in common were identified with the FDR 0.10 list and include FBLN1, FMOD, MFAP4, MN1, RBP4, SPON1 (P < .0001) (data not shown).⁸ One additional gene can be included using a more relaxed FDR 0.15 (PDZRN3) (P < .0001) (Fig 5). Our study identifies 25% of their significantly changed genes, which supports the importance of these genes in suture patency and that these downstream changes seem to be irrespective of the primary cause of craniosynostosis in these cases (Fig 5). There was no overlap between gene lists between any of the previously published studies.

Figure 5. Venn diagram. Overlap of the number of genes identified as significantly differentially expressed (FDR 0.15) in this study and other studies.7-9

Wnt signaling

Wnt signaling plays an important role in bone development and it is well known that an increase in Wnt signaling increases bone mineralization (reviewed in Baron and Kneissel¹⁹). In contrast to mature bone mineralization, when specifically investigating Wnt signaling in suture closure in mice, continual Wnt signaling is necessary for suture patency in sagittal sutures through the downregulation of TWIST.²⁰Wnt3a inhibition is predicted by Ingenuity Pathway Analysis to be upstream of gene expression changes from this study and seems to reinforce the findings of Behr et al.²⁰ Specifically, Wnt3a has been shown to inhibit chondrogenesis through the reorganization of the cytoskeleton and prevent osteogenic differentiation by promoting proliferation.^21-23

SFRP2, a secreted frizzled-related protein, was identified as an important gene downregulated in fused sutures as it is significantly decreased in fused tissue in all data sets (0.10 FDR, 0.15 FDR, and sagittal only 0.10 FDR and 0.15 FDR). Another Wnt/β-catenin pathway gene is also significantly downregulated, DKK2 (top 10 in FDR 0.15). DKK2 is involved in terminal osteoblast differentiation. Both SFRP2 and DKK2 were validated with qRT-PCR and have been identified as negative regulators of canonical and noncanonical Wnt/β catenin signaling pathway. One would expect that increased canonical Wnt signaling would result if negative regulators of Wnt signaling were downregulated. It is likely to be more complicated. DKK2 can be an agonist or antagonist to Wnt signaling. DKK2 can activate Wnt signaling through the LRP6 receptor.²⁴ In the vertebrate optic cup, SFRP1 and 2 was shown to positively modulate canonical Wnt signaling by helping diffuse Wnt through the tissue and bringing it into close proximity of its receptor.²⁵ Thus, downregulation of SFRP2 and DKK2 does not necessarily indicate an increase in Wnt signaling since SFRP2 and DKK2 can either antagonize or positively modulate Wnt signaling.

Additional differentially expressed genes and their involvement in craniosynostosis

Other top-10 downregulated genes include GPNMB, LRRC17, EDIL3, FAP, LRIG3, and GDF10 (FDR 0.15). These genes are known to be involved in osteogenesis, cell adhesion, and invasion or angiogenesis. For example, LRRC1 acts as a negative regulator of RANKL-induced (receptor activator of NF-κB ligand) osteoclast differentiation.²⁶ Blood vessel formation is known to be an important aspect of bone formation. The 4 additional downregulated genes in the top-10 list (SCN2B, FOXP2, PCDH19, and RBP4) all are found to be expressed in the brain or have brain phenotypes. Interestingly, signals from the dura and brain development are also thought to affect suture biology.

The top upregulated genes in the FDR 0.10 and FDR 0.15 lists are not known to be important in bone formation or suture biology. They are involved in the endothelial cells, hematological tissues, and important in cancer. These are candidate genes that may be found to be important in suture biology.

Sagittal suture synostosis

SPHKAP was identified as a significantly downregulated gene in the analysis of only fused sagittal sutures that may indicate the importance of this gene in preserving sagittal suture patency. This gene codes for a type I A-kinase anchoring protein and has been shown to interact with and modulate sphingosine kinase activity (SPK1).²⁷SPK1 mediates TGFβ signaling of extracellular matrix deposition and maintenance.²⁸ The extracellular matrix contains components that can signal through receptors that control proliferation, migration, differentiation, and apoptosis. In a recent genome-wide association study (GWAS) using trios with children with sagittal craniosynostosis, 2 regions of the genome that were associated with this phenotype were BMP2 and BBS9.²⁹BMP2 is a TGFβ family member. Although no known connection between SPHKAP and BMP2 has been observed, it will be interesting to determine whether SPHKAP has any effect on the signaling of BMP2 specifically in sagittal sutures.

Genes identified across studies

There have been a variety of studies that have identified differentially expressed genes between patent and fused sutures. Most studies included tissue from individuals with known syndromic mutations in their study population; others cultured the osteoblasts prior to gene expression analysis. In addition, different experimental platforms and different criteria were used to generate gene lists. This makes it difficult to compare gene lists across studies. We looked at gene lists from 3 additional human studies^7-9 and there was no overlap using their short gene lists between the previous studies but we had significant overlap with 7 genes with the study that also used primary suture tissue as their starting material⁸ and 2 genes with a study with cultured osteoblast. These 9 genes, GEM, MAB21L2, FBLN1, FMOD, MFAP4, MN1, RBP4, SPON1, and PDZRN3, are likely very important across suture type and irrespective of the cause of craniosynostosis.

FBNL1, FMOD, MFAP4, and SPONL are all proteins that are found in the extracellular matrix. The extracellular matrix not only is an exceedingly important component of bone formation as a structural component but also plays a role in cell signaling. For example, FMOD is a small proteoglycan that can bind collagen and inhibit fibrillogenesis as well as has the ability to sequester grow factors such as TGFβ in the extracellular matrix.^30,31 Two of the genes, MN1 and PDZRN3, are known to be important in osteoblast induction or function. Gene ontology revealed that MN1 has a critical function in intramembranous ossification. MN1 is the gene that was validated only in 1 out of 2 samples by qRT-PCR that would not preclude its importance in craniosynostosis in some cases. MN1 is a transcription coregulator and studies show that MN1 is necessary to maintain appropriate osteoblast proliferation, differentiation, and function, including support of osteoclastogenesis.³²PDZRN3 is necessary for the negative feedback loop of BMP2-induced osteoblast differentiation. This feedback loop utilizes the canonical Wnt signaling pathway.³³

Signaling networks and activators predicted to be involved in craniosynostosis

To identify upstream activators or other genes related to our genes lists, IPA analysis was used to create de novo networks and to identify putative upstream activators (Figs 3 and 4 and Table 5). Signaling molecules such as p38/MAPK, PI3K/Akt, ERK1/2, PDGF-BB and TGFβ were included with genes from our gene list to create our top gene networks. Constitutively active FGFRs are known to cause craniosynostosis and FGFRs also signal through p38/MAPK, PI3K/Akt and ERK1/2. Our samples have no known FGFR mutations, but similar downstream pathways seem involved in this craniosynostosis cohort. We do not know if there are changes in any FGF molecules but other upstream signals can affect these molecules. Activation of PI3K/Akt and ERK1/2 can be achieved through TGFβ. PDGF-BB can also signal through p38/MAPK, PI3K/Akt, and ERK1/2. TGFβ1 and 3 and Wnt3a were identified as putative upstream molecules (Table 5). Both are known to be important in bone formation.