Novel Genetic Variants Associated with Lumbar Spondylosis in Koreans : A Genome-Wide Association Study.

OBJECTIVE: The aim of this study was to identify the susceptibility genes responsible for lumbar spondylosis (LS) in Korean patients.
METHODS: Data from 1427 subjects were made available for radiographic grading and genome wide association studies (GWAS) analysis. Lateral lumbar spine radiographs were obtained and the various degrees of degenerative change were semi-quantitatively scored. A pilot GWAS was performed using the AffymetrixGenome-Wide Human single-nucleotide polymorphisms (SNPs), 500K array. A total of 352228 SNPs were analyzed and the association between the SNPs and case-control status was analyzed by stepwise logistic regression analyses.
RESULTS: The top 100 SNPs with a cutoff p-value of less than 3.7×10-4 were selected for joint space narrowing, while a cutoff p-value of 6.0×10-4 was applied to osteophytes and the Kellgren-Lawrence (K-L) osteoarthritis grade. The SNPs with the strongest effect on disc space narrowing, osteophytes, and K-L grade were serine incorporator 1 (rs155467, odds ratio [OR]=17.58, p=1.6×10-4), stromal interaction molecule 2 (STIM1, rs210781, OR=5.53, p=5×10-4), and transient receptor potential cation channel, subfamily C (rs11224760, OR=3.99, p=4.8×10-4), respectively. Leucine-rich repeat-containing G protein-coupled receptor 4 was significantly associated with both disc space narrowing and osteophytes (rs1979400, OR=2.01, p=1.1×10-4 for disc space narrowing, OR=1.79, p=3×10-4 for osteophytes), while zinc finger and BTB domain containing 7C was significantly and negatively associated with both osteophytes and a K-L grade >2 (rs12457004,OR=0.25, p=5.8×10-4 and OR=0.27, p=5.3×10-4, respectively).
CONCLUSION: We identified SNPs that potentially contribute to the pathogenesis of LS. This is the first report of a GWAS in an Asian population.

INTRODUCTION

Low back pain (LBP) is an important public health problem and is associated with substantial societal costs in industrialized countries. Degenerative disease of the lumbar spine (lumbar spondylosis, LS), which is characterized radiologically by the presence of osteophytes, endplate sclerosis, and disc space narrowing13), is thought to be related to LBP. The prevalence of LS in Asian countries has scarcely been reported. However, ethnic differences in the prevalence of osteophytosis and disc degeneration have previously been noted22). Although it is widely known that the prevalence of LS increases with age, there exists a paucity of data on the risk factors for LS. LS shares many features with osteoarthritis (OA) of the peripheral joint, the most common form of arthritis affecting the elderly. Although the pathogenesis of OA has been associated with aging and environmental factors such as work activity, twin studies have revealed a significant genetic contribution, with a heritability estimate of 65% for hand OA20). Recent reports have demonstrated the significant role that genetic contribution also plays in LS. A twin study showed that the heritability of LS was 74% for the lumbar spine and 73% for the cervical spine2). Other studies also revealed significant heritability on the basis of magnetic resonance imaging (MRI) of the spine9,16). The search for specific genes associated with OA and LS has since ensued. However, due to variations in the phenotype, as well as the severity of the disease and patient characteristics such as gender and age, the definitions of the disease have varied from study to study, leading to tremendous difficulty in identifying a causative genetic factor.

A large number of genetic association studies of OA have been conducted over the last decade, and a number of significant associations have been verified by a systematic review. For example, aggrecan 1 (AGC1) is associated with hand OA, asporin (ASPN) with hip and knee OA, and vitamin D receptor (VDR) with knee and spine OA15). However, the majority of the genes that were reported to have significant associations have never been replicated, making it difficult to discern the significance of these isolated significant associations in the pathogenesis of OA.

Complex conditions such as OA and LS are caused by numerous genetic and environmental factors, any one of which can have a relatively minor effect. Genome-wide association studies (GWAS) have provided a breakthrough in deciphering the genetic influence in complex diseases, by examining hundreds of thousands of single-nucleotide polymorphisms (SNPs), thus enabling identification of SNPs implicated in a large number of robustly replicated loci of common traits.

Although GWAS has begun to unravel the genetic influence of peripheral joint OA25), there has been little attempt as yet to do the same for LS. In addition, the results of an LS GWAS among an Asian population have not yet been reported. In this study, we sought to identify the susceptibility genes for LS in a GWAS among a community-based population of Korean adults.

