Variants of ACAN are associated with severity of lumbar disc herniation in patients with chronic low back pain

Introduction Disc herniation is a complex spinal disorder associated with disability and high healthcare cost. Lumbar disc herniation is strongly associated with disc degeneration. Candidate genes of the aggrecan metabolic pathway may associate with the severity of lumbar disc herniation. Objectives This study evaluated the association of single nucleotide variants (SNVs) of the candidate genes of the aggrecan metabolic pathway with the severity of lumbar disc herniation in patients with chronic mechanical low back pain. In addition, we assessed the in-silico functional analysis of the significant SNVs and association of their haplotypes with the severity of lumbar disc herniation. Methods A descriptive cross sectional study was carried out on 106 patients. Severity of disc herniation and disc degeneration were assessed on T2-weighted mid sagittal lumbar MRI scan. Sixty two exonic SNVs of ten candidate genes of aggrecan metabolic pathway (ACAN, IL1A, IL1B, IL6, MMP3, ADAMTS4, ADAMTS5, TIMP1, TIMP2 and TIMP3) were genotyped on a Sequenom MassARRAY iPLEX platform. Multivariable linear regression analysis was carried out using PLINK 1.9 software adjusting for age, gender, body mass index and severity of disc degeneration. Four online bioinformatics tools (Provean, SIFT, PolyPhen and Mutation Taster) were used for in-silico functional analysis. Results Mean age was 52.42 ± 9.42 years and 69.8% were females. The mean severity of disc herniation was 2.81 ± 1.98. The rs2272023, rs35430524, rs2882676, rs2351491, rs938609, rs3825994, rs1042630, rs698621 and rs3817428 variants and their haplotypes of ACAN were associated with the severity of lumbar disc herniation. However, only the rs35430524, rs938609 and rs3817428 variants of ACAN were detected as pathogenic by in-silico functional analysis. Conclusions SNVs of ACAN and their haplotypes are associated with the severity of lumbar disc herniation. Functional genetic studies are necessary to identify the role of these significant SNVs in the pathogenesis of disc herniation.


Introduction
Disc herniation is a complex spinal disorder associated with disability and high healthcare cost. Lumbar disc herniation is strongly associated with disc degeneration. Candidate genes of the aggrecan metabolic pathway may associate with the severity of lumbar disc herniation.

Objectives
This study evaluated the association of single nucleotide variants (SNVs) of the candidate genes of the aggrecan metabolic pathway with the severity of lumbar disc herniation in patients with chronic mechanical low back pain. In addition, we assessed the in-silico functional analysis of the significant SNVs and association of their haplotypes with the severity of lumbar disc herniation.

Methods
A descriptive cross sectional study was carried out on 106 patients. Severity of disc herniation and disc degeneration were assessed on T2-weighted mid sagittal lumbar MRI scan. Sixty two exonic SNVs of ten candidate genes of aggrecan metabolic pathway (ACAN, IL1A, IL1B, IL6, MMP3, ADAMTS4, ADAMTS5, TIMP1, TIMP2 and TIMP3) were genotyped on a Sequenom MassARRAY iPLEX platform. Multivariable linear regression analysis was carried out using PLINK 1.9 software adjusting for age, gender, body mass index and severity of disc degeneration. Four online bioinformatics tools (Provean, SIFT, PolyPhen and Mutation Taster) were used for in-silico functional analysis. PLOS

