The study protocol was approved by the ethical committee for human subjects of the Peking University Third Hospital. Informed consent was provided by all participating individuals. A total of 30 unrelated northern Chinese Han T-OPLL patients with myelopathy and/or neurological dysfunction [14 men, (mean age, 51.71±6.38 years); 16 women (mean age, 52.63±5.97 years)] and 5 unrelated healthy controls (mean age, 52.20±1.30 years, use for analysis of susceptibility gene expression levels in peripheral blood) were enrolled in this study from February 2010 to July 2016. Diagnosis of T-OPLL was performed by specialists based on clinical symptoms and radiologic examinations (including CT and magnetic resonance imaging) of the thoracic spine. The appearance of OPLL observed in radiographs was classified into four subtypes: i) Segmental; ii) continuous; iii) mixed; and iv) local (17). Neurological status was evaluated by the Japanese Orthopedic Association (JOA) score for thoracic myelopathy (a total of 11 points). Individuals who had lumbar spondylolisthesis, ankylosing spondylitis, diffuse idiopathic skeletal hyperostosis, and disc herniation of the thoracic spines were excluded in this study and did not take any drugs known to affect bone or calcium metabolism.

Genomic DNA was extracted from peripheral blood leukocytes using a standard method. To discover genetic variations, we performed WGS on 30 unrelated northern Chinese Han patients. Quality genomic DNA from the 30 samples was fragmented using a Covaris ultrasonicator. By adjusting shearing parameters, DNA fragments from each sample were concentrated in 500-bp peaks. These fragments were purified, end blunted, ‘A’ tailed, and adaptor ligated. DNA templates with adapters were then selectively enriched using polymerase chain reaction (PCR) to obtain a sufficient amount for the DNA library. The concentration of DNA for the libraries was quantified using a bioanalyser (Agilent Technologies, Inc., Santa Clara, CA, USA) and qPCR. Each qualified DNA library was sequenced on an Illumina HiSeq Xten platform using paired-end reads according to the manufacturer's instructions. Sequencing-derived raw image files were processed by Illumina base calling Software with default parameters, and sequence data for each individual were generated as paired-end reads, defined as ‘raw data’.

Bioinformatics analysis began with sequencing data (raw data from the Illumina machine). First, clean data were produced by data filtering of raw data. All clean data for each sample were mapped to the human reference genome (GRCh37/HG19). Burrows-Wheeler Aligner (18,19) software was used to do the alignment. To ensure accurate variant calling, we followed recommended best practices for variant analysis with the Genome Analysis Toolkit. Local realignment around InDels and base quality score recalibration were performed using GATK (20,21), with duplicate reads removed by Picard tools. The sequencing depth and coverage for each individual were calculated based on the alignments. In addition, a strict data analysis QC system throughout the entire pipeline was built to guarantee high-quality sequencing data.