MATERIALS AND METHODS

Study population

For this study, we selected a rural farming community (Ansung) in South Korea from among the populations incorporated into the ongoing prospective Korean Genome and Epidemiology Study (KoGES). The study methods have been described previously4). Briefly, the eligibility criteria included an age of 40–79 years, residence within the borders of the survey area for at least 6 months before testing, and the mental and physical ability to participate. Data relating to LBP were available for 4181 subjects for the years 2006 and 2007. Baseline demographic information was collected using a standard questionnaire during a face-to-face interview and included educational attainment, occupation, exercise, and co-morbidities.

Of the total number of eligible subjects, 2000 were randomly selected for spine radiography. No difference was found in the prevalence of LBP between those who underwent radiography and those who did not. After excluding 488 subjects who were unable to be evaluated due to poor film quality, and 85 patients whose genomic DNA was not obtained, data for 1427 subjects were available for radiograph grading and GWAS analysis. The study protocol was approved by the Ethics Committee of the KoGES, and written informed consent was obtained from each participant.

Radiographic evaluation of the lumbar spine

Lateral lumbar spine radiographs were taken according to a standard protocol with the film centered on the second lumbar vertebra. Lumbar radiographs were evaluated by a single observer. Each vertebral level from L1/2 to L4/5 was reviewed for the presence of radiographic features relating to degenerative change. Semi-quantitative scores were given for the following features using a reference atlas8) : presence and severity of anterior osteophytes (grade 0=none; grade 1=barely visible; grade 2=definite; grade 3=large), endplate sclerosis (grade 0=none; grade 1=present), and disc space narrowing (grade 0=none; grade 1=barely visible; grade 2=definite; grade 3=severe, bone to bone).

Additionally, the Kellgren-Lawrence (K-L) grading system was used for each vertebral level (grade 0=normal disc with no osteophytes; grade 1=slight anterior wear and osteophyte formation; grade 2=definite anterior wear and mild disc space narrowing with osteophyte formation; grade 3=moderate disc space narrowing with osteophytes and sclerosis; grade 4=large osteophytes, marked disc space narrowing, and sclerosis of vertebral end plates). LS was defined using one of the following criteria : 1) joint space narrowing ≥grade 2; 2) osteophytes ≥grade 2; 3) K-L grade ≥2. Radiographs were read by a single reader who was an academically-based rheumatologist. Intra-observer reproducibility was assessed by re-evaluating 50 films within 1 week of the first reading. The reproducibility of intra-reader assessments was high (for osteophyte grading, κ=0.89–0.93; for endplate sclerosis, κ=0.71–84; for joint space narrowing, κ=0.81–0.89; and for K-L grading, κ=0.69–0.80, for various vertebral levels). Films allocated different grades at each of the two readings were adjudicated by consensus between the original reader and a second reader.

Genome-wide association study

Genomic DNA was isolated from peripheral blood mononuclear cells. We performed a pilot GWAS, typing cases and controls on a single platform using the Affymetrix Genome-Wide Human SNP 500K array chip (Affymetrix, Inc., Santa Clara, CA, USA). Genotype calls were determined by Bayesian robust linear modeling using the Mahalanobis distance algorithm14). We sequentially discarded 38364 markers with a Hardy-Weinberg equilibrium p-value <10−6, 17926 with genotype call rates below 95%, and 92050 with a minor allele frequency (MAF) of 0.01. This left 352228 SNPs available for subsequent analysis

Statistical analysis

The evaluation of the association between case-control status and each individual SNP was based on the odds ratio (OR) and p-values. Covariates used for multivariable-adjustment were age, sex, education, and body mass index. To optimize the joint effect of the SNPs, we conducted stepwise logistic regression analyses and also analyzed the main effects on the selected variables using a generalized linear model implemented with the R statistical software (ver. 2.14.2; R Core Development Team, Austria, Vienna).

RESULTS

Table 1 displays the clinical characteristics of the study population, with cases of spine OA defined as disc space narrowing, osteophytes, or K-L grade all ≥2. Subjects with spine OA were older, more likely to be male, and had a lower level of education compared to the control subjects.

Table 1

Demographic characteristics of the study subjects

	Disc space narrowing			Osteophyte			K-L grade
	Case	Control	Total	Case	Control	Total	Case	Control	Total