Introduction
Lumbar disc herniation is a complex multifactorial spinal condition associated with disability, work time loss and high health related costs [1]. Lumbar disc degeneration is strongly associated with disc herniation [2]. In traditional view, lumbar disc degeneration and herniation are mainly related to age, gender, body mass index (BMI), smoking, physical activity and heavy loading of the spine. Significant associations were identified in our previous study between single nucleotide variants (SNVs) of candidate genes of the aggrecan metabolic pathway and lumbar disc degeneration [3]. Therefore, it is worthwhile to explore the association of SNVs of respective candidate genes (ACAN, IL1A, IL1B, IL6, MMP3, ADAMTS4, ADAMTS5, TIMP1, TIMP2 and TIMP3) with the severity of disc herniation. Aggrecan is the main proteoglycan of the intervertebral disc and provides osmotic properties which assist the disc in resisting mechanical compressive loads transmitted along the spine [4]. With aging, the disc undergoes dehydration due to loss of proteoglycan with aging and the degenerative process is further augmented by repetitive damages. This transforms the disc into a degenerative state and probably disc herniation [5]. Interleukins (ILs), matrix metalloproteinases (MMPs), a disintegrin and metalloproteinase with thrombospondin motifs (ADAMTSs) and tissue inhibitor of metalloproteinases (TIMPs) are important molecules in the aggrecan metabolic pathway [6].
Abnormal genetic variants of the candidate genes in the aggrecan metabolic pathway may alter the structural and functional properties of the aggrecan molecule. Although there are several studies which have identified associations of a Variable Number of Tandem Repeats (VNTR) polymorphism of ACAN with lumbar disc degeneration and symptomatic lumbar disc herniation [7], there is a lack of studies which have explored the association of SNVs of ACAN with the severity of disc herniation. The "T" allele of the rs1042631 variant of ACAN is associated with signal intensity of the disc and increased the frequency of disc herniation among the Finnish population [8]. Although the MMP1 and MMP3 are associated with disc herniation [9,10], relationships of SNVs of the ADAMTS4, ADAMTS5, TIMP1, TIMP2 and TIMP3 with the severity of disc herniation need to be explored.
Early detection of genetically predisposed individuals for disc herniation might provide opportunity to take necessary precautions in lifestyle related risk factors. Identification of significant SNVs/haplotypes provides probable therapeutic targets for gene therapy in the future.
Considering the findings of the previous genetic study on lumbar disc degeneration [3] and current evidence available in the literature, we extended our research to investigate the associations of candidate genes of the aggrecan metabolic pathway with the severity of lumbar disc herniation. The main objective of this study was to determine the associations of single nucleotide variants of the candidate genes of the aggrecan metabolic pathway with the severity of disc herniation in patients with chronic mechanical low back pain. In addition, we assessed in-silico functional analysis of significant SNVs of the candidate genes and associations of their haplotypes with the severity of lumbar disc herniation.

Study design, setting and participants
A descriptive cross-sectional study was carried out on patients with chronic mechanical low back pain who attended the rheumatology clinic, National Hospital of Sri Lanka, Colombo from May 2012 to October 2014. 368 consecutive patients with chronic mechanical low back pain who fulfilled the eligibility criteria were assessed with lateral x-rays of lumbar spine after obtaining written informed consent. Male and female patients (20 to 69 years) having pain/muscle tension or stiffness localised below the costal margin and inferior gluteal folds [11], during day time and worsening in the latter part of the day due to movements [12] for at least three months duration [13] were included in the study. Patients with back pain due to inflammatory causes, visceral origin, systemic infections affecting spine, metabolic bone diseases, fractures in the vertebral column, past surgeries in the spine, and spinal tumours were excluded. 120 patients were selected to undergo MRI scan of lumbar spine and genotyping based on the severity of disc related degenerative changes (disc space narrowing and anterior osteophyte) in the x-ray of lumbar spine [3,14]. The study was carried out in accordance with the Declaration of Helsinki and with the approval of the Ethics Review Committee of the Faculty of Medicine, University of Colombo.

Clinical evaluation
An interviewer administered questionnaire was used to record age, gender and clinical examination was used to record body mass index (BMI). Height (cm) and weight (kg) of the patients were assessed with light clothing and without shoes to the nearest 0.1 cm and 0.1 kg, respectively, and BMI was calculated (kg/m 2 ) [15]. International cut off values were used for categorisation of BMI (normal, overweight and obese) [16].

