Whole-genome methylation analysis reveals novel epigenetic perturbations of congenital scoliosis.

Congenital scoliosis (CS) is a congenital disease caused by malformations of vertebrae. Recent studies demonstrated that DNA modification could contribute to the pathogenesis of disease. This study aims to identify epigenetic perturbations that may contribute to the pathogenesis of CS. Four CS patients with hemivertebra were enrolled and underwent spine correction operations. DNA was extracted from the hemivertebrae and spinal process collected from the specimen during the hemivertebra resection. Genome-wide DNA methylation profiling was examined at base-pair resolution using whole-genome bisulfite sequencing (WGBS). We identified 343 genes with hyper-differentially methylated regions (DMRs) and 222 genes with hypo-DMRs, respectively. These genes were enriched in the mitogen-activated protein kinase (MAPK) signaling pathway, calcium signaling pathway, and axon guidance in Kyoto Encyclopedia of Genes and Genomes (KEGG) pathways and were enriched in positive regulation of cell morphogenesis involved in differentiation, regulation of cell morphogenesis involved in differentiation, and regulation of neuron projection development in Biological Process of Gene Ontology (GO-BP) terms. Hyper-DMR-related genes, including IGHG1, IGHM, IGHG3, RNF213, and GSE1, and hypo DMR-related genes, including SORCS2, COL5A1, GRID1, RGS3, and ROBO2, may contribute to the pathogenesis of hemivertebra. The aberrant DNA methylation may be associated with the formation of hemivertebra and congenital scoliosis.

Introduction

Congenital scoliosis (CS) is a form of spinal deformity affecting 0.05%–0.1% of newborns., In the clinic, CS can be classified into three types based on the causes, including failure of formation, failure of segmentation, and mixed deformity. Hemivertebra is the most common consequence of failure of formation, which is characterized by the absence of one side of the vertebral body and one vertebral pedicle.

Genetic defects have been reported to be associated with congenital scoliosis.4, 5, 6, 7 Besides, environmental factors, including hypoxia and high altitude, could also increase the risk for CS. Although studies exploring the pathogenesis of CS have continuously boosted for decades, the disease mechanism of the majority of CS is still unclear.

DNA methylation is an epigenetic modification that usually gives rise to 5-methylcytosine (5-mC) by targeting the fifth carbon of the pyrimidine ring of cytosine. In mammals, transcriptional silencing is the most pivotal function of DNA methylation. It is postulated that the aberrant DNA methylation plays an important role in congenital diseases, such as congenital renal agenesis and congenital heart disease. Recent studies also indicated that differential methylation of key genes or the CpG site was related to scoliosis.12, 13, 14 However, there is no known methylation region associated with CS.

In this study, we enrolled four CS patients with hemivertebra. We compared their genome methylation difference between the hemivertebra body and spinal process, and explored the epigenetic perturbation to the pathogenesis of CS.

Results

Clinical information

Two of the CS patients were female, and the mean age of patients was 7.75 (from 2 to 14) years old (Table 1). Three of them were classified into type 1 CS, and one was classified into type 3 CS. The average angle of main Cobb angle was 66 degrees. All of the patients had hemivertebra (Figure 1).

Table 1

Characteristics of the enrolled CS patients

ID	Gender	Age (y)	CS classification	Main Cobb angle	Hemivertebra location
DISCO552	F	12	1	73	T4
DISCO644	M	14	3	75	T7
DISCO799	F	3	1	36	L3
DISCO831	M	2	1	80	L1

CS, congenital scoliosis; F, female; M, male.

Figure 1

Spine images of the four CS patients with hemivertebra

DNA methylation differences

The distribution of global methylation level (Figure 2A) and sequencing depth (Figure 2B) between the hemivertebra and spinal process in the four CS patients were similar. We compared the methylation level distribution of the whole-genome DMRs between the hemivertebra body and the spinal process in each patient (Figure 2C). To characterize the distribution of DMRs, we assigned them to known genomic features, which include promoter, 5′ UTR, 3′ UTR, and intronic and intergenic regions. Among all of the hyper-DMRs, 36% were located in intergenic regions, whereas 14% were located in promoter regions. Similarly, for hypo-DMRs, 40% were located in intergenic regions, and a slightly lower percent (9%) was located in promoter regions (Figure 3). Finally, we identified 343 genes with hyper-differentially methylated regions (DMRs) and 222 genes with hypo-DMRs that are consistent among the four patients (Figure 4). We performed an unsupervised hierarchical clustering analysis using DNA methylation level of 343 hyper-DMR-related genes and 222 hypo-DMR-related genes to characterize the global methylation alterations between the hemivertebra group and spinal process group (Figure 5). We enriched these DMR-related genes and then subjected them to Kyoto Encyclopedia of Genes and Genomes (KEGG) and Biological Process of Gene Ontology (GO-BP) terms for pathway analysis. It revealed that DMR-related genes were enriched in several pathways. Three pathways with the most significant association were mitogen-activated protein kinase (MAPK) signaling pathway, calcium signaling pathway, and axon guidance in KEGG pathways. The top three pathways with the most significant association in GO-BP terms were positive regulation of cell morphogenesis involved in differentiation, regulation of cell morphogenesis involved in differentiation, and regulation of neuron projection development (Figure 6), indicating that these pathways could be involved in the development of hemivertebra in CS patients.