Age (years)	64.7±0.3*	58.5±0.4	59.3±0.2	62.4±0.3*	56.3±0.3	59.5±0.2	61.2±0.3*	55.3±0.4	59.5±0.2
Female	132 (57.9)	699 (58.4)	831 (58.3)	369 (48.4)*	462 (69.6)	831 (58.2)	536 (52.1)*	293 (74.2)	829 (58.2)
BMI	24.1±0.2	24.4±0.1	24.3±0.1	24.3±0.1	24.3±0.1	24.3±0.1	24.3±0.1	24.2±0.2	24.3±0.1
Smoker	45 (19.7)	222 (18.6)	267 (18.7)	168 (22.1)*	100 (15.1)	268 (18.8)	215 (20.9)*	52 (13.2)	267 (18.8)
Alcohol	79 (34.7)*	505 (42.2)	584 (41.0)	323 (42.3)	261 (39.3)	584 (40.9)	429 (41.7)	155 (39.2)	584 (41.0)
Employed	216 (94.7)	1124 (93.8)	1340 (94.0)	723 (94.8)	617 (92.2)	1340 (93.9)	974 (94.7)	364 (92.2)	1338 (94.0)
Education ≤6 years	64 (28.1)*	196 (16.4)	260 (18.2)	159 (20.8)*	101 (15.2)	260 (18.2)	207 (20.1)*	51 (12.9)	258 (18.1)

Values are presented as mean±standard deviation or number (%). Obesity=BMI≥27.

Denotes significant difference (p<0.05) compared to control.

K-L : kellgren-lawrence, BMI : body mass index

Genome-wide association

We used logistic regression analysis to identify statistically significant associations with LS. For multiple comparisons, a Bonferroni-corrected p-value of 1.5×10−7 (0.05/320942) was used. Since none of the SNPs reached the very conservative Bonferroni corrected value, the top 100 SNPs with a cutoff p-value less than 3.7×10−4 were selected for use in this GWAS for joint space narrowing while a cutoff p-value of 6.0×10−4 was applied for osteophytes and K-L grade. Tables 2–4 present the detailed characteristics of the SNPs, including gene name, rs identification number, position, MAF of the cases and controls, OR, and p-value. Each feature of LS was found to be associated with distinct SNPs. The SNPs that had the strongest effect on disc space narrowing included serine incorporator 1 (rs155467, OR=17.58, p=1.6×10−4), heat shock transcription factor 2 (rs563084, OR=14.48, p=1.4×10−4), cysteine-rich hydrophobic domain (rs1568512, OR=11.57, p=4×10−5), akirin 2 (rs2787938, OR=9.86, p=8×10−5), and fibroblast growth factor receptor 2 (rs11200052, OR=9.72, p=9×10−5), while those that had the most significant effect on osteophytes included stromal interaction molecule 2 (STIM1, rs210781, OR=5.53, p=5×10−4), protein kinase C and casein kinase substrate in neurons 2 (PACSIN2, rs738379, OR=5.37, p=3.6×10−4), and ubiquinol-cytochrome c reductase complex chaperone (rs6060373, OR=3.05, p=8×10−5). Transient receptor potential cation channel, subfamily C (rs11224760, OR=3.99, p=4.8×10−4), hypocretin (orexin) receptor 1 (rs3753613, OR=2.62, p=2×10−4), and coiled-coil domain containing 60 (rs10849640, OR=2.59, p=3.8×10−4) had the strongest effect on K-L grades >2. Leucine-rich repeat-containing G protein-coupled receptor 4 was significantly associated with both disc space narrowing and osteophytes (rs1979400, OR=2.01, p=1.1×10−4 for disc space narrowing, OR=1.79, p=3×10−4 for osteophytes), while zinc finger and BTB domain containing 7C was significantly and negatively associated with both osteophytes and a K-L grade >2 (rs12457004, OR=0.25, p=5.8×10−4 for osteophytes, OR=0.27, p=5.3×10−4 for K-L grade). No SNP was found to be significantly associated with LS defined using all three criteria.

Table 2

Results of the genome-wide association study showing top 30 SNPs having the highest OR associated with ≥2 K-L grade

