Two GWAS-identified variants are associated with lumbar spinal stenosis and Gasdermin-C expression in Chinese population

The aim of this study is to investigate the expression levels of genome-wide association studies (GWAS)-identified variants near Gasdermin-C (GSDMC) and its association with lumbar disc degeneration (LDD) in a Chinese population. In accordance with previously reported findings, our study involved the top 4 variants; rs6651255, rs7833174, rs4130415, and rs7816342. A total of 800 participants, 400 LDD patients and 400 controls were involved in the study. The LDD patients were divided into two mutually exclusive subgroups: subgroup 1: lumbar disc herniation; subgroup 2: lumbar spinal stenosis. Genotyping were performed using TaqMan assay, and Enzyme-Linked Immunosorbent Assay (ELISA) used to measure the plasma GSDMC levels, while quantitative reverse-transcription (qRT)-PCR and immunohistochemistry (IHC) were used to evaluate the GSDMC expression levels. Among the studied variants, there were no statistically significant differences in allelic and genotypic frequencies between LDD patients and their controls (all P > 0.05). However, the subgroup analysis revealed a significant association between rs6651255 and rs7833174 in patients with lumbar spinal stenosis (subgroup 2). Furthermore, the max-statistic test revealed that the inheritance models of two variants of lumbar spinal stenosis were represented by the recessive model. The plasma and mRNA expression levels of GSDMC were significantly higher in patients with lumbar spinal stenosis compared with the control group (P < 0.05). Furthermore, the CC genotypes of rs6651255 and rs7833174 were significantly associated with increased plasma expression levels of GSDMC in patients with lumbar spinal stenosis (P < 0.01). Two GWAS-identified variants (rs6651255 and rs7833174) near GSDMC were associated with a predisposition to lumbar spinal stenosis. GSDMC protein and mRNA expression levels may have prognostic qualities as biomarkers for the existence, occurrence or development of lumbar spinal stenosis.

candidate genes for LDD 21,22 . Many GWAS-associated variants are located in non-regulatory regions of the genome, thus, increasing the difficulty in assessing the underlying molecular mechanism of the nucleotide variation 23 . The replication of the association in different ethnic groups is important to validate the results of the GWAS 24 . Based on the existing GWAS findings, we focused on rs6651255, rs7833174, rs4130415 and rs7816342, which are the top 4 variants at 8q24.21 near GSDMC gene. To the best of our knowledge, this is the first study to investigate the relationship between of these GWAS-identified variants and LDD in a Chinese population. In this study, Enzyme-Linked Immunosorbent Assay (ELISA), quantitative reverse-transcription (qRT)-PCR and immunohistochemistry (IHC) were used to evaluate the plasma and intervertebral disc expression levels of GSDMC in LDD patients and healthy controls.

