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-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-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.

Methods

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.

Fig 1

Assessment of the x-ray features of lumbar disc degeneration—lateral x-ray of lumbar spine.

Arrows—A—no disc space narrowing / anterior osteophytes (grade 0 LDD); B—mild disc space narrowing and small anterior osteophytes (grade 1 LDD); C—small anterior osteophytes and moderate disc space narrowing (grade 2 LDD).

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/m2) [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.

Fig 2

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.

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.

Fig 3

Type I Modic changes in lumbar spine.

Hypointense changes in T1-weighted images (A) and hyperintense changes in T2 weighted MRI images (B).

Fig 5

Type III Modic changes in lumbar spine.

Hypointense changes in T1 (A) and T2 (B) weighted images.

Type II Modic changes in lumbar spine.

Hyperintense changes in both T1 (A) and T2 (B) weighted images.

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 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.

Table 1

Selected genes which are involved in aggrecan and its metabolic pathway.

Component	Molecule	Gene symbol	Function
Structural	Aggrecan	ACAN	Main proteoglycan of the extracellular matrix and provides osmotic properties with the help of glycosaminoglycan side chains which absorbs and retains water
Inflammatory	Interleukin 1α	IL1A	Reduces the synthesis of collagen type I, II and aggrecan
	Interleukin 1β	IL1B	Inhibits nucleus pulposus cell proliferation and increases the synthesis of catabolic enzymes of aggrecan
	Interleukin 6	IL6	Protective role in non degenerated cells—controls the actions of IL1-β and TNF-α; Inflammatory role in degenerated cells—increases the synthesis of IL1-β and TNF-α
Catabolic	Matrix metalloproteinases 3	MMP3	Digests protein core of the aggrecan and small non collagenous proteins
	A disintegrin and metalloproteinase with thrombospondin motifs 4	ADAMTS4	Major aggrecanase which is involved in degradation of aggrecan core protein
	A disintegrin and metalloproteinase with thrombospondin motifs 5	ADAMTS5	Major aggrecanase which is involved in degradation of aggrecan core protein
Anti catabolic	Tissue inhibitors of metalloproteinase 1	TIMP1	Inhibits the catabolic actions of MMPs and ADAMTS
	Tissue inhibitors of metalloproteinase 2	TIMP2	Inhibits the catabolic actions of MMPs and ADAMTS
	Tissue inhibitors of metalloproteinase 3	TIMP3	Inhibits the catabolic actions of MMPs and ADAMTS

Table 2

Allele frequencies and Hardy–Weinberg equilibrium of genetic variants of candidate genes in aggrecan metabolic pathway.

No.	CHR	Position GRCh37/hg19	Gene	Region/function	SNV	Minor allele	Major allele	MAF (%)	HWE (p value)
1	1	161160872	ADAMTS4	3’ UTR	rs34884997	C	T	0.12	0.68
2	1	161161284	ADAMTS4	nonsynonymous	rs41270041	C	G	0.18	0.52
3	1	161163037	ADAMTS4	nonsynonymous	rs4233367	T	C	0.18	0.76
4	1	161168004	ADAMTS4	synonymous	rs33941127	T	C	0.40	0.85
5	1	161168189	ADAMTS4	nonsynonymous	rs34448954	T	C	0.13	0.42
6	2	113532083	IL1A	3’ UTR	rs2856836	G	A	0.27	0.06
7	2	113532236	IL1A	3’ UTR	rs1304037	C	T	0.28	0.06
8	2	113537223	IL1A	nonsynonymous	rs17561	A	C	0.27	0.06
9	2	113542960	IL1A	5’ UTR	rs1800587	A	G	0.27	0.06
10	2	113587121	IL1B	downstream	rs2853550	A	G	0.29	0.65
11	2	113590390	IL1B	synonymous	rs1143634	A	G	0.13	0.69
12	7	22766246	IL6	upstream	rs1800796	C	G	0.46	0.59
13	7	22766645	IL6	upstream	rs1800795	C	G	0.13	0.10
14	7	22771156	IL6	synonymous	rs2069849	T	C	0.05	1.00
15	11	102709425	MMP3	synonymous	rs520540	A	G	0.33	0.10
16	11	102713465	MMP3	synonymous	rs602128	A	G	0.33	0.10
17	11	102713620	MMP3	nonsynonymous	rs679620	T	C	0.33	0.21
18	15	89382027	ACAN	synonymous	rs372041880	C	A	0.07	0.10
19	15	89382129	ACAN	nonsynonymous	rs16942318	A	C	0.05	0.05
20	15	89386652	ACAN	nonsynonymous	rs34949187	A	G	0.11	1.00
21	15	89388894	ACAN	nonsynonymous	rs148070768	G	A	0.14	0.26
22	15	89388905	ACAN	synonymous	rs16942341	T	C	0.07	1.00
23	15	89391160	ACAN	synonymous	rs2272023	C	A	0.34	0.42
24	15	89392689	ACAN	nonsynonymous	rs144501729	A	C	0.14	0.26
25	15	89398105	ACAN	synonymous	rs2351491	T	C	0.37	0.84
26	15	89398553	ACAN	nonsynonymous	rs35430524	A	C	0.05	0.05
27	15	89398553	ACAN	nonsynonymous	rs3743399	G	A	0.28	0.65
28	15	89398631	ACAN	nonsynonymous	rs938609	A	T	0.37	0.70
29	15	89400339	ACAN	nonsynonymous	rs2882676	A	C	0.47	0.46
30	15	89400680	ACAN	nonsynonymous	rs28407189	G	A	0.07	1.00
31	15	89400963	ACAN	nonsynonymous	rs79925540	T	G	0.16	0.48
32	15	89401109	ACAN	nonsynonymous	rs4932439	A	G	0.28	0.65
33	15	89401615	ACAN	synonymous	rs3825994	G	T	0.44	0.85
34	15	89401616	ACAN	nonsynonymous	rs76282091	C	G	0.16	0.48
35	15	89402051	ACAN	nonsynonymous	rs1042630	A	G	0.41	0.45
36	15	89402239	ACAN	synonymous	rs1042631	T	C	0.35	1.00
37	15	89402596	ACAN	synonymous	rs698621	G	T	0.38	0.85
38	15	89415247	ACAN	nonsynonymous	rs3817428	G	C	0.05	1.00
39	15	89417238	ACAN	nonsynonymous	rs1126823	G	A	0.36	0.84
40	17	76867017	TIMP2	synonymous	rs2277698	T	C	0.24	0.80
41	21	28291455	ADAMTS5	3’ UTR	rs1444269	G	A	0.26	1.00
42	21	28291846	ADAMTS5	3’ UTR	rs2298657	C	T	0.05	1.00
43	21	28292581	ADAMTS5	3’ UTR	rs3746836	A	G	0.25	0.81
44	21	28293095	ADAMTS5	3’ UTR	rs229072	T	A	0.48	0.27
45	21	28293117	ADAMTS5	3’ UTR	rs229073	G	A	0.48	0.27
46	21	28293924	ADAMTS5	3’ UTR	rs11700721	T	C	0.12	1.00
47	21	28294090	ADAMTS5	3’ UTR	rs16979423	G	T	0.14	0.70
48	21	28294143	ADAMTS5	3’ UTR	rs9978597	G	T	0.04	1.00
49	21	28296135	ADAMTS5	3’ UTR	rs229078	T	G	0.22	1.00
50	21	28296324	ADAMTS5	3’ UTR	rs151065	A	G	0.19	0.56
51	21	28296389	ADAMTS5	synonymous	rs3746839	G	A	0.08	1.00
52	21	28302355	ADAMTS5	nonsynonymous	rs226794	A	G	0.10	0.60
53	21	28338298	ADAMTS5	nonsynonymous	rs457947	G	C	0.07	0.46
54	21	28338423	ADAMTS5	synonymous	rs55933916	G	C	0.08	0.50
55	22	33253280	TIMP3	synonymous	rs9862	T	C	0.50	0.46
56	22	33253292	TIMP3	synonymous	rs11547635	T	C	0.07	0.42
57	22	33257322	TIMP3	3’ UTR	rs1427384	C	T	0.19	0.53
58	22	33258050	TIMP3	3’ UTR	rs2267184	T	C	0.16	0.48
59	22	33258288	TIMP3	3’ UTR	rs1065314	C	T	0.17	0.32
60	23	47444879	TIMP1	nonsynonymous	rs5953060	C	G	0.43	0.66
61	23	47444985	TIMP1	synonymous	rs4898	C	T	0.43	0.66
62	23	47445286	TIMP1	nonsynonymous	rs6609533	G	A	0.43	0.66

CHR—chromosome number, SNV—single nucleotide variant, MAF—minor allele frequency, HWE—Hardy Weinberg equilibrium, UTR—untranslated region

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® 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 quantitative variable and coded 0, 1 or 2 to represent the number of variant allele, consistent with an additive genetic model. Age, gender and BMI were used as the additional covariates. Permutation tests with 10,000 permutations were carried out to generate the significance levels empirically 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. In-silico functional analysis was carried out using four free online bioinformatics prediction tools: Provean [42], SIFT [43], PolyPhen [44] and Mutation taster [45].

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).