X-ray assessment
Lateral lumbar x-rays were carried out while patients were in lateral recumbent position on the table flexing the knees and hips just enough to achieve comfortable position [17]. The intervertebral disc spaces (L1/L2 to L5/S1) in lateral lumbar x-rays were assessed for disc space narrowing and anterior osteophytes by a consultant radiologist blinded to clinical details of the patients and overall lumbar disc degeneration was calculated (grade 0-2) according to a scoring system defined by Lane et al. 1993 (Fig 1) [14].

MRI assessment
Patients with grade 2 or more lumbar disc degeneration in the lateral lumbar x-rays were selected to undergo MRI scan of the lumbar spine. In addition age and gender matched samples from patients with grade 0 and grade 1 lumbar disc degeneration were selected to undergo MRI scans of the lumbar spine. T2 weighted sagittal MRI scans of lumbar spine were carried out in supine position using a GE 1.5T MRI Scanner (Signa, General Electric, Milwaukee, Wisconsin) and assessed for the severity of lumbar disc degeneration and disc herniation.

Severity of lumbar disc degeneration
The severity of lumbar disc degeneration was graded on T2-weighted mid sagittal MRI images using modified Pfirrmann grading system and each lumbar level was graded from grade 1 to 5 (Table 1) [18,19]. Grade 2 and 3 was defined as mild lumbar disc degeneration, while grades 4 and 5 were defined as severe lumbar disc degeneration.The grades of lumbar disc degeneration of the five lumbar levels were summed to calculate the severity of lumbar disc degeneration.

Severity of lumbar disc herniation
Disc herniation was considered present when the disc material is displaced beyond the disc space. It was further categorised as disc protrusion when the distance between the edges of the herniated disc material extending outside the disc space is less than the distance between the Arrows-A-no disc space narrowing/anterior osteophyte (grade 0 lumbar disc degeneration), B-mild disc space narrowing and small anterior osteophyte (grade 1 lumbar disc degeneration), C-small anterior osteophyte and moderate disc space narrowing (grade 2 lumbar disc degeneration) [3]. edges of the base of that disc material. Disc extrusion was considered present when the distance between the edges of the herniated disc material extending outside the disc space was more than the distance between the edges of the base of that disc material (Fig 2) [20]. Presence of disc protrusion in a disc was given a score of 1 and presence of disc extrusion in a disc was given a score of 2. The scores of each lumbar level were summed to calculate the severity of disc herniation.