Results
Characteristics of subjects. The characteristics of the study and control groups are summarized in Table 1. There were no statistically significant differences in age, gender, BMI, smoking history, and alcohol consumption between the study and control groups. The median age of the control group was 50.1 years, and ranged between 20.0 and 60.0 years, while that of the case group was 49.6 years, ranging from 18.0 to 69.0 years. There was no statistically significant difference between the case and control groups (P > 0.05). The disc samples from 94 patients (32 from subgroup 1; 28 from subgroup 2; and 34 from the control group) were assigned for qRT-PCR and IHC analysis. No statistically significant differences in gender, age and BMI were observed among three groups (all P > 0.05).
Case-control association analysis. The genotype frequencies of the top 4 variants at 8q24.21 locus (rs6651255, rs7833174, rs4130415 and rs7816342) are shown in Table 2. There was no significant deviation of allelic frequencies from the Hardy-Weinberg equilibrium (HWE) in control group (all P > 0.05) among the studied variants. With regards to the rs6651255 and rs7833174 variants, there were no statistically significant differences in genotypic frequency between the LDD groups and the control (P > 0.05) ( Table 2), however, the subgroup analysis revealed a significant association between them and lumbar spinal stenosis (Table 3). Furthermore, the use of the max-statistic test and genetic model selection analysis revealed that the inheritance models of two variants of lumbar spinal stenosis were represented by the recessive model. Thus, it was proposed that the CC genotypes of rs6651255 and rs7833174 possibly increased the susceptibility to lumbar spinal stenosis, on the contrary, the TT genotypes of rs6651255 and rs7833174 were protective. There were no statistically significant difference between the allelic or genotypic frequencies of rs4130415 and rs7816342 between the study and control groups. Furthermore, there were no statistically significant difference in rs4130415 and rs7816342 between the 2 subgroups (all P > 0.05).
Plasma GSDMC levels and LDD. The plasma expression levels of GSDMC for the case and control groups are shown in Fig. 1. The mean plasma GSDMC level for the control group was 592.33 ng/L (range 198.35-1270.26 ng/L), while that of the case groups was 726.77 ng/L (range 187.27-1833.20 ng/L). In case groups, the levels for subgroups 1 and 2 were 671.44 ng/L (range 187.27-1536.31 ng/L) and 837.42 ng/L (range 280.72-1833.20 ng/L) respectively. There was a statistically significant difference in plasma GSDMC levels between case and control (P < 0.01) as well as subgroup 2 and the control (P < 0.01), but not between subgroup 1 and the control (P > 0.05) (Fig. 1a,b). Thus, indicating that increased plasma GSDMC levels were significantly observed in patients with lumbar spinal stenosis. It was also observed that the plasma expression levels of GSDMC in subgroup 2 with rs6651255 or rs7833174 CC genotypes were significantly higher than those with TT + TC (all P < 0.05) (Fig. 1c,d).
GSDMC expressions and LDD. The results of the qRT-PCR and IHC were consistent with our previously reported plasma findings. It was observed that the GSDMC mRNA level were increased 1.8 fold in subgroup 2 ( Fig. 2e) when compared with the control, and the difference was statistically significant (P < 0.05). However, there were no statistically significant difference between the GSDMC mRNA level in subgroup 1 and control group (P > 0.05). The results of the IHC analysis also revealed that the GSDMC expressions levels were significantly higher in subgroup 2 than the control group (immunopositive cells: 62.19 ± 6.14% vs. 34.87 ± 3.34% respectively; P < 0.001) (Fig. 2a,c,d). However, the difference in expression levels between subgroup 1 and the Table 1. Characteristics of the study subjects between healthy controls and all cases including two subgroups. a All the cases composed of the subgroup 1 and subgroup 2.

Characteristics
Healthy controls All

