A pair of twins (females, 15 years old) were enrolled into the present study from the Department of Orthopedic Surgery, The People's Liberation Army General Hospital in June 2017. The twin-patient (T-P) individual suffered from CS, and the twin-healthy (T-H) individual did not (Fig. S1). The inclusion criteria were: i) Spinal scoliosis diagnosed by radiological imaging and clinical symptoms; ii) T1-L5 vertebral formation defects including hemivertebra, butterfly vertebra and wedge vertebra; iii) T1-L5 vertebral defects of segmentation including vertebral bar and block vertebrae; and iv) with or without rib abnormalities. The exclusion criteria were: i) Presence of frequently encountered syndromes associated with congenital vertebral malformation (CVM), including Alagille syndrome, Goldenhar's syndrome, Jarcho-Levin syndrome, Klippel-Feil syndrome, Sotos syndrome and VACTERL association conditions; ii) vertebral malformations in association with other renal, cardiac or spinal cord malformations; and iii) cervical CVM. Informed consent was obtained from the patients and their parents. The study (approval no. 3012017114) was approved by the Ethics Committee of The People's Liberation Army General Hospital and written informed consent was obtained from all participants for their DNA samples to be used in the experimental procedures and their images to be published. Any personally identifiable information of the participants were removed. The study was conducted in accordance with the Declaration of Helsinki.