rs number	Chr	Gene	MAF	OR (95%CI)	Bonf p-value
rs11224760	11q22.1	NA	0.2299	3.99 (1.835–8.674)	0.00048
rs3753613	1p35.2	HCRTR1	0.2949	2.618 (1.577–4.347)	0.0002
rs10849640	12q24.23	NA	0.301	2.595 (1.534–4.39)	0.00038
rs11120305	1q41	PTPN14	0.3232	2.565 (1.552–4.24)	0.00024
rs2271933	1p35.2	HCRTR1	0.2824	2.476 (1.476–4.152)	0.00059
rs6569814	6q23.2	TAAR2 TAAR3	0.3197	2.385 (1.488–3.821)	0.0003
rs1164894	9q34.3	NA	0.3685	2.335 (1.522–3.582)	0.0001
rs10774756	12q24.21	LOC105369998	0.3997	2.212 (1.483–3.301)	0.0001
rs1473047	5q14.3	NA	0.3982	2.173 (1.46–3.234)	0.00013
rs10072084	5q14.3	NA	0.3966	2.135 (1.434–3.178)	0.00019
rs7966636	12q24.21	LOC105369998	0.3749	2.092 (1.374–3.183)	0.00057
rs197457	6q24.2	HIVEP2	0.3891	2.085 (1.403–3.1)	0.00028
rs3794214	12q24.31	ACADS	0.3825	2.039 (1.367–3.043)	0.00048
rs6868338	5q34	NA	0.392	1.984 (1.349–2.918)	0.0005
rs2834443	21q22.11	NA	0.4292	1.969 (1.368–2.835)	0.00027
rs9533738	13q14.11	LOC105370182	0.4776	1.89 (1.36–2.627)	0.00015
rs2378931	14q12	LOC105370438	0.4329	1.876 (1.326–2.654)	0.00038
rs10026693	4p16.1	SORCS2	0.4699	1.835 (1.325–2.542)	0.00026
rs2878620	4p16.1	SORCS2	0.4703	1.834 (1.324–2.54)	0.00027
rs3857194	4p16.1	SORCS2	0.471	1.794 (1.297–2.481)	0.00041
rs9884489	4p16.1	SORCS2	0.4706	1.794 (1.297–2.481)	0.00041
rs6833329	4p16.1	SORCS2	0.4706	1.794 (1.297–2.481)	0.00041
rs2937545	5p13.2	NA	0.4675	1.767 (1.28–2.441)	0.00055
rs873471	8p22	PSD3	0.488	0.6036 (0.4528–0.8047)	0.00058
rs2063076	8p22	PSD3	0.4941	0.6007 (0.451–0.8001)	0.00049
rs4667789	2q24.3	SCN2A SCN3A	0.4888	0.6002 (0.4491–0.8022)	0.00056
rs751217	8p22	PSD3	0.4941	0.599 (0.45–0.7974)	0.00045
rs7556825	2q24.3	SCN2A SCN3A	0.4968	0.5984 (0.4516–0.7929)	0.00035
rs2914908	5q23.1	NA	0.478	0.5889 (0.4398–0.7886)	0.00038
rs7576705	2q37.3	NA	0.486	0.5839 (0.4373–0.7796)	0.00026

SNPs : single-nucleotide polymorphisms, OR : odds ratio, K-L : kellgren-lawrence, CI : confidence interval, MAF : minor allele frequency, NA : not available

Table 3

Results of the genome-wide association study showing top 30 SNPs having the highest OR associated with ≥2 osteophyte