Discussion
To the best of our knowledge, this is the first study to investigate the possible relationship between the top 4 GWAS-identified variants (rs6651255, rs7833174, rs4130415, rs7816342) and LDD in a Chinese population. Our results revealed that 2 of the top 4 variants, i.e. rs6651255 and rs7833174 were associated with increased susceptibility to lumbar spinal stenosis. The CC genotypes of both rs6651255 and rs7833174 were found to assert a high genetic predisposition to lumbar spinal stenosis. Furthermore, we identified significantly increased plasma GSDMC levels and related gene expressions in patients with lumbar spinal stenosis.  17 . According to the genome-wide significance thresholds, their findings identified 4 variants among the 37 variants (rs6651255, p = 5.61 × 10 -12 ; rs7833174, p = 5.97 × 10 -12 ; rs4130415, p = 5.99 × 10 -12 ; and rs7816342, p = 7.81 × 10 -12 ) to be highly correlated with LDH. A subsequent meta-GWAS on chronic back pain involving 158,000 individuals of European ancestry further confirmed rs7833174 to be strongly associated with chronic back pain (P = 4.4 × 10 -13 ) 25 . A recent study by Wu et al. 26 evaluated the relationship between other GSDMC gene variations (rs4527833, rs77681114, rs4285452, rs4733741 and rs4509280) and LDH risk, and identified that rs77681114 was significantly associated with a decreased risk of LDH. However, these findings were inconsistent with our study as our results found neither the allelic nor the genotypic frequencies of these variants (rs6651255, rs7833174, rs4130415 and rs7816342) to be significantly different between our LDD group and control. Our analysis of the overall cohort of LDD patients resulted in their classification into two mutually exclusive subgroups in accordance with their presented clinical condition; i.e. disc herniation or spinal stenosis. This approach allowed for better definition and investigation Table 3. Genotype frequencies of rs6651255, rs7833174, rs4130415 and rs7816342 in controls and subgroups. Subgroup 1 patients with lumbar disc herniation, Subgroup 2 patients with lumbar spinal stenosis, OR odds ratio, CI confidence interval, MAX3 max-statistic test, GMS genetic model selection. The www.nature.com/scientificreports/ into which of the relationships among the observed variants were true for all LDD patients or our predefined subgroups and how our subgroups of particular lumbar spine pathologies interacted. Interestingly, the results of our subgroup analysis of the variants rs6651255 and rs7833174, revealed statistically significant differences in both the allelic and genotypic distributions between subgroup 2 (lumbar spinal stenosis) and the control. The genetic heterogeneity among different ethnic populations may be a major explanation for different results. The previously reported studies involved subjects from European backgrounds 17,25 , whiles the current study focused on an ethnic Chinese population. The ethnic and genetic diversity coupled with social and cultural differences may reflect in both genotypic and phenotypic heterogeneity between Chinese populations and other nationalities. Furthermore, the term LDD as a general phenotype has been used synonymously with a multitude of terminologies, such as disc herniation, spinal stenosis and spondylolisthesis. LDH happens as a result of aging and the break down that occurs within the intervertebral disc 27 . In contrast, spinal stenosis not only is caused by intervertebral disc degeneration but may also arise because of degenerative changes in the facet joints, and/ or ligamentum flavum 28 . Clinically, spinal stenosis was usually defined as a specific phenotype secondary to disc herniation, which at this point represented the advanced stage of LDD. Bjornsdottir et al. 17 found that these variants were associated with LDH. In this study, significant association was only observed in spinal stenosis subgroup, but not in LDH subgroup. This discrepancy could be explained by the fact that LDD has a complex and heterogeneous pathogenesis [29][30][31] . Therefore, from our results, we postulated that these variants may exert more influence in the advanced stages of LDD development in Chinese population. The closest protein coding gene for rs6651255, rs7833174, rs4130415 and rs7816342 at 8q24.2 is GSDMC, and as a member of the Gasdermin superfamily, it is primarily expressed in the trachea and spleen and thought to be involved in the regulation of apoptosis and immune-related functions 18,19,32 . The immune system has been reported to be activated in the development of pathogenic disc degeneration, and numerous studies have confirmed the migration of immune cells into the nucleus pulpous and increased levels of inflammatory cytokines in degenerated discs 33,34 . However, the relationship between rs6651255, rs7833174 variants and GSDMC expression has not been previously reported. In the current study, we observed the plasma levels GSDMC were significantly higher in lumbar spinal stenosis subgroup compared to the control. Correspondingly, the results of the qRT-PCR and IHC revealed elevated expression levels of GSDMC in degenerated discs. Thus, these results indicate that www.nature.com/scientificreports/ elevated expression levels of GSDMC to could signify an increased risk for lumbar spinal stenosis. Moreover, our investigation into whether these variants could affect plasma levels of GSDMC revealed that the CC genotypes of rs6651255 and rs7833174 were significantly associated with higher expression levels of GSDMC compared to the TT + TC genotypes. In summary, our findings identified the CC genotypes of rs6651255 and rs7833174 as a possible genetic risk factor indicator for development of LDD, and probably involved in increasing the gene expression levels of GSDMC. However, the exact molecular methodology by which it influences LDD development requires further investigations. There are several limitations to this study. Firstly, the sample size is relatively small, hence diminishing the statistical power to detect the much smaller effects of the studied genetic variations in our subgroup analysis. Secondly, we could since only participants who consented to the study were from a single institution were involved in the study, there may be some level of the potential selection bias. Future studies in other populations, preferably multicentered with lager study population could be beneficial in confirming our findings.
In conclusion, two GWAS-identified variants near GSDMC (rs6651255 and rs7833174) may be related with a predisposition to lumbar spinal stenosis. And their CC genotypes may be significantly associated with increased plasma GSDMC expression level in patients with lumbar spinal stenosis, and thus could be useful as diagnostic or prognostic biomarkers for the existence, occurrence or as a risk factor for the development of lumbar spinal stenosis.