To decrease noise in the sequencing data, data filtering was carried out, including: i) removal of reads containing the sequencing adapter; ii) removal of reads with a low-quality base ratio (base quality less than or equal to 5) that was more than 50%; and iii) removal of reads with an unknown base (‘N’ base) ratio that was more than 10%. Statistical analysis of data and downstream bioinformatics analyses were performed on the filtered, high-quality data, referred to as the ‘clean data’. Variant frequency was compared with the 1000G SNP database (http://www.1000genomes.org/). Potential deleterious effects of identified sequence variants were assessed by various algorithms, such as SIFT (http://sift.jcvi.org/) (22), PolyPhen-2 (http://genetics.bwh.harvard.edu/pph2/) (23), MutationTaster (http://www.mutationtaster.org/) (24), and GERP++ (http://mendel.stanford.edu/SidowLab/downloads/gerp/) (25).

Sanger DNA sequencing (ABI 3730XL; Applied Biosystems; Thermo Fisher Scientific, Inc., Waltham, MA, USA) was used to confirm the accuracy of variants identified by WGS. The PCR fragments were submitted for Sanger sequencing at the Beijing Genomics Institute. Details of the two studied SNPs and the primer sequences are listed in Table I. PCR was performed with 20 ng genomic DNA per 15 µl reaction mixture, containing 0.2 µM of each primer, 200 µM of deoxyribonucleotides, 50 mM KCl, 10 mM Tris HCl (pH 8.3), 1.5 mM MgCl₂ and 0.5 units of Taq DNA polymerase in a DNA Gradient PCR machine (Bio-Rad Laboratories, Inc., Hercules, CA, USA). The thermocycling conditions were as follows: initial denaturation at 95°C for 10 min; followed by 35 cycles of 95°C for 30 sec, annealing at an assay-specific temperature (48 to 65°C) for 45 sec, and elongation at 72°C for 45 sec; and a final terminal elongation step at 72°C for 5 min. The PCR products were analyzed by direct sequencing using a BigDye Terminator v3.1 Cycle Sequencing kit (Thermo Fisher Scientific, Inc.) with POP-7™ Polymer in a 3730XL DNA Analyzer with Sequencing Analysis Software version 5.2 (Thermo Fisher Scientific, Inc.).

CLUSTAL W (http://www.genome.jp/tools/clustalw/) was used to compare homologous protein sequences among multiple species to analyse the consistency of amino acid sequences, particularly mutation sites, in seven representative species.

Plasma collection and storage were carried out using standard methods. Plasma COL6A1 (cat. no. HG21134) and IL17RC (cat. no. HG1762) levels were quantified using commercially available ELISA kits (Trust Specialty Zeal, San Francisco, CA, USA). All samples were assayed according to the manufacturer's instructions and were run in duplicate. The optical density of each well was determined using a microplate reader at 450 nm. No interference and no cross reactivity were expected based on the manufacturer's instructions.

Total RNA was purified from blood using the SK1321 RNA Blood Mini kit (Sangon Biotech Co., Ltd., Shanghai, China). An on-column DNase digest (Sangon Biotech Co., Ltd.) was performed before the clean-up step to eliminate residual genomic DNA. cDNA was synthesized from total RNA (2 µg) using a RevertAid Premium Reverse Transcriptase kit from Thermo Fisher Scientific, Inc. Relative quantitative RT-PCR was applied to quantify the mRNAs levels of COL6A1, IL17RC and GAPDH using SYBR-Green Real-Time PCR master mix on the LightCycler480 Real-Time System from Roche (Basel, Switzerland). All experiments were performed in triplicate and normalized to GAPDH.

All statistical analyses were performed using SPSS v17.0 software (SPSS, Inc., Chicago, IL, USA). Descriptive data for continuous variables are presented as the mean ± standard deviation. Student's t-tests were used to determine the age and JOA score differences between patients with or without OPLL gene mutations. The OPLL subtype differences between patients with or without OPLL gene mutations were determined using one-way analysis of variance with a post hoc Fisher's test. P<0.05 was considered to indicate a statistically significant difference.

After standard SNP quality control, we performed WGS on 30 DNA samples with an average of 119,583.62 Mb raw bases from the Illumina Hiseq Xen sequencer. After removing low-quality reads, we obtained an average of 793,363,211 clean reads (119,004.48 Mb). The clean reads of each sample had high Q20 and Q30, indicating high sequencing quality. The average GC content was 41.01% (Table II). Strict data quality control (QC) was performed during the entire analysis pipeline to obtain accurate data (Table III).

To prioritize potential pathogenic variants, we focused on the identification of rare [multiple allele frequency (MAF) ≤0.005, based on the BGI database] and damaging variants predicted by at least two algorithms (e.g., SIFT, Polyphen-2, MutationTaster, and GERP++). Two deleterious variants were identified in seven unrelated patients (Table IV), and these findings were further confirmed by directional Sanger sequencing (Fig. 1). Indicating that the test results for WGS are accurate.

c.1534G>A (p.Gly512Ser) was located in the exon regions of 42.9% of patients with OPLL. Moreover, this variant was found to be evolutionarily conserved in three primates. c.2275C>A (p.Leu759Ile) was located in intron regions of 57.1% of patients with OPLL. However, this variant was not found in six other vertebrate species (Fig. 2). These two heterozygous mutations were found to be relatively conserved, suggesting highly genetic heterogeneities and complex pathogenesis.

To elucidate the functional roles of loci and genes associated with T-OPLL, we measured the expression levels of genes in three patients carrying p.Gly512Ser mutations in COL6A1 and four patients carrying p.Leu759Ile mutations in IL17RC compared with 5 healthy controls by enzyme-linked immunosorbent assay (ELISA). We found that plasma COL6A1 concentrations (20.90±0.64 µg/l) were significantly higher (P<0.05), approximately 4-fold higher than those in healthy controls (4.79±1.18 µg/l). Moreover, plasma IL17RC concentrations (6.13±0.36 µg/l) were also significantly higher (P<0.05), approximately 2-fold higher than those in healthy controls (2.89±0.30 µg/l). We performed RT-qPCR analysis using peripheral blood cells and found that COL6A1 mRNA levels in three patients carrying the p.Gly512Ser mutation were approximately 6-fold higher than those in healthy controls, the difference was statistically significant (P<0.05). The IL17RC mRNA levels in four patients carrying the p.Leu759Ile mutation were approximately 4-fold higher than those in the healthy controls, the difference was also statistically significant (P<0.05). Compared with healthy controls group, we observed that these two mutations significantly increase their respective gene expressed, suggesting that these two potential pathogenic loci have potential effect on their respective gene expressed cells.

Phenotype-genotype correlations were analysed among the seven patients with missense mutations and 23 patients without significant mutations (Table V). No differences were found between these two groups in terms of sex or age at diagnosis. Two-dimensional computed tomography (2D-CT) scans of seven patients with OPLL harbouring deleterious variants and one healthy control are shown in Figs. 3 and 4. More specifically, one patient belonged to the Local subtype (Fig. 3A), one patient belonged to the Mixed subtype (Fig. 3B), one patient showed the Segmental subtype (Fig. 3C), and four patients showed the Continuous subtype (Figs. 3D and 4A-C), Fig. 4D showed the healthy control. However, the JOA score for thoracic myelopathy was significantly lower in patients with rare missense mutations compared with patients without mutations (P=0.001). Additionally, radiological analysis of OPLL morphology (17) revealed that the frequency of the continuous subtype was significantly higher in patients with rare missense mutations than in patients without mutations (P=0.033). Moreover, four mutation-positive patients (57.1%) showed the continuous subtype, one patient (14.3%) showed the local subtype, one patient (14.3%) showed the segmental subtype, and one patient (14.3%) showed the mixed subtype. By contrast, three (13%) mutation-negative patients showed the continuous subtype, six patients (26.1%) showed the local subtype, 10 patients (43.4%) showed the segmental subtype, and four patients (17.4%) without mutations showed the mixed subtype.