Single Nucleotide Variants of Candidate Genes in Aggrecan Metabolic Pathway Are Associated with Lumbar Disc Degeneration and Modic Changes

Introduction Lumbar disc degeneration (LDD) is genetically determined and severity of LDD is associated with Modic changes. Aggrecan is a major proteoglycan in the intervertebral disc and end plate. Progressive reduction of aggrecan is a main feature of LDD and Modic changes. Objectives The study investigated the associations of single nucleotide variants (SNVs) of candidate genes in the aggrecan metabolic pathway with the severity of LDD and Modic changes. In-silico functional analysis of significant SNVs was also assessed. Methods A descriptive cross sectional study was carried out on 106 patients with chronic mechanical low back pain. T1, T2 sagittal lumbar MRI scans were used to assess the severity of LDD and Modic changes. 62 SNVs in ten candidate genes (ACAN, IL1A, IL1B, IL6, MMP3, ADAMTS4, ADAMTS5, TIMP1, TIMP2 and TIMP3) were genotyped on Sequenom MassARRAY iPLEX platform. Multiple linear regression analysis was carried out using PLINK 1.9 in accordance with additive genetic model. In-silico functional analysis was carried out using Provean, SIFT, PolyPhen and Mutation Taster. Results Mean age was 52.42±9.42 years. 74 (69.8%) were females. The rs2856836, rs1304037, rs17561 and rs1800587 variants of the IL1A gene were associated with the severity of LDD and Modic changes. The rs41270041 variant of the ADAMTS4 gene and the rs226794 variant of the ADAMTS5 gene were associated with severity of LDD while the rs34884997 variant of the ADAMTS4 gene, the rs55933916 variant of the ADAMTS5 gene and the rs9862 variant of the TIMP3 gene were associated with severity of Modic changes. The rs17561 variant of the IL1A gene was predicted as pathogenic by the PolyPhen prediction tool. Conclusions SNVs of candidate genes in ACAN metabolic pathway are associated with severity of LDD and Modic changes in patients with chronic mechanical low back pain. Predictions of in-silico functional analysis of significant SNVs are inconsistent.