rs number	Chr	Gene	MAF	OR (95%CI)	Bonf p-value
rs210781	4p15.2	NA	0.1943	5.533 (2.107–14.53)	1.00E+00
rs738379	22q13.2	PACSIN2	0.1738	5.368 (2.133–13.51)	1.00E+00
rs2284097	22q13.2	PACSIN2	0.1737	5.358 (2.129–13.48)	1.00E+00
rs2038062	22q13.2	PACSIN2	0.1761	4.532 (1.901–10.8)	1.00E+00
rs6060373	20q11.22	UQCC1	0.2566	3.053 (1.756–5.306)	1.00E+00
rs6088791	20q11.22	UQCC1	0.2526	3.053 (1.754–5.312)	1.00E+00
rs1539581	1p13.2	NA	0.211	2.991 (1.653–5.411)	1.00E+00
rs6060369	20q11.22	UQCC1	0.2523	2.989 (1.714–5.211)	1.00E+00
rs2425062	20q11.22	UQCC1	0.2565	2.988 (1.716–5.204)	1.00E+00
rs8127664	21q22.3	NA	0.2411	2.852 (1.664–4.89)	1.00E+00
rs761166	22q13.31	PARVB	0.2511	2.707 (1.557–4.705)	1.00E+00
rs10998893	10q22.1	NA	0.2519	2.648 (1.52–4.612)	1.00E+00
rs17705721	4q34.3	NA	0.2458	2.638 (1.557–4.467)	1.00E+00
rs4911178	20q11.22	UQCC1	0.2573	2.549 (1.501–4.33)	1.00E+00
rs4911496	20q11.22	UQCC1	0.2575	2.544 (1.498–4.322)	1.00E+00
rs1570004	20q11.22	UQCC1	0.2568	2.542 (1.496–4.317)	1.00E+00
rs10503404	8p23.1	MSRA	0.2913	2.47 (1.546–3.946)	1.00E+00
rs12038162	1q44	SMYD3	0.2435	2.451 (1.513–3.97)	1.00E+00
rs761165	22q13.31	PARVB	0.2776	2.429 (1.491–3.957)	1.00E+00
rs930140	2q36.1	PAX3	0.3026	2.277 (1.466–3.535)	1.00E+00
rs11970088	6p22.3	NA	0.2904	2.262 (1.454–3.52)	1.00E+00
rs10958163	8q21.13	NA	0.308	2.198 (1.442–3.35)	1.00E+00
rs11578091	1q41	NA	0.3061	2.185 (1.426–3.347)	1.00E+00
rs10848193	12q24.33	NA	0.2966	2.128 (1.418–3.195)	1.00E+00
rs16934897	10p11.22	NA	0.3385	2.028 (1.373–2.995)	1.00E+00
rs10194645	2p25.3	NA	0.4203	1.932 (1.412–2.642)	1.00E+00
rs1380255	1q41	NA	0.3742	1.908 (1.342–2.711)	1.00E+00
rs2049164	2q34	ERBB4	0.3791	1.898 (1.359–2.651)	1.00E+00
rs12564579	1p36.21	KAZN	0.4042	1.85 (1.352–2.531)	1.00E+00
rs9360980	6q14.1	NA	0.4209	1.844 (1.346–2.526)	1.00E+00

SNPs : single-nucleotide polymorphisms, OR : odds ratio, MAF : minor allele frequency, CI : confidence interval, NA : not available

Table 4

Results of the genome-wide association study showing top 30 SNPs having the highest OR associated with ≥2 disc space narrowing

rs number	Chr	Gene	MAF	OR (95% CI)	Bonf p-value
rs155467	6q22.31	PKIB	0.08182	17.58 (3.979–77.71)	1.00E+00
rs155458	6q22.31	PKIB	0.08117	17.54 (3.97–77.47)	1.00E+00
rs563084	6q22.31	NA	0.08351	14.48 (3.653–57.41)	1.00E+00
rs1568512	4q12	CHIC2	0.08514	11.57 (3.61–37.08)	1.00E+00
rs2787938	6q15	NA	0.08875	9.858 (3.176–30.6)	1.00E+00
rs2754273	6q15	NA	0.08875	9.858 (3.176–30.6)	1.00E+00
rs6915593	6q15	NA	0.09005	9.826 (3.165–30.51)	1.00E+00
rs11200052	10q26.13	NA	0.1031	9.717 (3.126–30.21)	1.00E+00
rs3870374	8q24.13	NA	0.1226	7.942 (3.087–20.44)	1.00E+00
rs10459466	14q32.33	NA	0.1281	7.551 (2.54–22.45)	1.00E+00
rs11845269	14q32.33	NA	0.1296	6.731 (2.353–19.26)	1.00E+00
rs17113276	14q32.33	NA	0.1296	6.731 (2.353–19.26)	1.00E+00
rs10163015	15q26.2	MCTP2	0.1387	6.506 (2.595–16.31)	1.00E+00
rs4119133	4q32.1	NA	0.1618	6.003 (2.418–14.9)	1.00E+00
rs12499551	4q32.1	NA	0.1618	6.003 (2.418–14.9)	1.00E+00
rs4684126	3p25.2	IQSEC1	0.1069	5.891 (2.347–14.79)	1.00E+00
rs10133227	14q32.33	NA	0.1412	5.779 (2.273–14.69)	1.00E+00
rs8022729	14q32.33	NA	0.1428	5.769 (2.269–14.67)	1.00E+00
rs6536428	4q32.1	NA	0.143	5.763 (2.262–14.68)	1.00E+00
rs4446584	6q22.31	NA	0.1446	5.666 (2.328–13.79)	1.00E+00
rs6767561	3p25.2	IQSEC1	0.1071	5.615 (2.265–13.92)	1.00E+00
rs3742689	14q31.3	KCNK10	0.127	5.583 (2.305–13.52)	1.00E+00
rs279626	8q24.13	NA	0.1205	5.452 (2.236–13.29)	1.00E+00
rs4690943	4q32.1	NA	0.1433	5.413 (2.159–13.57)	1.00E+00
rs10108494	8q24.13	NA	0.1409	5.129 (2.22–11.85)	1.00E+00
rs10095460	8q24.13	NA	0.1412	5.123 (2.216–11.84)	1.00E+00
rs254411	5q14.1	NA	0.1445	5.066 (2.392–10.73)	1.00E+00
rs4540278	6q22.31	NA	0.1432	5.029 (2.126–11.9)	1.00E+00
rs279617	8q24.13	NA	0.1191	4.933 (2.056–11.84)	1.00E+00
rs7761112	6q22.31	NA	0.1432	4.85 (2.071–11.36)	1.00E+00