Materials and methods
Ethic statement. This study was approved by the Ethics Committee of the First Affiliated Hospital of Guangxi Medical University (2018-KY-NSFC-025). Written informed consent was obtained from all participants involved, and all procedures were performed in accordance with the guidelines of the institutional research committees and the Helsinki Declaration.
Subjects. This study involved 800 participants, including 400 sporadic LDD patients (188 males, 212 females; mean age 49.6 ± 18.3 years, range 18-69 years) and 400 healthy controls (194 males, 206 females; mean age 50.2 ± 16.3 years, range 20-60 years). All participants were recruited from the Department of Spine and Osteopathic Surgery of the First Affiliated Hospital of Guangxi Medical University. All patients were examined and diagnosed by two recognized spine surgeons via physical examination and MRI. Participants were included in the study if: (1) Have history of low back pain for at least 3 months; (2) MRI scans revealed degenerative changes to the lumbar spine; (3) Have no previous history of spinal surgery or other treatments that could cause deformation of lumbar vertebrae. The control group included healthy participants with matched age and sex without history of back pain who were randomly recruited from the Physical Examination Centre of the First Affiliated Hospital of Guangxi Medical University. In accordance with the modified Pfirrmann grading system 35 , healthy subjects with Pfirrmann Grade 1 were included in the study as the control group. The LDD patients were divided into two mutually exclusive subgroups based on their MRI phenotypes (Fig. 3), with subgroup 1 involving 286 patients with lumbar disc herniation, and subgroup 2 involving 114 patients with lumbar spinal stenosis. Additionally, herniated disc tissues (n = 32), degenerative disc tissues (n = 28) and normal disc tissues (n = 34) were collected from patients with lumbar disc herniation (subgroup 1), lumbar spinal stenosis (subgroup 2) and traumatic lumbar vertebral fracture patients respectively. The patients with traumatic lumbar fracture had no previous preoperative history of low back pain, and their MRI scans revealed no disc degeneration. According to the Schneiderman's classification 36 (Table S1). These samples were used to evaluate the GSDMC gene expression via qRT-PCR and IHC. All participants were of Zhuang ethnicity from Guangxi Zhuang Autonomous Region of China.
DNA isolation and genotyping. In accordance with established protocols (MagNA Pure LC, Roche Applied Science, Indianapolis, IN, USA), Genomic DNA of each participant were isolated from 5 mL peripheral blood leukocytes. Qubit 3.0 fluorometer (Invitrogen) was used for DNA quantification while NanoDrop2000 spectrophotometer (Thermo Fisher Scientific, FL, USA) was used to assess DNA purity. The variants were selected based on their strong associations with LDD in GWAS 17 , and the top 4 variants among the 37 highly correlated variants i.e. rs6651255, rs7833174, rs4130415 and rs7816342 selected. The primers, probes and reaction conditions used in this study are available upon request (Table S2). TaqMan variant Genotyping Assays in an ABI 7900 sequence detection system (Applied Biosystems, Foster City, CA, USA) was used for variant genotyping, and all investigators involved blinded to the status of the subjects. Ten percent of the samples were randomly selected as duplicate quality control samples. ELISA detection of plasma GSDMC level. The blood samples were obtained from both LDD groups and the control group, and the plasma stored at − 80 °C after centrifugation. The serum GSDMC concentrations were measured using ELISA (Abcam, Cambridge, MA, USA) in accordance with the manufacturer's procedural instructions. The sensitivity of the assay method was less than 2 pg/mL, and there was no cross-reactivity with other cytokines. The intra-assay and inter-assay coefficients of variation were no more than 10%, and all samples were analyzed twice.
RNA extraction and qRT-PCR analysis. The intervertebral disc tissues were lysed in TRIzol (Invitrogen Inc, Carlsbad, CA, USA) and Rneasy Mini Kit (Qiagen, Valencia, CA, USA) used for total RNA extraction in accordance with manufacturer's protocol. The reverse transcriptions (RT) were performed using PrimeScript RT Master Mix kit (Takara, Japan), with 1 μg total RNA used for the synthesis of the complementary DNA (cDNA) Statistical analysis. Differences between allelic frequencies, genotype distributions and gene expressions between cases and controls were analyzed using the standard χ 2 -analysis, with Hardy-Weinberg equilibrium used for goodness-of-fit χ 2 test. The differences in allelic frequencies were evaluated by calculating the odd ratios and 95% confidence intervals. In order to avoid multiple comparisons by fitting three genetic models and determining the best-fitting model among them, we tested the global P value via two robust tests, max-statistic test (MAX3) 37 and genetic model selection (GMS) 38 . For ELISA, qRT-PCR and IHC experiments, unpaired Student t-test was used for comparing the means of two groups, and one-way analysis of variance (ANOVA) with Turkey's post-hoc test used for comparing the means of multiple groups. SPSS software Version 17.0 (Chicago, USA) and R software (Version 3.3.3, https ://www.r-proje ct.org/) with SNPassoc package 39 were used for the data analysis, and all results considered statistically significant when P < 0.05. Patient consent. Obtained.