Figure 2

DNA methylation differences in different samples

(A) Violin plot of the methylation level distribution in different samples. (B) Sequencing depth in different samples. (C) Boxplot of DMR methylation level.

Figure 3

Characterization of identified DMRs between hemivertebra group and spinal process group

(A) Co-localization of hyper-DMRs with known genomic features. (B) Co-localization of hypo-DMRs with known genomic features.

Figure 4

Matrix layout for visualization of all intersections of DMR-related genes between four patients, sorted by degree and size

Dark circles in the matrix indicate sets that are part of the intersection. (A) All intersections of hyper-DMR-related genes. (B) All intersections of hypo-DMR-related genes.

Figure 5

The unsupervised hierarchical clustering heatmap of hemivertebra and spinal process in four patients based on DNA methylation patterns of 343 hyper-DMR-related genes and 222 hypo-DMR-related genes

Figure 6

Functional enrichment of DMR-related genes

(A) Bar chart of top 10 DMR-related genes enrichment with GO-BP. (B) Scatterplot of top 10 DMR-related genes enrichment with KEGG pathway.

Validation of candidate DMR-related genes

To identify the specific pathogenic epigenetic variants, we ranked the recurrent DMR-related genes according to the arithmetic mean value of the methylation difference between hemivertebra body and normal spinal process in each group (Table S1). The top five candidate DMR-related genes with greatest hyper-methylation differences were IGHG1, IGHM, IGHG3, RNF213, and GSE1 (Table 2). The top five candidate DMR-related genes with greatest hypo-methylation differences were SORCS2, COL5A1, GRID1, RGS3, and ROBO2 (Table 2). We hypothesize that the different methylation of DMR-related genes may influence the gene expression by regulating the transcription and be associated with the phenotype of hemivertebra of CS.

Table 2

Top 5 common DMR-related genes with greatest mean hyper- and hypo-methylation differences in the CS patients

DMR-related genes	Gene ID	DNA methylation direction in HV	Mean methylation difference
IGHG1	ENSG00000211896	hyper	0.39
IGHM	ENSG00000211899	hyper	0.33
IGHG3	ENSG00000211897	hyper	0.33
RNF213	ENSG00000173821	hyper	0.28
GSE1	ENSG00000131149	hyper	0.24
SORCS2	ENSG00000184985	hypo	−0.25
COL5A1	ENSG00000130635	hypo	−0.23
GRID1	ENSG00000182771	hypo	−0.21
RGS3	ENSG00000138835	hypo	−0.21
ROBO2	ENSG00000185008	hypo	−0.20

SORCS2 encodes Sortilin Related VPS10 Domain Containing Receptor 2 (SORCS2), which is one family member of the vacuolar protein sorting 10 (VPS10) domain-containing receptor proteins. It is indicated that SorCS2 was prominently expressed in the vertebrae of mice embryo. In a genome-wide association study (GWAS), a locus of SORCS2 is significantly associated with insulin-like growth factor-binding protein-3 (IGFBP-3), which is involved in bone metabolism. We hypothesized that the abnormal methylation could influence the expression of SORCS2, which could induce the malformation of vertebrae.

COL5A1 encodes an alpha chain for one of the low-abundance fibrillar collagens. Mutations of this gene could lead to Ehlers-Danlos syndrome (EDS). EDS is a kind of connective tissue disorder, characterized by severe muscle hypotonia at birth, joint hypermobility, hyperelastic skin, osteopenia, increased bone fragility, and progressive scoliosis., Although the patients in our study manifested with only hemivertebra, the abnormal methylation of COL5A1 could lead to abnormal expression of collagen, which may increase the risk of scoliosis development.

RGS3 encodes Regulator of G-protein signaling 3, which is an antagonism of G protein and interacts with ephrin. It is postulated that the dislocation of RGS3 could impair the function of Ephrin-B-RGS cell fate signaling complex. The ephrin signaling within the osteoclast and osteoblast lineages could promote or inhibit cell differentiation in bone remodeling. Different methylation level of RGS3 could influence bone remodeling and have an effect on the formation of hemivertebra.

ROBO2 encodes Roundabout Guidance Receptor 2, which belongs to the ROBO family and functions in axon guidance and cell migration. The protein could function as receptor for SLIT2 and regulate bone morphogenetic proteins (BMPs) activity. Previous studies indicate that ROBO2 could promote chondrocyte maturation. In a rat animal model, mRNA of robo2 was expressed during the differentiation of osteoblast. The aberrant methylation of ROBO2 could regulate the differentiation of osteoblast and chondrocyte maturation, which may contribute to the manifestation of hemivertebra. Further functional studies are needed to explore the mechanism.

Discussion