SNPs : single-nucleotide polymorphisms, OR : odds ratio, MAF : minor allele frequency, CI : confidence interval, NA : not available

DISCUSSION

In this study, we investigated the association of clinical and genetic factors for LS defined by three features found on simple radiographs. The previously identified genetic risk factors for LS using a candidate approach included genes encoding extracellular matrix proteins expressed in the nucleus pulposus (inner structure) and annulus fibrosus (outer layer) of the disc, such as type IX collagen (COL9A2 and COL9A3), AGC1, and cartilage intermediate layer protein (CILP)15). These findings suggest that LS is caused by changes in the structural integrity of the intervertebral disc. Given that OA and LS are both degenerative diseases of the skeletal joints, and because articular cartilage and intervertebral discs share similar patterns of gene expression, it is expected that subjects with OA and LS have a similar genetic susceptibility. A systematic review of genetic associations for peripheral joint OA and degenerative disease of the spine revealed the difficulties related to complex phenotypes of these diseases. While many logical and reporting problems, including missing population details, multiple testing, and an over-reliance on subgroup analysis, were found, cases in which significant associations were replicated in independent studies were also identified15). Some of the genes thus identified are of functional importance. For example, ASPN, a member of the small leucine-rich proteoglycan family, inhibits in vitro chondrogenesis and the expression of COL2A1 and AGC1 through the inhibition of TGF-b signaling19).

A recent systematic search of the literature, including 52 studies, identified ASPN (D-repeat), COL11A1 (rs1676486), growth differentiation factor 5 (GDF5) (rs143383), SKT (rs16924573), THBS2 (rs9406328), and MMP9 (rs17576) as genes associated with LS in humans defined by MRI with a moderate level of evidence5). However, the phenotype definition of lumbar disc degeneration was highly variable among the studies, including a decrease in disc signal intensity or disc height, disc bulges, disc herniations without specification of the symptoms, Modic changes, osteophytes, and lumbar spinal stenosis. In addition, the phenotype of disc degeneration varied between the initial and replication studies. As a result, the replications were inconsistent, and most of the associations were presented with a weak level of evidence.

We used simple x-ray for diagnose and grade lumbar spondylosis. The K-L grade is most popular grading system with classification into five grade scales (0–4) where K-L grade ≥2 is the conventional standard of the diagnosis8). Epidemiological studies showed that K-L grade was associated with the degree of low back pain in elderly subjects4,11).