Introduction
Chronic low back pain is the leading cause of years lived with disability in most regions of the world including South Asia [1]. Most of the chronic low back pain is due to mechanical causes such as strains, sprains and complications related to lumbar disc degeneration (LDD) [2]. Lumbar intervertebral disc consists of central nucleus pulposus, surrounding annulus fibrosus and cartilage end plate. LDD is characterised by the structural failure and cell mediated response to the changes in the composition of the disc [3]. Modic changes are types of degenerative changes which involve the vertebral end plate and bone marrow [4] and are also associated with the severity of LDD and chronic low back pain [5]. There is a wide variation in the severity of LDD/Modic changes in the population as some individuals develop early and severe types of LDD/Modic changes compared to other individuals in the same age category. With the recent advancement in genetic research, it has been suggested that these degenerative changes are genetically determined and are modified to some degree by behavioural and environmental factors [6,7].
Genetics can determine the size, shape and mechanical properties of the spinal structures including the intervertebral disc. LDD is characterised by structural damage to the disc matrix and increased activity of matrix degrading enzymes [8]. In addition, persistent inflammation and subsequent vascular and neural responses are key features of painful LDD and Modic changes [9,10]. Therefore the structural components and molecules in their degradation pathways are ideal candidates for genetic association studies. Aggrecan (coded by the ACAN gene) is a large molecular weight proteoglycan with numerous glycosaminoglycan side chains and is the core proteoglycan in the disc matrix and vertebral end plate. A Variable Number of Tandem Repeat (VNTR) polymorphism in exon 12 of the ACAN gene is associated with LDD where individuals with shorter number of tandem repeats are at higher risk for severe LDD [11][12][13][14]. There are a few studies which have assessed the associations of single nucleotide variants (SNVs) of the ACAN gene with the severity of LDD, however most of them are intronic variants [15,16].
Candidate genes involved in metabolic pathways of LDD and Modic changes (eg. interleukins, matrix metalloproteinases, aggrecanases and tissue inhibitors of metalloproteinases) are also implicated with the onset and degree of mechanical failure of the intervertebral disc and vertebral end plate. Interleukins (IL1α, IL1β and IL6) are a group of local cytokines involved in the regulation of the inflammatory response. The "T" allele of the rs1800587 variant of the IL1A gene (codes the IL1α molecule) is associated with early LDD in girls [17] and Modic changes [18]. Matrix metalloproteinase 3 (MMP3 and is coded by the MMP3 gene) is the most frequently studied metalloproteinase and the 5A allele of the rs3025058 variant of the MMP3 gene is associated with progression and severity of LDD [19][20][21]. A disintegrin and metalloproteinase with thrombospondin motifs 4 and 5 (ADAMTS 4 and ADAMTS 5 and are coded by the ADAMTS4 and ADAMTS5 genes, respectively) are the main aggrecanases involved in degrading the protein core of aggrecan [22]. The "T" allele of the rs4233367 variant of the ADAMTS4 has a lower risk for LDD [23]. In addition, the "A" allele of the rs151058 variant, "A" allele of the rs229052 variant and "A" allele of the rs162502 variant of the ADAMTS5 gene (intronic variants) are associated with LDD [24]. Tissue inhibitors of metalloproteinases (TIMPs) are a group of protease inhibitors which inhibit the catabolic actions of MMPs and ADAMTSs [25] and the C124T variant of the TIMP1 gene is associated with radiographic progression of LDD [20].
Intensity of pain and severity of disability of the chronic mechanical low back pain are correlated with severity of LDD [26] and Modic changes [5]. Severity of LDD and Modic changes might be affected by the SNVs of the molecules involved in the metabolic pathway of the ACAN gene. However most of the identified SNVs associated with LDD/Modic changes are located at the intronic regions of the DNA sequence and evidence for SNVs in exonic regions of the respective candidate genes is limited. ADAMTSs and TIMPs are relatively recently found molecules and there are very few studies which have investigated the associations of SNVs of genes of these molecules with the severity of LDD/Modic changes [25]. Although the associations of SNVs of candidate genes are identified, evidence of their functional involvement in LDD/Modic changes are not well established. The main objective of this study was to determine the associations of single nucleotide variants of candidate genes in the aggrecan metabolic pathway, with the severity of LDD and Modic changes in patients with chronic mechanical low back pain. We also assessed the in-silico functional analysis of the significant SNVs associated with the severity of LDD/Modic changes.

Study design, setting and participants
A descriptive cross-sectional study was conducted on patients with chronic mechanical low back pain who attended the rheumatology clinic, National Hospital of Sri Lanka, Colombo. Both male and female patients of Sri Lankan origin with chronic mechanical low back pain aged 20 to 69 years were included in the study based on the severity of disc related degenerative changes in the x-rays of lumbar spine. Low back pain was defined as pain, muscle tension, or stiffness localized below the costal margin and above the inferior gluteal folds, with or without radiating pain to the leg [27]. Low back pain during day time worsening in the latter part of the day due to movements such as bending, lifting, walking, running, standing and sitting were considered as being due to a mechanical cause [28,29]. Chronicity was defined as pain on most days of the week of at least three months duration [30]. 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. Pregnant females and patients who refused to participate in the study were also excluded.
Three hundred and sixty eight patients with chronic mechanical low back pain who attended the rheumatology clinic, National Hospital of Sri Lanka, Colombo from May 2012 to October 2014 underwent lateral x-ray of lumbar spine and x-rays were evaluated by a consultant radiologist blinded to the clinical details of the patients. Disc space narrowing and anterior osteophytes (L1/L2 to L5/S1) were graded and overall LDD was derived using a scoring system defined by Lane et al. 1993 (Fig 1). Each lumbar level was graded from grade 0 to grade 2 LDD [31]. Hundred and twenty patients were selected to undergo MRI scan of the lumbar spine based on following criteria. Patients with grade 2 LDD on lateral x-ray of lumbar spine were further assessed using mid-sagittal MRI scans of the lumbar spine. In addition, age and gender matched patients from other two groups of LDD (grade 0 LDD and grade 1 LDD) were also assessed using MRI scans of the lumbar spine for comparison purposes. 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. Written informed consent was obtained from the eligible patients before recruitment.

