DOI: https://doi.org/10.21203/rs.3.rs-41775/v1
Degenerative lumbar spinal stenosis (DLSS) is a common lumbar disease requiring surgery. Previous studies indicate that genetic mutation may play a part in DLSS. However, there is scarcely specific study regarding the genetic study for DLSS. Whole exome sequencing (WES) is an effective strategy to search for disease-causing genes. To identify probable susceptible genes, we performed a WES in a cohort of 50 patients and 25 controls with DLSS. By bioinformatics analyses, we successfully identified several novel variants including HLA-DRB1, PARK2, ACTR8, AOAH, BCORL1, MKRN2, NRG4, NUP205, etc. To the best of our knowledge, this is the first genetic study of susceptibility genes with DLSS by WES in Chinese patients. In conclusion, our findings revealed that deleterious mutations in these genes may contribute to the development of DLSS.
Lumbar spinal stenosis (LSS) is defined as any type of narrowing of the lumbar spinal canal, the lateral recesses, or the intervertebral foramina. It may be congenital and acquired or degenerative (1, 2). Congenital stenosis is rare and typically caused by achondroplasia and hypochondroplasia. Degenerative lumbar spinal stenosis (DLSS) is usually attributable to degenerative changes in the intervertebral discs, facet joints, and ligamentum flavum, spondylolisthesis. DLSS typically affects subjects over 50 years of age, as one of the most common spinal conditions worldwide that often requires surgery. Neurogenic claudication is the classic clinical manifestation associated with lumbar spinal stenosis causing diminished function and impaired quality of life (3, 4).
To date, the etiology and pathomechanism of DLSS remain unclear, and previous studies indicate that genetic factors may play a part in DLSS. A previous Finnish study identified a splice site variant in COL9A2 was associated with DLSS (5). Seung-Jae Hyun reported a COL9A2 haplotype (HAP2) was significantly associated with DLSS in the Korean population, whereas another haplotype (HAP4) may play a protective role against DLSS development (6). However, there is inadequate study regarding the susceptibility gene mutation for DLSS. It is urgent to identify pathogenic genes to provide novel approaches for future investigation and intervention of DLSS.
Whole-exome sequencing (WES) represents a cost-effective, reproducible, and robust strategy for the sensitive and specific identification of variants causing protein-coding changes in individual human genomes. Previous genetic studies were limited to implementing one or several candidate genes approach rather than the more established measure of WES study, whereby hundreds of thousands of single nucleotide polymorphisms (SNPs) can be screened (7). For the population with pathogenic gene mutations, preventive measures such as lifestyle modification, and better monitoring for stenosis development may be applied. In this study, we analyzed data of WES on fifty DLSS patients and twenty-five healthy controls to attempt the identification of pathogenic genes. We used Sanger DNA sequencing to validate the variants disclosed by WES. Single nucleotide polymorphisms were analyzed by bioinformatics analysis. Our study reported several pathogenic variants in genes associated with DLSS in the Chinese population.
A total of 50 unrelated Chinese patients with symptoms consistent with DLSS and 25 unrelated healthy controls were consecutively recruited between the year 2016 January and 2017 December from the Orthopaedic Department in China-Japan Friendship Hospital. This study was conducted with the approval of the ethics committee of the China-Japan Friendship Hospital Institutional Review Board. Written informed consent was obtained from all patients whose specimens and clinical information were used for this study. All patients were examined both clinically and radiologically to confirm the degenerative nature of the stenosis. All DLSS patients met the all following three criteria: (1) diagnosis of DLSS by magnetic resonance imaging (MRI); (2) conservative treatment and monitoring for more than 3 months by spinal surgeons; and (3) a history of typical LSS symptoms, consisting of self-reported intermittent neurogenic claudication, stenotic symptoms extending to the lower extremities upon extension of the lumbar spine, or numbness or weakness of the lower extremities. Patients with neurological disorders, spinal fracture, spondylolisthesis, spinal tumor, presentation with trauma, and diagnosis with an infectious disease were excluded. 25 biologically unrelated healthy individuals, with a similar ethnic background, were enrolled as the control group.
Genome DNA was extracted from whole blood with TIANamp Blood DNA kit (TIANGEN BIOTECH,Beijing, China). Exome capture was performed using Agilent SureSelect Human All Exon V6 kit (Agilent Technologies, Santa Clara, CA) as per the manufacturer’s instructions. A whole exome sequencing was performed on an Illumina HiSeq Xten platform. Sequencingderived 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’ The experiment was performed according to manufacturer’s protocol. One sample was loaded per flow cell lane to obtain a minimum of 10 × read depth across ~ 96% of the target regions.
Bioinformatics analysis started with the raw sequencing data from the Illumina pipeline. Raw data were processed and filtered and low-quality reads were discarded. Low-quality sequencing reads were filtered by the following criteria:1) removal of reads containing the sequencing adapter, 2) removal of reads with a low-quality base ratio (base quality less than or equal to 5 that was more than 50%; 3) removal of reads with an unknown base (‘N’ base) ratio that was more than 10%. Variants located in exon regions or alternative splicing regions were retained. Variant frequency was compared with the 1000G SNP database (http://www.1000genomes.org/). Exome Aggregation Consortium (ExAC) (http://exac.broadinstitute.org;) Analysis for potential deleterious mutations was performed using various algorithms, such as PolyPhen2 (http://genetics.bwh.harvard.edu/pph2/) (8), SIFT (http://sif.jcvi.org/) (9), Mutation Taster (http://www.mutationtaster.org/) (10) as well as CADD (Combined Annotation Dependent Depletion). We used the identification of rare combined with at least two of the four algorithms above to prioritize potential pathogenic variants. R software (http://www.bioconductor.org/packages/release/bioc/html/maftools.html) was used to perform the clustering analysis to determine statistical significance in genotype and allele frequencies between patients and controls (FDR, false discovery rate < 0.05). Phenolyzer (Phenotype Based Gene Analyzer, http://phenolyzer.wglab.org/) was used to discover genes based on user-specific disease/phenotype terms.
Toppgene suite was used for gene list enrichment analysis and candidate gene prioritization based on functional annotations and protein interactions network (https://toppgene.cchmc.org/). A web-based portal named Metascape (http://metascape.org/) was used to do Gene Ontology (GO) enrichment analysis. GSEA software (version 3.0) with “c2.cp.kegg.v6.1symbols.gmt” gene sets database was used to do Gene Set Enrichment Analysis (GSEA) (http://www.broadinstitute.org/gsea/index.jsp), to identify significantly enriched or depleted genes that show statistically significant, concordant differences between two given clusters.
A Fisher's exact test was used for evaluating the association between rare variants and disease phenotype (case/control data) with Bonferroni correction to P-value.
Twenty-one male and 29 female patients were included in this study, with an average age of 52.4 years ranging from 46 to 70 years. There are 10 male and 15 female controls in this study, with an average age of 52.4 years ranging from 29 to 70 years. After Principal Component Analysis (PCA), three records were obviously outliers. Without the three outliers, all the 72 cases are in random distribution. The outliers three (bijinfang, bijinfen, bijinhua) have closer genetic relationship. (Supplementary Fig. 1 and Fig. 2)
We performed WES on seventy-five DNA samples from the Illumina Hiseq Xen sequencer. To filter potential pathogenic variants, we focused on the identification of rare 1000G_EAS ≤ 0.005, based on BGI database) and damaging variants predicted by at least two of four algorithms (i.e. SIFT, Polyphen2, Mutation Assessor, Mutation Taster, and CADD). Detrimental variants are classified into missense_mutation, frame_shift_del, nonsense_mutation, frame_shift_ins, splice_site, in_frame_ins, in_frame_del, nonstop_mutation. The basic variant status in the 75 cases were shown in Supplementary Fig. 3 and Fig. 4.
Genes related to LDSS by clustering analysis by R software. We filtered allele mutations in controls to retain only mutations in patients and got multiply SNPs on thirty-four genes. (FDR, false discovery rate < 0.05.) (Fig. 1 and Table 1)
GENE |
Feature_1 |
Feature_2 |
N_of_Feature1 |
N_of_Feature2 |
FDR |
---|---|---|---|---|---|
GPRIN2 |
Patient |
Normal |
45 of 50 |
0 of 25 |
2.70983E-15 |
MYOT |
Patient |
Normal |
45 of 50 |
0 of 25 |
2.70983E-15 |
PDE4DIP |
Patient |
Normal |
45 of 50 |
0 of 25 |
2.70983E-15 |
EPPK1 |
Patient |
Normal |
44 of 50 |
0 of 25 |
1.40008E-14 |
SYT15 |
Patient |
Normal |
39 of 50 |
0 of 25 |
1.14246E-11 |
CGREF1 |
Patient |
Normal |
31 of 50 |
0 of 25 |
3.2011E-08 |
IRF5 |
Patient |
Normal |
25 of 50 |
0 of 25 |
3.3725E-06 |
SIRPB1 |
Patient |
Normal |
24 of 50 |
0 of 25 |
9.2172E-06 |
HLA-DQA1 |
Patient |
Normal |
23 of 50 |
0 of 25 |
1.06146E-05 |
MUC6 |
Patient |
Normal |
19 of 50 |
0 of 25 |
0.000141694 |
PROSER3 |
Patient |
Normal |
18 of 50 |
0 of 25 |
0.000315482 |
IGFN1 |
Patient |
Normal |
17 of 50 |
0 of 25 |
0.000375298 |
PABPC1 |
Patient |
Normal |
16 of 50 |
0 of 25 |
0.000728519 |
KRT79 |
Patient |
Normal |
9 of 50 |
0 of 25 |
0.02508378 |
FLNA |
Patient |
Normal |
9 of 50 |
0 of 25 |
0.02508378 |
PHRF1 |
Patient |
Normal |
9 of 50 |
0 of 25 |
0.02508378 |
ALS2CL |
Patient |
Normal |
10 of 50 |
0 of 25 |
0.025483837 |
KNSTRN |
Patient |
Normal |
10 of 50 |
0 of 25 |
0.025483837 |
COL6A5 |
Patient |
Normal |
8 of 50 |
0 of 25 |
0.046170366 |
TSHZ2 |
Patient |
Normal |
8 of 50 |
0 of 25 |
0.046170366 |
AADAC |
Patient |
Normal |
8 of 50 |
0 of 25 |
0.046170366 |
ADAT1 |
Patient |
Normal |
8 of 50 |
0 of 25 |
0.046170366 |
CMA1 |
Patient |
Normal |
8 of 50 |
0 of 25 |
0.046170366 |
DEAF1 |
Patient |
Normal |
8 of 50 |
0 of 25 |
0.046170366 |
DLAT |
Patient |
Normal |
8 of 50 |
0 of 25 |
0.046170366 |
DST |
Patient |
Normal |
8 of 50 |
0 of 25 |
0.046170366 |
FLG |
Patient |
Normal |
8 of 50 |
0 of 25 |
0.046170366 |
KIAA0556 |
Patient |
Normal |
8 of 50 |
0 of 25 |
0.046170366 |
LPA |
Patient |
Normal |
8 of 50 |
0 of 25 |
0.046170366 |
P2RX5 |
Patient |
Normal |
8 of 50 |
0 of 25 |
0.046170366 |
PANK2 |
Patient |
Normal |
8 of 50 |
0 of 25 |
0.046170366 |
RASSF4 |
Patient |
Normal |
8 of 50 |
0 of 25 |
0.046170366 |
WARS |
Patient |
Normal |
8 of 50 |
0 of 25 |
0.046170366 |
ZNF572 |
Patient |
Normal |
8 of 50 |
0 of 25 |
0.046170366 |
Frequency of allele mutations were zero in controls and FDR < 0.05. |
To decipher the specified characteristics between patients and control groups, we analyzed the differences by Gene Ontology (GO) and pathway analysis. (Fig. 2 and Supplementary Table 1)
Fifty cases and 2245 control cases (Novegene database) were used to do the sites-based association. (Fig. 3) Forty-three genes were filtered up to P-value ≤ 0.001, FDR_BH_Allele ≤ 0.01. After toppgene test, the top two genes were HLA-DRB1 and PRKN. The top four SNP related to lumbar disease were shown in Table 2. We perform pathway and Phenolyzer to discover genes based on lumbar disease. Genes related to LDSS (Phenolyzer score ≥ 0.01) were showed in Supplementary Fig. 5.
Chromosome ID |
Ref |
Alt |
Alt Case |
Alt Control |
Ref Case |
Ref Control |
P-value |
Function |
|
---|---|---|---|---|---|---|---|---|---|
6 |
rs113663708 |
G |
T |
14 |
199 |
78 |
4151 |
0.0001 |
exonic |
6 |
rs201219028 |
G |
A |
29 |
724 |
63 |
3616 |
0.000602 |
splicing |
6 |
rs9269955 |
G |
C |
40 |
1110 |
58 |
3328 |
0.000883 |
exonic; splicing |
6 |
rs12190823 |
G |
A |
42 |
1179 |
50 |
3289 |
0.000109 |
intronic |
After filtering the data by 1000G_EAS ≤ 0.005 and damaging variants predicted by at least two of four algorithms (i.e. SIFT, Polyphen2, Mutation Assessor, Mutation Taster, and CADD), there are 322 SNP on 60 candidate genes which P-value ≤ 0.001. We perform pathway and Phenolyzer to discover genes based on lumbar disease. As ‘lumbar degeneration’ is not in the disease item, we input lumbar disc degeneration as the disease item in Phenolyzer. Genes related to LDSS (Phenolyzer score ≥ 0.01) were showed in Table 3 and Supplementary Fig. 6.
Gene |
Alt_Case |
Ref_Case |
Alt_Control |
Ref_Control |
P-value |
---|---|---|---|---|---|
PLEC |
9 |
91 |
0 |
4490 |
7.69E-16 |
SHC4 |
5 |
95 |
31 |
4459 |
0.000979097 |
NRG4 |
3 |
97 |
4 |
4486 |
0.00033 |
TNKS |
2 |
98 |
0 |
4490 |
0.00047 |
DKK4 |
2 |
98 |
0 |
4490 |
0.00047 |
RHPN1 |
2 |
98 |
0 |
4490 |
0.00047 |
Phenolyzer score ≥ 0.01 |
As there are genetic relations among bijingfang, bijinfen, and bijinhua, not all the patients are sporadic cases. We tried to find SNPs that not only in the three cases but also in the other 47 cases as more as possible. Six genes were identified to have relations with LDSS with high possibility. (Table 4 and Table 5)
bijinfang |
|||
---|---|---|---|
bijinfang |
bijinfen |
||
bijinfang |
bijinfen |
bijinhua |
|
nonsynonymous cSNP,splice site variatant or coding indel(NS/SS/I) |
1060 |
657 |
87 |
NS/SS/I not dbSNP |
133 |
76 |
30 |
NS/SS/I not 1000g2015aug_all |
133 |
76 |
30 |
NS/SS/I not 1000g2015aug_Chinese |
133 |
76 |
30 |
NS/SS/I not esp6500siv2_all |
133 |
76 |
30 |
NS/SS/I not NovoDb_WES_SNP |
113 |
63 |
19 |
NS/SS/I not ExAC |
113 |
63 |
19 |
29 |
15 |
6 |
47 Case |
50 case |
|||
---|---|---|---|---|
Gene |
1000g2015aug_all < 0.01 |
ExAC_ALL < 0.01 |
1000g2015aug_all < 0.01 |
ExAC_ALL < 0.01 |
ACTR8 |
13 |
1 |
16 |
4 |
AOAH |
13 |
10 |
16 |
13 |
BCORL1 |
3 |
3 |
6 |
6 |
MKRN2 |
3 |
3 |
6 |
6 |
NRG4 |
6 |
6 |
9 |
9 |
NUP205 |
2 |
2 |
5 |
5 |
To date, several lines of evidence have shown that genetic factors may play a part in DLSS. We performed a WES study design on 50 patients and 25 controls to detect genes contributing to DLSS. In our study, we discover several novel candidate genes that have not been previously identified in relation to DLSS.
Prior studies have identified several candidate genes in association with lumbar disc degeneration. Jason et al. carried out a genome-wide association study, and multiple SNPs identified suggested a multi-factorial origin of DLSS (11). HLA-DRB1 genotyping was reported to increase the risk of developing pain after surgery or lumbar disc herniation (12). A variant in the PARK2 gene associated with Lumbar disc degeneration by influencing methylation of the PARK2 (13). In this study, HLA-DRB1 and PARK2 are identified as susceptibility genes associated with a predisposition to DLSS. Lumbar disc degeneration causes intervertebral collapse, which may accelerate the development of DLSS. Obviously, there could be some overlapped D:\æœç‹—émechanism between the two degenerative processes. To some extent, our results are in accordance with the previous studies. The specific mechanism of the two genes contributing to DLSS needs further investigation.
In this study, ACTR8, AOAH, BCORL1, MKRN2, NRG4, NUP205 were novel genes reported in correlation with DLSS. ACTR8 is located on chromosome 3p21.1 and has 14 exons. Mutation in the ACTR8 was associated with lineage-specific expression in primates (14). AOAH is on chromosome 7p14.2, encoding both the light and heavy subunits of acyloxyacyl hydrolaseregion. AOAH polymorphisms have been reported to associate with asthma, chronic rhinosinusitis, and bronchial hyperreactivity (15, 16). Located on chromosome Xq26.1, BCORL1 encoding a transcriptional corepressor that is found tethered to promoter regions by DNA-binding proteins. Pathogenic variants in BCORL1 underlie a newly identified X-linked epigenetic syndrome (17). MKRN2, on chromosome 3p25.2, encodes a probable E3 ubiquitin ligase containing several zinc finger domains. MKRN2 is involved in the regulation of inflammatory response, non-small-cell lung cancer (18, 19). NRG4 is located on chromosome 15q24.2. NRG4 is a member of the epidermal growth factor (EGF) family of extracellular ligands, highly expressed in adipose tissues, enriched in brown fat, and markedly increased during brown adipocyte differentiation (20). NUP205 on chromosome ch7q33 encoding a nucleoporin, which is a subunit of the nuclear pore complex that functions in the active transport of proteins, RNAs, and ribonucleoprotein particles between the nucleus and cytoplasm. Mutations in NUP205 are associated with steroid-resistant nephrotic syndrome (21).
Several other candidate genes associated with DLSS were identified including GPRIN2, MYOT, PDE4DIP, etc. Some pathways were enriched through analysis of the differentially expressed genes between the patient and control groups. These genes had never been reported to be associated with DLSS,
We recognize several limitations of this study. First, the number (50) is still inadequate for genetic study. Concerning the inadequate number, we enroll three family cases in this study. Numerous candidate genes identified in this study warrant further investigation. Functional study should be conducted to determine the mechanism and pathways in which the genes involved affect the development of the disease. As a next step, we will continue increasing the sample size to reveal more susceptible genes.
In conclusion, in this study, we identifed mutations with DLSS in Chinese patients for the first time using WES. From our analysis, we found two potential pathogenic loci for DLSS on HLA-DRB1 and PARK2, which have previously been reported to be associated with lumbar disc degeneration, and SNPS on several novel genes including ACTR8, AOAH, BCORL1, MKRN2, NRG4, NUP205, etc, which had not previously been reported to be DLSS. Further replication and functional studies are required to confirm these positive findings, helping to detect early and to provide novel therapeutic approaches of DLSS.
Not applicable
The authors declare that they have no conflict of interest.
The clinical study was approved by the ethics committee of China-Japan Friendship Hospital and was conducted in accordance with the provisions of the Declaration of Helsinki. Written informed consent was obtained from all participants before enrolment.
The data used to support the findings of this study are available from the corresponding author upon request.
Not applicable
Xin Jiang and Dong Chen Conceived, coordinated, designed, drafted and edited the manuscript.