The variation in the identified genetic associations may reflect ethnic diversity as well as phenotypic heterogeneity. While ASPN was identified as a candidate gene for lumbar disc degeneration with a moderate level of epidemiological evidence among Asians, GDF5 was instead identified among the Northern European population24). In addition, flaws inherent in a candidate gene approach, such as isolated analyses of disparate potential associations, may merely add to the growing repertoire of weak evidence. The genome-wide association study is based on the ‘common disease-common variation’ theory, which proposes that multiple common polymorphisms with a MAF of 0.5–1%, and with small effect sizes, might be predisposed to common disorders1). Because more complex modes of inheritance involving multiple genes rather than a simple monogenic Mendelian pattern appear to be in effect for LS, the discovery of its associated genes is more likely with a GWAS than by an SNP analysis. A meta-analysis of 4 GWAS that addressed LS among 4683 individuals with European ancestries was reported by Williams et al.23) LS was defined as disc space narrowing and the formation of osteophytes, as in our study. Among the four markers identified (rs17034687, rs2187689, rs7767277, and rs926849), the rs926849 marker located in the intronic region of Parkinson protein 2, the E3 ubiquitin protein ligase (PARK2) gene, on chromosome 6 remained strongest after adjusting for age and gender. The rs2187689 and rs7767277 markers were in strong linkage disequilibrium with proteasome subunit beta type 9 (large multifunctional peptidase 2) (PSMB9). Although we used the same 500K Affymetrix kit, the genes identified in our study did not overlap with those reported by Williams et al.24). Again, the discrepancy may arise from ethnic differences between the study populations and the difference in phenotype definition.

The SNP of leucine-rich repeat-containing G protein-coupled receptor 4 (LGR4) was significantly associated with both disc space narrowing and osteophytes in our study. LGR4 encodes a receptor for R-spondins (RSPOs), which play a pleiotropic role in normal development, and the development of cancer, as well as in the survival of adult stem cells through potentiation of Wnt signaling7). Knockout of LGR4 or RSPOs in mice presents with severe developmental abnormalities, causing neonatal/embryonic lethality10). Recent studies showed that the RSPO-LGR4 axis elevates the levels of Wnt receptors through direct inhibition of the ubiquitination of Wnt receptors or though interaction with intracellular signaling proteins to potentiate the Wnt pathways3). Increased Wnt/β-catenin signaling is observed in the nucleus pulposus of dogs that suffer from premature intervertebral disc degeneration18), suggesting the role of Wnt/β-catenin signaling, and possibly the RSPO-LGR4 axis, in the pathogenesis of LS.

The SNP of the heat shock transcription factor 2 gene (HSTF2) was one of the SNPs most significantly associated with LS as defined by disc space narrowing. HSTF2 specifically binds to the heat-shock promoter element and activates transcription for heat-shock response genes under conditions of heat or other stress. Previous studies showed that the formation of heterocomplexes between HSTF1 and HSTF2 leads to enhanced activity, which activates the hsp70 promoter, and that HSTF2 was able to modulate the HSTF1-mediated expression of major heat shock protein genes12). A previous study showed that HSTF expression was more frequent in clustered cells in both the annulus fibrosus and nucleus pulposus of herniated discs17). Because the K-L grade is accounted for more by osteophytes than by joint space narrowing, it is not surprising that genes associated with LS, as defined by the K-L grade, are different from those defined by disc space narrowing. The strongest association with K-L grade was observed with the SNP for transient receptor potential cation channel (TRPC), subfamily C, which encodes a receptor-activated calcium channel. In a study using OA chondrocytes, a correlation between the appearance of TRPC6 and the state of de-differentiation of chondrocytes was identified6). Because the loss of a differentiated chondrocyte phenotype is one of the hallmarks of OA, as well as LS, the role of TRPC6 in the process of disc degeneration is plausible. On the other hand, TRPC6 was expressed in bone-derived cells and inhibitors of the TRP channel inhibited the effects of bradykinin-induced Ca2+-influx, suggesting a role in bone metabolism21).

This is the first reported GWAS in an Asian population. In terms of phenotyping, we used disc space narrowing as well as osteophytosis and K-L grading to cover all aspects of degenerative change observable in LS. Limitations included a sample size too small to detect the true effect of SNPs. Although we did not identify the same loci reported in previous studies using a candidate gene approach, differences in ethnicity and phenotype definition may have affected the results. Further studies using larger Asian samples are warranted to establish if the gene variants identified in this study are associated with an increase in the risk of LS in Asian populations.

CONCLUSION

A GWAS was conducted to identify the susceptibility genes responsible for lumbar spondylosis assessed by simple radiography. We identified SNPs that potentially contribute to the pathogenesis of LS. This is the first report of a GWAS in an Asian population.

PATIENT CONSENT

The patient provided written informed consent for the publication and the use of their images.