Clinical evaluation
Age, gender and body mass index (BMI) were recorded using a pretested interviewer administered questionnaire and clinical examination. Height (cm) and weight (kg) of the patients were recorded with light clothing and without shoes to the nearest 0.1 cm and 0.1 kg, respectively, and BMI was calculated (kg/m 2 ) [32]. International cut off values were used for categorisation of BMI [33].

Radiological assessment
T1 and T2 weighted sagittal MRI scans of lumbar spine were performed in supine position using a GE 1.5T MRI Scanner (Signa, General Electric, Milwaukee, Wisconsin). Patients underwent fast spin-echo sagittal T1 and T2 weighted imaging with slice thickness of 4 mm and interslice gap of 1 mm. MRI scans were assessed for the presence and severity of LDD and Modic changes from L1/L2 to L5/S1 levels by a consultant radiologist blinded to the clinical details of the patients. The severity of LDD was assessed on T2-weighted MRI images using the modified Pfirrmann grading system (Fig 2) and were graded as grade 1 (normal disc shape, no horizontal bands, distinction of nucleus and annulus is clear), grade 2 (nonhomogeneous shape with horizontal bands, some blurring between nucleus and annulus), grade 3 (nonhomogeneous shape with blurring between nucleus and annulus, shape of the annulus still recognizable), grade 4 (nonhomogeneous shape with hypointensity, distinction between nucleus and annulus impossible, disc height usually decreased), and grade 5 (same as grade 4 but with a collapsed disc space) [34,35]. A particular grade of LDD was identified for each of the lumbar levels, and scores of five lumbar levels were summed to calculate the severity of LDD.
Modic changes were recorded using both T1 and T2 weighted MRI images where vertebral end plate and bone marrow changes were scored as absent (no modic changes), type 1 (hypointense changes in T1-weighted images and hyperintense changes in T2 weighted images), type II (hyperintense changes in T1 and T2 weighted images) and type III (hypointense changes in T1 and T2 weighted images) (Figs 3-5) [4,36]. Presence of type I, type II and type III Modic changes were given a score of 1, 2 and 3, respectively. The scores of each lumbar level were summed to get the severity of Modic changes.

Selection of single nucleotide variants of candidate genes
The ACAN gene and candidate genes involved in the metabolic pathway of aggrecan were selected for the study (IL1A, IL1B, IL6, MMP3, ADAMTS4, ADAMTS5, TIMP1, TIMP2 and TIMP3). Aggrecan and molecules involved in aggrecan metabolic pathway are summarised with their functions and respective gene symbols in Table 1. The SNVs of these genes are involved in regulation of aggrecan metabolism and other regulatory functions of the extracellular matrix (ECM). Common SNVs (minor allele frequency ! 0.05) in the exonic and Grading system for assessment of MRI features of lumbar disc degeneration. Grade 1 (A)-homogeneous disc structure, normal disc height and distinction of nucleus and annulus is clear; grade 2 (B)-nonhomogeneous shape with horizontal bands, some blurring between nucleus and annulus; grade 3 (C)-nonhomogeneous shape with blurring between nucleus and annulus, shape of the annulus still recognizable; grade 4 (D)-nonhomogeneous shape with hypointensity, distinction between nucleus and annulus impossible, disc height usually decreased; grade 5 (E)-same as grade 4 but with a collapsed disc space.
doi:10.1371/journal.pone.0169835.g002 untranslated regions (UTRs) of selected genes in twenty five Sri Lankan exomes were identified using an in-house variant calling pipeline at the Human Genetics Unit, Faculty of Medicine, University of Colombo. In addition, common variants in exonic and UTRs of the selected candidate genes among Gujarati Indians in Houston (GIH) were extracted from the HapMap database [37] (HapMap database provides human genetic variations among several global ancestry groups and GIH ancestry is the genetically most similar ancestry to the Sri Lankan population). Sixty two SNVs in exonic regions and UTRs from 10 genes were selected for genotyping ( Table 2). Although there is a VNTR polymorphism and other intronic variants which are associated with the severity of LDD and Modic changes, this study was confined to exonic and UTR SNVs of the candidate genes.  Sample collection, DNA extraction and quantification Venous blood was drawn from a superficial forearm vein into an EDTA tube and stored at -20˚C at the Human Genetics Unit, Faculty of Medicine, University of Colombo until processing. DNA was extracted from 200ml of venous blood using QIAamp DNA Mini Kit according to the blood and body fluid protocol (spin protocol) [38] summarised in S1 Protocol. Extracted DNA was quantified using the Quantus fluorometer (Promega) with QuantiFluor 1 Double stranded DNA (dsDNA) system [39] summarised in S1 Protocol. Quantified DNA was diluted in water and normalised and optimized for the subsequent Sequenom MassARRAY system.