In general, genetic variants and epigenetic modifications are pivotal causes to embryonic somatogenesis failure., Environmental factors, including hypoxia, are validated principal regulators of fibroblast growth factor (FGF) signaling, which is an essential pathway in the pathogenesis of CS. Up until now, the epigenetic modification of CS is still unclear. In this study, we present the epigenetic landscape of CS, or more specifically hemivertebra. According to previous studies, DNA methylation can predispose to a wide range of congenital developmental diseases., Some studies also indicated that DNA methylation is related to scoliosis.12, 13, 14 However, there is no study investigating the relation between the whole-genome DNA methylation pattern and CS.

Therefore, we performed whole-genome bisulfite sequencing (WGBS) to detect the whole-genome DNA methylation between the hemivertebra body and normal spinal process of CS patients. Potential pathological signal pathways and DMR-related genes were identified. DMR-related genes with differently methylated level were enriched in several signaling pathways, including the MAPK signaling pathway, calcium signaling pathway, and axon guidance in KEGG pathways, and positive regulation of cell morphogenesis involved in differentiation, regulation of cell morphogenesis involved in differentiation, and regulation of neuron projection development in GO-BP terms.

The MAPK signaling pathway is a critical pathway regulating the proliferation and differentiation of multiple cell types, including osteoprogenitor cells. Mutations of the genes involved in the MAPK signaling pathway are associated with multiple diseases manifesting scoliosis, such as Noonan syndrome, Costello syndrome, and cardiofaciocutaneous (CFC) syndrome. Transcriptional regulation, including miRNA and DNA methylation in genes involved in the MAPK signaling pathway, was also associated with scoliosis., Thus, we hypothesized that the abnormal DNA methylation of the MAPK signaling pathway-related genes could lead to the hemivertebra development.

Calcium metabolism is very important in bone formation., The abnormal calcium signaling pathway could lead to osteoblast dysfunction and scoliosis., Previous studies also proved that the abnormal transport of calcium ion is involved in scoliosis., Previous studies indicated that the proteins of axon guidance could guide growing axons during development. The defective axon guidance could lead to several diseases manifesting scoliosis., Thus, the abnormalities of these pathways could have potential to increase the risk for hemivertebra.

We also found several common DMR-related genes manifested with the same methylation tendency. Hyper-DMR-related genes, including IGHG1, IGHM, IGHG3, RNF213, and GSE1, and hypo-DMR-related genes, including SORCS2, COL5A1, GRID1, RGS3, and ROBO2, were the common DMR-related genes with greatest mean methylation differences.

There are several limitations in this study. First, the abnormal methylation can cause transcriptome disorder in target tissues. In our study, the transcriptome data were not detected; further studies are needed to explore the relationship between abnormal methylation modifications and unfaithful activation of key genes. Second, although this is the pilot, and to our knowledge the first study elucidating the relation between epigenetic modification and congenital scoliosis, the sample size of the study was relatively small. Further studies with larger sample size are necessary to validate our findings.

Materials and methods

Patients and materials

Four patients with CS were recruited from Peking Union Medical College Hospital (PUMCH) from July 2015 to January 2016, as a part of the Deciphering Disorders Involving Scoliosis and Comorbidities (DISCO) study (http://www.discostudy.org/). All the patients had no relationship to each other. All the CS patients manifested with hemivertebra and underwent hemivertebra resection. Parts of the hemivertebra body and normal spinal process were collected from the specimen during the surgery. Genomic DNA was extracted from those specimens using DNeasy Blood & Tissue Kits (QIAGEN, Eastwin Scientific, Beijing, China) according to the manufacturer’s instructions. Written informed consent was obtained from all the participants or their parents. The Ethical Review Board of Peking Union Medical College Hospital approved this study.

WGBS and data analysis

After the extraction of genomic DNA, WGBS was performed and data analysis was conducted as described in our previous study. DMRs between the hemivertebra body and the spinal process were identified using MethPipe (v.3.4.3) with significant differentially CpGs ≥5 and a minimal number of 10 CpGs that the DMR spans.

According to the methylation level difference, all the DMRs were classified into two groups, the hypermethylation group and hypomethylation group. ChIPseeker (v.1.14.0) Bioconductor package was used to evaluate the genomic position distribution of DMRs. Every DMR was assigned to a related RefSeq gene from 3 kb downstream to 3 kb upstream of the TSS. The DMR-related genes were defined as genes with at least one DMR mapped to it. In the case of multiple DMRs being mapped to the same gene, we retained the DMR with the largest absolute methylation difference to represent the gene. Hypo/hyper-DMR-related genes observed among four patients were identified as recurrent DMR-related genes. Hierarchical clustering analysis was performed to characterize the global methylation alterations between the hemivertebra group and spinal process group using the R package. Then, the recurrent DMR-related genes were ordered by the arithmetic mean value of the methylation difference of corresponding DMRs in four groups. All of the common DMR-related genes were enriched at Gene Ontology (GO) and KEGG using clusterProfiler (v.3.6.0) Bioconductor package.