DNA was extracted from peripheral blood of monozygotic twins who were discordant for CS. To confirm that both twins were monozygotic, microsatellite analysis (Excel Microsatellite Toolkit; Version 3.1.1., http://courses.washington.edu/ fish543/Software.htm) was performed as previously described (15). The common SNPs were filtered out from the total SNPs using the SNPinProbe database (dbSNP v129, https://skydrive.live.com/?uc=2&cid=fb2a64e 541add2be) and 1,000 Genomes Project (16).

The genomic DNA samples that met all the inclusion criteria (from the twins) were randomly fragmented and the fragments of 150–200 bp length were purified as described in a previous study (17). Adaptors were ligated to both ends of the fragments, followed by ligation-mediated polymerase chain reaction (PCR), purification, and hybridization using a SureSelect Biotinylated RNA Library with baits (Thermo Fisher Scientific, Inc.). The DNA fragments unlinked to streptavidin beads (Thermo Fisher Scientific, Inc.) were washed out following 24 h hybridization. Each captured library was sequenced using the Illumina Hiseq 2000 sequencing system (Illumina, Inc.).

Raw image files were processed by Illumina Basecaller Software 1.7 for base-calling with default parameters (Illumina, Inc.). The exome sequences of each individual were generated as 90-bp paired-end reads. The map of the reads onto the Genome Reference Consortium Human Build 37 (GRCh37; University of California Santa Cruz Genome Reference Consortium; http://genome.ucsc.edu/cite.html) was performed using SOAPaligner/SOAP2 (Release 2.21, http://soap.genomics.org.cn/down/soap2.21release.tar.gz), with 3 mismatches allowed per read. The minimum and maximum insert sizes allowed were 90 and 600 bp, respectively. Coverage and depth calculation were based on all mappable reads in the exome regions. The mean, 2.1 Gb of mappable sequences were generated per individual and ~60% of the reads were mapped to the exome. Database for Annotation, Visualization and Integrated Discovery (DAVID) Bioinformatics Resources 6.8 (https://david.ncifcrf.gov/), with P<0.01 was determined to confer statistical significance.

Genomic DNA was fragmented to 100–500 bp by sonication (frequency, 18–20 KHz). Adaptors were ligated to both ends of the fragments. The end repair, base addition and adaptor ligation steps were performed using the Illumina Paired-End DNA Sample Prep kit (Illumina, Inc.) according to the manufacturer's protocol. The adaptor-ligated DNA was immunoprecipitated with anti-5mC antibody (Abcam; cat. no. ab73938; 5MC-CD, 1:200) at room temperature for 1 h, as described previously (18), captured by protein A/G agarose beads (Applied Biosystems; Thermo Fisher Scientific, Inc.), followed by the MeDIP products were validated by qPCR. qPCR was performed using an ABI7500 qPCR instrument (Applied Biosystems; Thermo Fisher Scientific, Inc.). In total, 5 µl SYBR Premix Ex Taq II (Takara Biotechnology, Co., Ltd.), 0.4 µl forward primer (10 µM), 0.4 µl reverse primer (10 µM), 0.2 µl ROX Reference Dye (Takara Biotechnology, Co., Ltd.), 1.0 µl cDNA template and 3.0 µl ddH₂O were mixed in the reaction solution. qPCR was performed as follows: Initial denaturation for 10 sec at 95°C, followed by 45 cycles of 5 sec at 95°C and 34 sec at 60°C. Relative gene expression was calculated using the comparative quantification cycle (Cq) method. Fold changes were calculated using the 2^−ΔΔCq method (19). Primers were designed as follows: RAB3C forward, 5′-GATTCCACGCTTTCCCTCCAG-3′ and reverse, 5′-AAATTCTCAAGTGCCCTCGCTACA-3′; and β-actin forward, 5′-GATCATTGCTCCTCCTGAGC-3′ and reverse, 5′-ACTCCTGCTTGCTGATCCAC-3′. For the selection of the genomic DNA libraries that were suitable for the Illumina Genome Analyser II, DNA was amplified, purified and size-selected (200–300 bp including the adaptor sequence) using an agarose gel extraction kit (Tiangen Biotech Co., Ltd.), according to the manufacturer's instructions. The sequences of each individual were generated as 44-bp paired-end reads.

The sequencing results were mapped to the human reference genome GRCh37. The mean, 4.4 Gb of mappable sequences were generated per individual and ~89% of the reads could be uniquely mapped to the reference genome. The peak calling of the mappable reads was performed using the MACS software (version 1.4, http://miltenyi.ebiomall.com/18-1.html). The MACS parameters for peak calling of each individual were-g hs-nolambda-nomodel-shiftsize 150, whereas the parameters for differential peak calling were-gsize hs-nomodel-shiftsize 150. In addition, the parameters for T-H specific peak calling were-t T-H.bed-c T-P.bed, and for T-P specific peak calling-t T-P.bed-c T-H.bed. As the sequencing depth of the MeDIP-seq experiment was high (~95 million reads per individual), P<10⁻²⁰ was used as the cutoff for peak calling. The allelic bias in the MeDIP-seq read coverage for each individual was determined by χ² test (P<0.01). SPSS software (version 18.0; IBM Corp.) was used for statistical analysis.

To determine the DNA sequence variation in the disease-discordant twins, WES was performed using the Illumina Hiseq 2000 sequencing system. A total of ~2.6 (with 78.4% mappable reads) and 2.5 Gb (83.6% mappable reads) sequences for T-H and T-P were generated, respectively, of which 60% were annotated as exons. The mean sequencing depth of the target region of the exome was 34.09- and 36.01-fold, respectively, which was considered sufficient for accurate SNP detection. Subsequently, 25,384 SNPs distinct from the reference human genome GRCh37, but shared by the twins, were identified using SOAPsnp (Release 1.03, http://soap.genomics.org.cn/soapsnp.html). Among them, 1,931 and 3,144 SNPs were T-H and T-P specific, respectively (Table SI). In addition, 1,194 (shared), 188 (T-H specific) and 319 (T-P specific) SNPs were located in CDS and annotated as non-synonymous mutations (Fig. 1).

Although common variants, defined as those present in >5% of the population according to the 1,000 Genomes Project (http://www.1000genomes.org/page.php) (16), may be involved in common diseases or traits, rare variants are considered more likely to be involved in rare diseases (20). Therefore, in the present study, minor allele frequency variants were identified (≤1%). In brief, the common SNPs were filtered out from the total SNPs using the SNPinProbe database (dbSNP v129, http://skydrive.live.com/?uc=2&cid=fb2a64e541add2be) and 1,000 Genomes Project (16). The analysis revealed 1,194 shared rare SNPs in the twins, 188 specific to the T-H and 319 specific to the T-P. Among these rare variants, 665 (shared), 87 (T-H specific) and 143 (T-P specific) were located in CDS regions, of which 367, 60 and 86, respectively, were non-synonymous (Fig. 1 and Table SI).

Subsequently, functional analysis of the genes with non-synonymous changes in the CDS was performed. Given that the basic traits of CS are developed during gestation and that cytoskeletal mechanisms are considered important for proper spinal development, the present study specifically focused on these categories. However, although >3,000 genes with non-synonymous CDS changes were identified using the total number of SNPs, none of these were demonstrated to be enriched during the embryonic development or cytoskeleton classifications. In contrast, when candidate genes (n=329) with rare non-synonymous CDS changes were considered, the cytoskeleton and development functional terms were significantly enriched (P<0.00032; Table I). Taken together, these results supported the hypothesis that rare mutations were more likely to induce the development of rare diseases. Furthermore, the T-P-specific mutations, according to functional enrichment terms, appeared to be more likely to affect CS development than T-H-specific mutations.

Several studies have proposed that the environment may contribute to CS pathogenesis by effecting DNA methylation (5,7,21,22). Therefore, the methylation profiles of the twins were investigated using MeDIP-seq. Following fragmentation, immunoprecipitation and high-throuhgput sequencing, ~95 million reads (2.6 Gb) per individual were generated, with 89% of the reads mapped to the genome. Overall, the DNA methylation profiles of T-H and T-P were identical at both the genomic level and in gene-associated regions (promoter and gene body), indicating that the twins did not exhibit any significant epigentic differences (Fig. S2). As the sequencing depth was sufficiently high, ~95 million reads per individual, stringent criteria were employed in order the MeDIP-seq peaks to be indetified (P<10⁻²⁰) by using MACS (23). The analysis identified 132,130 and 125,907 peaks for T-H and T-P, respectively. In addition, the same peak calling criteria were used in order the individual-specific peaks to be identified. The differential DNA methylation regions in the twins were 0 and 7 in T-H and T-P specific peaks, respectively. Among the 7 T-P specific peaks, 3 were aligned to RAB3C, member RAS oncogene family, ATPase phospholipid transporting 10B (putative) and leucine zipper protein 2 genes. Although the possibility that methylation changes in ≥1 of these genes may be associated with the pathogenicity of CS cannot be excluded, the small number of the differential regions indicated that this probability is low.

Previous studies in mice and humans have identified 57 genes that may be involved in CS pathogenesis, the majority of which have been demonstrated to be crucial for normal somite formation (6,10,24–32). The known CS candidate genes are summarized in Table SII. Therefore, analyses on these genes were performed so that differences in their genetic or epigenetic features could be identified.

The SNPs identified in T-H or T-P human Refseq genes were aligned. Among the 57 identified SNPs in these genes, 35 were present in either the T-P or the T-H individuals (Table II). In addition, 29 shared, and 4 T-H- and 2 T-P-specific SNPs were detected that were annotated as non-synonymous. None of these SNPs were considered rare variants. Finally, the methylation analysis did not reveal any significant differential methylation peaks (P<10⁻²⁰) in the promoter or gene body of any of the 57 known CS-associated candidate genes.

ASM may directly affect the pattern of gene expression. Several DNA regions acquire ASM during development and often in a tissue-specific manner, whereas others are methylated during gametogenesis, which is stably maintained throughout development (3,22,33–36). To determine the role of ASM in CS pathogenesis, the possible differences in ASM in the twins were investigated by assessing the methylation status of different alleles.

The exon-capture MeDIP-seq reads were used to identify the heterozygous variants between the twins. A total of 14,535 heterozygous variants were identified and subsequently the allelic bias in the MeDIP-seq read coverage for each individual was determined by χ² test (P<0.01). In the ASM regions, 68 and 54 SNPs were detected in T-H and T-P, respectively. In addition, the twins shared 11 of these SNPs (Table III). Among them, 1 SNP, located on chr10:29839864, was detected in the SVIL gene. The SVIL gene was of interest, as DAVID analysis (37,38) indicated that it was related to Gene Ontology (GO) terms associated with skeletal tissue development, including cytoskeleton organization (GO:0007010), muscle organ development (GO:0007517) and skeletal muscle tissue development (GO:0007519) (35,37–39). In T-H, one allele containing a C was highly methylated compared with the other allele, whereas in T-P the allele containing an A showed comparatively higher methylation. The aforementioned observation indicates that the pathogenesis of CS is associated with more limited alterations in the methylation status of alleles and not to wider ones. These minor changes in the methylation status lead to differential allelic expression between the twins that affected their development.