Genotyping
All selected genetic variants were genotyped using Sequenom iPLEX MassARRAY system (Sequenom, San Diego, CA) [40] at the Australian Genome Research Facility, Australia. The workflow of the Sequenom iPLEX reaction included locus specific polymerase chain reaction, locus specific primer extension reaction and annealing of oligonucleotide primer at the upstream of the polymorphic site being genotyped. Then the primer and amplified target DNA were incubated with mass modified dideoxynucleotide terminators. The primer extension was made according to the sequence of the variant site. Mass of the extended primer was determined through the use of matrix-assisted laser desorption ionization time of flight mass spectrometry. The alleles present at the polymorphic site of interest were assessed based on the sequence specified by the mass of the primer. The mass of the observed primers were translated into a genotype by the SpectroTYPER software of the Sequenom system [40].

Statistical analysis
Descriptive statistics were calculated to summarise the sample characteristics. Statistical analysis was carried out using PLINK 1.9 software [41]. Variants with Hardy-Weinberg equilibrium ! 0.05 were included into the analysis. Multiple linear regression analysis was used for assessing associations of quantitative outcomes (severity of LDD and severity of Modic changes) with SNVs of interest. The genotype of the respective SNV was treated as a

Results
Among hundred and twenty patients selected, only 106 patients attended for the MRI scanning. Sample characteristics of the 106 patients who underwent MRI scan of lumbar spine are summarised in Table 3. The majority of patients belonged to 50-59 years age category. Most of the patients were females. Twenty three of the females (31.1%) were obese while 4 males were obese (12.5%). There was moderate agreement between the x-ray scoring system and MRI scoring system of LDD (Cohen's kappa = 0.49). All 106 patients who underwent MRI assessment had atleast grade 2 LDD. Mean severity of LDD was 12.56 ± 2.88. Eighteen patients (17%) had Modic changes and type 2 Modic changes were common. Mean severity of Modic changes was 0.49 ± 1.22. Genotyping efficiency was more than 98% and minor allele frequency ranged from 0.04 to 0.5. All variants were in Hardy-Weinberg equilibrium ( Table 2).
Associations of genetic variants with the severity of lumbar disc degeneration IL1A, ADAMTS4 and ADAMTS5 genes were significantly associated with the severity of LDD (Table 4). Each additional "G" allele of the rs2856836 variant, "C" allele of the rs1304037 variant, "A" allele of the rs17561 variant and "A" allele of the rs1800587 variant of the IL1A gene were associated with progressive reduction in the severity of LDD. These four SNVs were in strong linkage disequilibrium with each other (R 2 ! 0.94, D' ! 0.98). The presence of each additional "C" allele of the rs41270041 variant of the ADAMTS4 gene was associated with progressive reduction in severity of LDD. The presence of the "A" allele of the rs226794 variant of the ADAMTS5 gene increased the severity of LDD.

Associations of genetic variants with the severity of Modic changes
IL1A, ADAMTS4, ADAMTS5 and TIMP3 genes were associated with the severity of Modic changes (Table 4). Each additional "G" allele of the rs2856836 variant, "C" allele of the rs1304037 variant, "A" allele of the rs17561 variant and "A" allele of the rs1800587 variant of the IL1A gene were associated with progressive reduction in severity of Modic changes. All four SNVs were in strong linkage disequilibrium as mentioned previously. The presence of each additional "C" allele of the rs34884997 variant of the ADAMTS4 gene was associated with progressive increase in the severity of Modic changes. Each additional "G" allele of the rs55933916 variant of the ADAMTS5 gene increased the severity of Modic changes. Furthermore, the presence of each additional "T" allele of the rs9862 variant of the TIMP3 gene was associated with progressive reduction in the severity of Modic changes. In-silico functional analysis of the significant variants associated with the severity of lumbar disc degeneration and Modic changes  (Table 5).

Discussion
In this study we assessed the associations of single nucleotide variants of the candidate genes in the metabolic pathway of the aggrecan gene with the severity of LDD and Modic changes.
In addition we performed the in-silico functional analysis of the variants significantly SNV-single nucleotide variant, LDD-lumbar disc degeneration, SD-standard deviation, A1 -minor allele, A2 -major allele, MAF-minor allele frequency, NA-not applicable, β-standardised regression coefficient, data were analysed by multiple linear regression on variant genotypes adjusting for age, gender and body mass index. *-p value < 0.05 a The tabulation of MRI features according to genotype of 62 SNVs and results of their multiple linear regression are summarised in S1 Table. doi:10.1371/journal.pone.0169835.t004 associated with the severity of LDD and Modic changes. We found that, the "G" allele of the rs2856836 variant, the "C" allele of the rs1304037 variant, the "A" allele of the rs17561 variant and the "A" allele of the rs1800587 variant of the IL1A gene to be associated with the severity of LDD and Modic changes. In addition the "C" allele of the rs41270041 variant of the ADAMTS4 gene and the "A" allele of the rs226794 variant of the ADAMTS5 gene were associated with the severity of LDD. The "C" allele of the rs34884997 variant of the ADAMTS4 gene, the "G" allele of the rs55933916 variant of the ADAMTS5 gene and the "T" allele of the rs9862 variant of the TIMP3 gene were associated with the severity of Modic changes. The rs17561 variant of the IL1A gene was predicted as pathogenic by the PolyPhen functional prediction tool.
LDD is a complex phenomenon and it was previously considered that LDD was totally age dependent. With the advancement of genetic research, genetic associations have become the main determining factor for initiation and progression of LDD [46]. Intervertebral discs with unfavourable genetic variations are weak and matrix can be easily disrupted due to activities of normal daily living [47]. Furthermore, functions of metabolic pathways (inflammatory, catabolic and anti-catabolic pathways) of structural components of the intervertebral disc can be affected by SNVs of candidate genes [7,48]. With LDD, normal balance of the biomechanical forces of the spine are disrupted and reflected on the nearby structures including the vertebrae and ligaments [3]. Modic changes are a type of degenerative changes which involve the end plate and bone marrow of the vertebrae and are well correlated with the severity of LDD [36]. In our sample, severity of Modic changes correlated positively with the severity of LDD (correlation coefficient (r) = 0.29, p < 0.01).
Aggrecan is a large aggregating proteoglycan and a main structural component of the intervertebral disc and vertebral end plate [49]. Any factor that alter the synthesis and function of the aggrecan in the extracellular matrix would make the disc vulnerable to biomechanical forces leading to increased secretion of catabolic enzymes and subsequent degeneration [50]. In our study, variants of the ACAN gene were not associated with the severity of LDD or Modic changes based on the additive genetic model. Most published studies have identified an association of VNTR polymorphism in exon 12 of the ACAN gene with LDD [11][12][13]. However, the studies which have focused on associations of SNVs of exonic regions of the ACAN gene with LDD are limited [15,16]. The "T" allele of the rs1042631 variant of the ACAN gene was associated with reduced signal intensity in a population based cohort of 588 Finnish males [16]. The Finnish study graded the LDD based on intensity of the disc signal, but the modified Pfirrmann grading system which we used includes additional variables to grade LDD such as homogeneity of the disc, distinction between nucleus and annulus and height of disc [36]. Therefore discrepancies in scoring systems may have influenced the results.
It has been postulated that persistent inflammation and subsequent vascular and nerve ingrowths are the key factors associated with painful LDD [9] and Modic changes [6]. SNVs of the IL1A gene may modify the function of the IL1α resulting altered synthesis of ECM molecules including aggrecan. According to published literature, the "T' allele of the rs1800587 variant of the IL1A gene is associated with Modic changes in middle aged men among Finnish occupational groups (machine drivers, carpenters and office workers) [18]. It is also associated with the early onset of LDD in Danish young girls (12-14 yrs) [51]. The "C" and "T" alleles are the recorded allele combination in the Finnish population, but our sample had "G" and "A" allele combination. However, the genotype frequencies are compatible with the GIH ancestry in HapMap database. According to our results, the presence of each additional "A" allele of the rs1800587 variant reduced the severity of LDD and Modic changes. The rs2856836, rs1304037, rs17561 and rs1800587 variants of the IL1A gene are in strong linkage disequilibrium and share similar minor allele frequency and functional properties. Similarly, the presence of each additional minor allele of the rs2856836, rs1304037 and rs17561 variants reduced the severity of LDD and Modic changes. These variants have been identified as associated factors in inflammatory arthritis. Specifically, the rs2856836 variant and the rs17561 variant of the IL1A gene are associated with ankylosing spondylitis [52]. The rs17561 variant and the rs1800587 variant of IL1A gene are associated with erosive/aggressive arthritis and the rs1304037 variant of the IL1A gene is associated with rheumatoid arthritis [53]. However there are no previous reports on the association between degenerative features of the spine and rs2856836, rs1304037, and rs17561 variants of the IL1A gene.
The rs2856836, rs1304037 and rs1800587 variants were located at UTRs of the IL1A gene and were predicted as non pathogenic by the mutation taster tool. The rs17561 variant is a non-synonymous variant located in exon 5 of the IL1A gene and causes an Alanine to Serine substitution at position 114 of the amino acid sequence. This substitution is a highly conservative substitution and is located at a non conserved region of the IL1A gene. Highly conservative substitutions are not considered to be pathogenic. Likewise variants located at non conserved regions of the genes are also less likely to be pathogenic. In our study, one prediction tool predicted the rs17561 variant as pathogenic (Provean-neutral; SIFT-tolerated; PolyPhen-possibly damaging; Mutation taster: polymorphism automatic). However, it involves in splice site regulation and alters the features of the final protein (IL1α) according to the Mutation Taster. As overall results are inconsistent, it is important to confirm its functional role with functional genetic studies.
ADAMTS4 and ADAMTS5 are main aggrecanases which degrade the protein core of the aggrecan molecule leading to matrix breakdown of the intervertebral disc and vertebral end plate [22,54]. According to animal studies increased expression of the ADAMTS4 gene accelerates the extracellular matrix breakdown and reduces synthesis of proteoglycan in the intervertebral disc [55,56]. Each additional "C" allele of the rs41270041 variant of the ADAMTS4 gene reduced the severity of LDD in our sample, but this SNV did not have a significant association with LDD in the Chinese Han Population [23]. The rs41270041 variant is a non-synonymous variant located in exon 9 of the ADAMTS4 gene and causes a Proline to Alanine substitution at position 720 of the amino acid sequence, but it does not affect the overall function of the ADAMTS4 gene according to in-silico functional analysis. The rs34884997 variant is a 3' UTR variant of the ADAMTS4 gene and the "C" allele was associated with increased severity of Modic changes, but there was no evidence with regard to this SNV in the literature.
The evidence for association of SNVs of exonic regions of the ADAMTS5 gene with LDD/ Modic changes is limited [24]. The presence of the "A" allele of the rs226794 variant of the ADAMTS5 gene increased the severity of LDD in our patients. This SNV has been studied for its association with Achilles tendon pathology and knee osteoarthritis, but none of the studies had strong evidence to support the hypothesis [57,58]. The rs226794 variant is a non-synonymous variant located in exon 7 of the ADAMTS5 gene and causes a Leucine to Proline substitution at position 692 of the amino acid sequence. Although it causes a non conservative substitution, in-silico functional analysis was unremarkable. In addition, the presence of each additional "G" allele of the rs55933916 variant of the ADAMTS5 gene increased the severity of Modic changes. Minor allele frequency of the "G" allele of the rs55933916 variant of the ADAMTS5 gene is 0.08, but there was only one patient with the "GG" genotype in our sample and this patient had a comparatively higher mean severity in LDD and Modic changes. The rs55933916 variant is a synonymous variant and there were no records on this SNV in the GIH ancestry in the HapMap database. Although it does not change the amino acid sequence, it is involved in splice site variation and causes subsequent changes in the protein features. However, the rs226794 and rs55933916 variants do not affect the overall function of the ADAMTS5 protein according to in-silico functional analysis. There were no previous records on both of these variants of the ADAMTS5 gene with regard to its association with LDD/ Modic change.
TIMPs are regulatory proteins which inhibit the local actions of MMP3, ADAMTS4 and ADAMTS5 molecules during the process of LDD. Although there are three TIMP molecules which are involved in LDD/Modic changes, only the "T" allele of the rs9862 variant of the TIMP3 gene was associated with a progressive reduction in the severity of Modic changes in our sample. Male oral cancer patients in Taiwan who carried the "T" allele of the rs9862 variant of the TIMP3 gene had increased plasma levels of TIMP3 molecule [59]. Therefore the "T" allele of the rs9862 variant of the TIMP3 gene may have increased the production of TIMP3 molecules resulting in increased inhibition of ADAMTS4 and ADAMTS5 molecules causing subsequent reduction of aggrecan degradation. However the rs9862 variant is a synonymous variant and results of the in-silico analysis were unremarkable.
Quantitative traits like severity of LDD and Modic changes may have cumulative/additive effects from multiple variants from several candidate genes. In addition there may be a significant effect from the environmental factors on the expression of the respective genetic variants [7]. We have not corrected the genetic association analysis for the multiple variants we already assessed and other probable variants in other respective genes. Only one of the significant variants predicted as pathogenic by the PolyPhen prediction tool and results of the other prediction tools were inconsistent. Therefore studies using a larger sample number with whole genome sequencing may be required to identify the significant SNVs associated with these quantitative traits and will need to be followed by functional genetic studies. However this study provides an initial basis for identifying probable SNVs which are associated with the severity of LDD and Modic changes.

Limitations of the study
This study was conducted on a specific group of patients with chronic mechanical low back pain and the results are not generalisable. Although the study was conducted at one centre, it represented patients from three districts of the Western Province of Sri Lanka. Low sample number is another limitation and the study had a 80% power to detect significant p values of regression coefficients ! 0.25 according to "Quanto" software programme [60]. However, with the use of permutation function in PLINK software tool and other covariates (age, gender, and BMI), the study detected significant p values of regression coefficients ! 0.2. There may be other confounding factors for LDD and Modic changes such as occupation, type of physical activity which we have not corrected for.

Conclusions
Variants of the IL1A, ADAMTS4 and ADAMTS5 genes are associated with the severity of LDD and Modic changes. In addition, a variant of the TIMP3 gene is associated with the severity of Modic changes. The degenerative findings of the spine are not a mere product of ageing process. Although nine SNVs have significant associations with LDD and Modic changes, only the rs17561 variant of the IL1A gene is predicted as pathogenic by in-silico functional analysis. Generally many of the other SNVs are acting as genetic markers for probable functional variants in the locations elsewhere in the same gene or nearby genes.