Table 3

Summary of the sample characteristics of 106 patients who underwent genetic association analysis.

Variable		N (%)
Total patients		106
Age	Mean	52.42 ± 9.42 yrs
	20–29 yrs	3 (2.8)
	30–39 yrs	7 (6.6)
	40–49 yrs	22 (20.8)
	50–59 yrs	51 (48.1)
	60–69 yrs	23 (21.7)
Gender	Female	74 (69.8)
	Male	32 (30.2)
BMI	Normal (18–24.9 kg/m2)	37 (34.9)
	Overweight (25–29.9 kg/m2)	42 (39.6)
	Obese (≥ 30 kg/m2)	27 (25.5)
LDD (maximum score)	Grade 1	0 (0.0)
	Grade 2	7 (6.6)
	Grade 3	56 (52.8)
	Grade 4	36 (34.0)
	Grade 5	7 (6.6)
Modic changes	Yes	18 (17.0)
	No	88 (83.0)
Types of Modic changes	I	1 (0.9)
	II	15 (14.2)
	III	2 (1.9)

LDD—lumbar disc degeneration, LDD (maximum score)–maximum grade of LDD in L1/L2 to L5/S1 levels of the respective spine, BMI—body mass index

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 (R2 ≥ 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.

Table 4

Severity of MRI features tabulated according to genotype and the results of multiple linear regression (significant SNVs).

No.	Gene	SNV	Variable	A1	MAF	A2/A2	A1/A2	A1/A1	Regression coefficient (β)	Adjusted p value
1	ADAMTS4	rs34884997		C	0.12
			Genotypes			T/T	C/T	C/C
			N			80	23	2
			LDD mean (SD)			12.34 (2.97)	13.13 (2.56)	14.00 (2.83)	0.12	0.14
			Modic mean (SD)			0.30 (1.05)	1.04 (1.58)	2.00 (0.00)	0.31	0.00*
2	ADAMTS4	rs41270041		C	0.18
			Genotypes			G/G	C/G	C/C
			N			69	35	1
			LDD mean (SD)			12.72 (2.73)	12.34 (3.10)	7.00 (0.00)	-0.22	0.01*
			Modic mean (SD)			0.51 (1.21)	0.49 (1.29)	0.00 (0.00)	-0.03	0.80
3	IL1A	rs2856836		G	0.27
			Genotypes			A/A	G/A	G/G
			N			58	38	9
			LDD mean (SD)			13.09 (2.93)	11.92 (2.72)	11.67 (2.83)	-0.21	0.01*
			Modic mean (SD)			0.69 (1.46)	0.32 (0.87)	0.00 (0.00)	-0.20	0.04*
4	IL1A	rs1304037		C	0.28
			Genotypes			T/T	C/T	C/C
			N			58	37	10
			LDD mean (SD)			13.09 (2.93)	11.97 (2.73)	11.50 (2.72)	-0.22	0.01*
			Modic mean (SD)			0.69 (1.46)	0.32 (0.88)	0.00 (0.00)	-0.20	0.04*
5	IL1A	rs17561		A	0.27
			Genotypes			C/C	A/C	A/A
			N			58	38	9
			LDD mean (SD)			13.09 (2.93)	11.92 (2.72)	11.67 (2.83)	-0.21	0.01*
			Modic mean (SD)			0.69 (1.46)	0.32 (0.87)	0.00 (0.00)	-0.20	0.04*
6	IL1A	rs1800587		A	0.27
			Genotypes			G/G	A/G	A/A
			N			59	38	9
			LDD mean (SD)			13.10 (2.91)	12.11 (2.82)	10.89 (2.03)	-0.23	0.00*
			Modic mean (SD)			0.68 (1.46)	0.32 (0.87)	0.00 (0.00)	-0.20	0.04*
7	ADAMTS5	rs226794		A	0.1
			Genotypes			G/G	A/G	A/A
			N			84	21	0
			LDD mean (SD)			12.14 (2.92)	14.14 (2.15)	NA	0.20	0.02*
			Modic mean (SD)			0.48 (1.26)	0.57 (1.12)	NA	0.03	0.75
8	ADAMTS5	rs55933916		G	0.08
			Genotypes			C/C	G/C	G/G
			N			90	13	1
			LDD mean (SD)			12.56 (2.92)	12.31 (3.87)	15.00 (0.00)	-0.02	0.82
			Modic mean (SD)			0.40 (1.08)	0.62 (1.26)	6.00 (0.00)	0.27	0.01*
9	TIMP3	rs9862		T	0.5
			Genotypes			C/C	T/C	T/T
			N			23	58	23
			LDD mean (SD)			13.26 (2.42)	12.64 (3.06)	11.61 (2.78)	-0.15	0.07
			Modic mean (SD)			1.04 (1.58)	0.40 (1.15)	0.22 (0.85)	-0.23	0.02*

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.

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

Among twenty four nonsynonymous variants, only three were associated with severity of LDD/Modic changes. Amino acid substitutions of the rs41270041 variant of the ADAMTS4 gene and the rs226794 variant of the ADAMTS5 gene are non conservative substitutions. Among the nine significant variants, only the rs17561 variant of the IL1A gene is predicted as pathogenic by the PolyPhen functional prediction tool. The rs55933916 variant of the ADAMTS5 gene and the rs9862 variant of the TIMP3 gene are located in conserved regions of the respective genes. However both these variants were predicted as non pathogenic by the four prediction tools (Table 5).

Table 5

Significant SNVs associated with severity of lumbar disc degeneration/Modic changes and their in-silico functional analysis.

No	Gene	SNV	Region/function	Significant associations	Amino acid change	Conservative substitution	Conserved region	Provean	SIFT	PolyPhen	Mutant Taster
1	ADAMTS4	rs34884997	3' UTR	Modic changes			No				Harmless
2		rs41270041	nonsynonymous	LDD	P720A	No	No	Neutral	Tolerated	Benign	Harmless
3	IL1A	rs2856836	3' UTR	LDD, Modic changes			No				Harmless
4		rs1304037	3' UTR	LDD, Modic changes			No				Harmless
5		rs17561	nonsynonymous	LDD, Modic changes	A114S	Yes	No	Neutral	Tolerated	Damaging	Harmless
6		rs1800587	5' UTR	LDD, Modic changes
7	ADAMTS5	rs226794	nonsynonymous	LDD	L692P	No	No	Neutral	Tolerated	Benign	Harmless
8		rs55933916	synonymous	Modic changes			Yes	Neutral	Tolerated	Benign	Harmless
9	TIMP3	rs9862	synonymous	Modic changes			Yes	Neutral	Tolerated	Benign	Harmless

SNV—single nucleotide variant, LDD—lumbar intervertebral disc degeneration, UTR—untranslated regions, Mutant taster (polymorphism automatic) = Harmless, AA—amino acid

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 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-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.

DNA extraction and quantification methodology.

(DOCX)

Click here for additional data file.

The tabulation of MRI features according to genotype of 62 SNVs and results of their multiple linear regression.

(DOCX)

Click here for additional data file.

Genotype data of 106 patients.

(XLSX)

Click here for additional data file.