SNV selection, DNA extraction and genotyping
The selected molecules of the aggrecan metabolic pathway are aggrecan, interleukin 1α, interleukin 1β, interleukin 6, matrix metalloproteinases 3, a disintegrin and metalloproteinase with thrombospondin motifs 4 and 5, tissue inhibitor of metalloproteinase 1, 2 and 3. Their functions and respective gene symbols are summarised in the Table 1 in the previous article [3]. Common SNVs in the exonic regions of these ten candidate genes in the aggrecan metabolic pathway (ACAN, IL1A, IL1B, IL6, MMP3, ADAMTS4, ADAMTS5, TIMP1, TIMP2 and TIMP3) were identified from 25 Sri Lankan exomes and Gujarati Indians in Houston (GIH) in the Hapmap database [3]. Sixty two SNVs in exonic regions of the candidate genes in aggrecan metabolic pathway were selected for genotyping (

Statistical analysis
Descriptive statistics were used to summarise the sample characteristics. Chi-square test was used to assess the association between two categorical variables while Pearson's correlation test was used to assess association between two continuous variables. SNVs with Hardy-Weinberg equilibrium (HWE) ! 0.05 were included in the genetic association analysis. Quantitative trait association of SNVs was carried out with the multivariable linear regression analysis using PLINK 1.9 software [24]. Severity of disc herniation was used as the quantitative outcome. The genotype of the respective SNV was used as the main independent variable and was treated as a quantitative variable and coded 0, 1 or 2 to represent the number of variant allele, consistent with an additive genetic model. Haplotype frequencies of significant SNVs were performed using PLINK 1.07 software [24], which uses the expectation maximization algorithm to determine common haplotypes. Multivariable linear regression model was used to assess the quantitative trait association of haplotypes of significant SNVs of the candidate genes with the severity of disc herniation. Age, gender, BMI and highest grade of lumbar disc degeneration out of five lumbar levels were used as the additional covariates. Permutation function in PLINK was used with 10,000 permutations to generate the significance levels empirically. It helps to relax the assumptions of normality and correct the errors due to small sample size. P value < 0.05 was used as the level of significance. Four free online bioinformatics prediction tools (Provean [25], SIFT [26], PolyPhen [27] and Mutation taster [28]) were used to conduct the in-silico functional analysis.

Results
Although 120 patients were selected to undergo MRI scans of lumbar spine, only 106 patients attended for MRI scans of lumbar spine. Characteristics of the patients who underwent both the MRI scan of lumbar spine and genotyping (106 patients) are summarised in Table 3. The proportion of patients increased up to the 50-59 years age group and declined after age of 60 years. The majority of patients were females. Most of the patients were overweight and obese. Mean severity of disc herniation was 2.81 ± 1.98. The majority of patients had disc herniations with 45 (42.4%) having disc extrusions. Among a total of 530 lumbar levels evaluated, 72 (13.6%) and 69 (13.0%) of the discs had severe lumbar disc degeneration and disc extrusions, respectively. Severe lumbar disc degeneration and disc extrusions were more frequent in lower lumbar levels (65.3% and 66.7%, respectively). Severity of disc herniation was positively correlated with the severity of disc degeneration (correlation coefficient = 0.49, p = <0.001). Furthermore, percentage of disc extrusions were higher in discs with severe disc degeneration (33.3%) compared to the discs with mild disc degeneration (9.8%) (χ 2 (1) = 30.36, p < 0.001).
Overall genotyping efficiency of these patients was 98.68%. Minor allele frequency ranged from 0.05 to 0.50 and all the variants were in HWE ( Table 2).

Associations of SNVs of ACAN with severity of lumbar disc herniation
Nine SNVs of ACAN were significantly associated with the severity of disc herniation ( Table 4). Presence of each additional "C" allele of the rs2272023 variant, "A" allele of the rs35430524 variant, "A" allele of the rs2882676 variant and "A" allele of the rs1042630 variant of ACAN were associated with a progressive increase in the severity of disc herniation. Furthermore, the presence of each additional "T" allele of the rs2351491 variant, "A" allele of the rs938609 variant, "G" allele of the rs3825994 variant and "G" allele of the rs698621 variant of ACAN were associated with progressive reduction in the severity of disc herniation. SNVs, rs2351491, rs938609 and rs698621 variants were in linkage disequilibrium (R 2 ! 0.93, D' ! 0.98). The rs3817428 variant of ACAN did not have homozygous genotype for the minor allele and the presence of "G" allele was associated with reduction in the severity of disc herniation. However, there was no significant association between the SNVs of the other candidate genes and the severity of disc herniation (Table 4).       In-silico functional analysis of the significant variants associated with the severity of lumbar disc herniation Among the nine significant SNVs only three (rs35430524, rs938609 and rs3817428) were predicted as pathogenic by the in-silico functional analysis ( Table 5). The rs3817428 variant of ACAN was identified as pathogenic by two prediction tools (Provean and PolyPhen) while the other two were detected by one of either Provean or PolyPhen prediction tools. Both the rs35430524 and rs3817428 variants were located in conserved regions of ACAN. In addition, the rs35430524 variant had a non-conservative amino acid substitution.

Associations of haplotypes of ACAN with the severity of lumbar disc herniation
Among the nine significant SNVs, rs938609 and rs2882676 formed three haplotypes while rs3825994 and rs1042630 also formed three haplotypes. The haplotypes frequencies and results of multivariable linear regression analysis are summarised in Table 6. In haplotype analysis, the T-A haplotype of rs938609 and rs2882676 loci and T-A haplotype of rs3825994 and rs1042630 loci increased the severity of disc herniation. In contrast A-C haplotype of rs938609 and rs2882676 loci and G-G haplotype of rs3825994 and rs1042630 loci reduced the severity of disc herniation after adjusting for confounders (age, gender, BMI and the severity of lumbar disc degeneration).

Discussion
The study assessed associations of SNVs of candidate genes of the aggrecan metabolic pathway with the severity of lumbar disc herniation in patients with chronic mechanical low back pain.
In addition, we assessed the in-silico functional analysis of significant SNVs and associations of their haplotypes with the severity of lumbar disc herniation. In our results, minor allele of the rs2272023 (C), rs35430524 (A), rs2882676 (A), rs2351491 (T), rs938609 (A), rs3825994 (G), rs1042630 (A), rs698621 (G) and rs3817428 (G) variants of ACAN were associated with the severity of lumbar disc herniation. Furthermore, the rs35430524, rs938609 and rs3817428 variants were detected as pathogenic by the in-silico functional analysis. In addition, haplotypes of rs938609 and rs2882676 loci, and haplotypes of rs3825994 and rs1042630 loci of ACAN were associated with the severity of lumbar disc herniation. In disc degeneration, the concentration of aggrecan is reduced due to increased breakdown and reduced synthesis of aggrecan molecules leading to disc dehydration [6]. This process is further augmented by concurrent changes of the proteins involved in the aggrecan metabolic pathway (interleukins, matrix metalloproteinases, a disintegrin and metalloproteinase with thrombospondin motifs and tissue inhibitors of metalloproteinases). Annular fissures appear in the disc with the loss of integrity of annulus fibrosus and disc dehydration. Nuclear material can leak through these fissures and displace beyond the disc space margins leading to disc herniation. Disc herniation is categorised as disc protrusion and extrusion based on the extension of herniated disc material beyond the disc space [20]. In our sample, almost all the patients had disc herniations with 42.4% having disc extrusions.
Out of sixty two selected SNVs of ten candidate genes, only nine SNVs of ACAN were associated with the severity of disc herniation which included five nonsynonymous variants (rs35430524, rs938609, rs2882676, rs1042630 and rs3817428) and four synonymous variants (rs2272023, rs2351491, rs3825994 and rs698621). A nonsynonymous variant alters the amino acid sequence of the respective protein. This altered amino acid sequence may affect the function of the respective protein. In contrast, synonymous variants do not alter the amino acid sequence of the respective protein, but it can cause changes in long non-coding RNA, microRNA and promoters causing the disease by affecting the protein expression, conformation and function. Furthermore any nonsynonymous or synonymous variant may act as a marker of a functional SNV in the same gene or nearby gene [29]. VNTR polymorphisms reported as significantly associated with the disease are likely to be found as such, because they are in linkage disequilibrium with pathogenic SNVs [30]. Therefore, the focus of this study was to assess the SNVs associated with the severity of disc herniation.
There are a limited number of studies which have explored the association of SNVs of the ACAN with the severity of disc herniation. The presence of each "T" allele of the rs1042631 variant of ACAN increased the frequency of disc herniation among the Finnish population [8]; it also increased the risk of annular tears in Indian patients who attended a spine unit of a tertiary care hospital [31]. However, this variant was not associated with the severity of disc herniation in our sample. The rs2351491, rs938609 and rs698621 variants of ACAN were in strong linkage disequilibrium and each additional minor allele progressively reduced the severity of disc herniation ( Table 4). The rs938609 variant is located in exon 12 of ACAN and causes a Serine to Threonine substitution at position 939 of the amino acid sequence. This substitution is a very highly conservative substitution, but only the Polyphen bioinformatics tool predicted it as pathogenic ( Table 5). The rs2351491 and rs698621 variants are synonymous variants and functionally unremarkable.
The presence of each additional minor allele of the rs35430524 (A), rs1042630 (A) and rs2882676 (A) variants of ACAN was associated with a progressive increase in the severity of disc herniation while the presence of an additional "G" allele of the rs3817428 variant of ACAN reduced the severity of disc herniation. All these four variants were nonsynonymous, but only the rs35430524 and rs3817428 variants were predicted as pathogenic. The rs35430524 variant is located in exon 12 of ACAN and causes a Proline to Threonine substitution at position 913 of the amino acid sequence. It is located at a highly-conserved region of ACAN and is not a conservative substitution. Therefore there is a high probability that the rs35430524 variant is pathogenic and which was identified as such by the Provean prediction tool ( Table 5). The rs3817428 variant located in exon 15 of ACAN causes an Aspartic acid to Glutamic acid substitution at position 2373 of the amino acid sequence. This is located in a highly-conserved area in ACAN, but its amino acid substitution is a highly conservative substitution. Two of the four bioinformatics tools (Provean and Polyphen) predicted that the rs3817428 variant is pathogenic ( Table 5). The rs3825994 and rs2272023 variants of ACAN are synonymous variants and the results of the functional prediction tools were unremarkable.
Our findings on the associations of SNVs of ACAN with the severity of lumbar disc herniation were further strengthened by the results of the haplotype analysis. In haplotype analysis, we found that T-A (rs938609-rs2882676) and T-A (rs3825994-rs1042630) haplotypes increased the severity of disc herniation while A-C (rs938609-rs2882676) and G-G (rs3825994-rs1042630) haplotypes reduced the severity of disc herniation. This suggests that the severity of disc herniation is influenced by single nucleotide variation and also by haplotype differences.
Several studies reported significant associations of SNVs of IL6, IL10 [32], MMP1 and MMP3 [9] with disc herniation. Some of these variants were intronic variants [32] and some variants were not polymorphic in our population. We have assessed only the exonic variants and the SNVs of the candidate genes of inflammatory, catabolic, anti-catabolic molecules were not associated with the severity of disc herniation (Table 4). In our previous study, SNVs of IL1A (inflammatory gene), ADAMTS4 and ADAMTS5 (catabolic genes) were found to be associated with the severity of disc degeneration [3], but these SNVs were not associated with the severity of disc herniation. Thus, findings of the two studies suggest that SNVs of candidate genes in inflammatory and catabolic pathways are responsible for the variation in the severity of disc degeneration while SNVs of ACAN (structural gene) are responsible for the variation in the severity of disc herniation.
Quantitative traits like disc herniation may have significant effects from variants of several candidate genes and environmental factors. However, we have not corrected the genetic association analysis for the multiple variants we already assessed and other probable variants in other respective genes. Although the rs3817428, rs35430524 and rs938609 variants of ACAN were predicted as pathogenic by the in-silico functional analysis, overall results of the four prediction tools were inconsistent. Considering these facts, it is important to carry out genetic association analysis with a larger sample number using whole genome sequencing followed by functional genetic studies to identify the exact effect of the significant SNVs.

Limitations of the study
The results of this study are not generalizable as it was conducted on a specific group of patients with chronic mechanical low back pain. However, there were patients representing the three districts of the Western Province of Sri Lanka. Low sample number is another limitation of the study, but we used permutation function of the PLINK software tool to relax the assumptions of normality and correct the errors due to small sample size. There may be other confounding factors that we have not corrected for.

Conclusions
Single nucleotide variants of ACAN and their haplotypes are associated with the severity of lumbar disc herniation. The rs35430524, rs938609 and rs3817428 variants of ACAN are predicted as pathogenic by in-silico functional analysis. However functional involvement of these three SNVs should be confirmed with functional genetic studies. This study provides basic insight of genetic markers which could be used as probable therapeutic targets for novel treatment strategies for lumbar disc herniation.