Genome-wide association study identifies new susceptibility loci for adolescent idiopathic scoliosis in Chinese girls

Abstract

Adolescent idiopathic scoliosis (AIS) is a structural deformity of the spine affecting millions of children. As a complex disease, the genetic aetiology of AIS remains obscure. Here we report the results of a four-stage genome-wide association study (GWAS) conducted in a sample of 4,317 AIS patients and 6,016 controls. Overall, we identify three new susceptibility loci at 1p36.32 near AJAP1 (rs241215, P_combined =2.95 × 10^-9 ), 2q36.1 between PAX3 and EPHA4 (rs13398147, P_combined =7.59 × 10^-13 ) and 18q21.33 near BCL-2 (rs4940576, P_combined =2.22 × 10^-12 ). In addition, we refine a previously reported region associated with AIS at 10q24.32 (rs678741, P_combined =9.68 × 10^-37 ), which suggests LBX1AS1, encoding an antisense transcript of LBX1, might be a functional variant of AIS. This is the first GWAS investigating genetic variants associated with AIS in Chinese population, and the findings provide new insight into the multiple aetiological mechanisms of AIS.

ARTICLE

Received 28 Apr 2015 | Accepted 13 Aug 2015 | Published 22 Sep 2015

DOI: 10.1038/ncomms9355 OPEN

Genome-wide association study identies new susceptibility loci for adolescent idiopathic scoliosis in Chinese girls

Zezhang Zhu1,2,*, Nelson Leung-Sang Tang2,3,4,5,*, Leilei Xu1,2,*, Xiaodong Qin1,2, Saihu Mao1,

Yueming Song6, Limin Liu6, Fangcai Li7, Peng Liu8, Long Yi9, Jiang Chang10, Long Jiang11, Bobby Kin-Wah Ng12, Benlong Shi1, Wen Zhang1, Jun Qiao1,2, Xu Sun1,2, Xusheng Qiu1,2,Zhou Wang1, Fei Wang1, Dingding Xie1, Ling Chen1, Zhonghui Chen1, Mengran Jin1, Xiao Han1, Zongshan Hu1, Zhen Zhang1, Zhen Liu1, Feng Zhu1, Bang-ping Qian1,2, Yang Yu1,2, Bing Wang1,2, K.M. Lee5, Wayne Y.W. Lee12, T.P. Lam12, Yong Qiu1,2,** & Jack Chun-Yiu Cheng2,5,12,**

Adolescent idiopathic scoliosis (AIS) is a structural deformity of the spine affecting millions of children. As a complex disease, the genetic aetiology of AIS remains obscure. Here we report the results of a four-stage genome-wide association study (GWAS) conducted in a sample of 4,317 AIS patients and 6,016 controls. Overall, we identify three new susceptibility loci at 1p36.32 near AJAP1 (rs241215, Pcombined 2.95 10 9), 2q36.1 between PAX3 and

EPHA4 (rs13398147, Pcombined 7.59 10 13) and 18q21.33 near BCL-2 (rs4940576, Pcombined 2.22 10 12). In addition, we rene a previously reported region associated with

AIS at 10q24.32 (rs678741, Pcombined 9.68 10 37), which suggests LBX1AS1, encoding an antisense transcript of LBX1, might be a functional variant of AIS. This is the rst GWAS investigating genetic variants associated with AIS in Chinese population, and the ndings provide new insight into the multiple aetiological mechanisms of AIS.

1 Department of Spine Surgery, the Afliated Drum Tower Hospital of Nanjing University Medical School, Nanjing 210008, China. 2 Joint Scoliosis Research Center of The Chinese University of Hong Kong and Nanjing University, Nanjing 210008, China. 3 Department of Chemical Pathology, Faculty of Medicine, The Chinese University of Hong Kong, Hong Kong, China. 4 School of Biomedical Sciences, Faculty of Medicine, The Chinese University of Hong Kong, Hong Kong, China. 5 Li Ka Shing Institute of Health Sciences, Faculty of Medicine, The Chinese University of Hong Kong, Hong Kong, China. 6 Department of Orthopaedic, The West China Hospital, Sichuan University, Chengdu 610000, China. 7 Department of Orthopaedic, The Second Afliated Hospital of Zhejiang University School of Medicine, Hangzhou 310000, China. 8 Department of Orthopaedic, China-Japan Union Hospital of Jilin University, Changchun 130022, China. 9 Jiangsu Key Laboratory for Molecular Medicine, Nanjing University Medical School, Nanjing 210008, China. 10 State Key Laboratory of Environment Health (Incubation), MOE (Ministry of Education) Key Laboratory of Environment & Health, Ministry of Environmental Protection Key Laboratory of Environment and Health (Wuhan), School of Public Health, Tongji Medical College, Huazhong University of Science and Technology, Wuhan 430000, China. 11 Department of Orthopaedic, Yixing People Hospital, Wuxi 214200, China. 12 Department of Orthopaedics and Traumatology, Faculty of Medicine, The Chinese University of Hong Kong, Hong Kong, China. * These authors contributed equally to this work. ** These authors jointly supervised this work. Correspondence and requests for materials should be addressed to Y.Q. (email: mailto:[email protected]

Web End [email protected] ) or to J.C.-Y.C. (email: mailto:[email protected]

Web End [email protected] ).

NATURE COMMUNICATIONS | 6:8355 | DOI: 10.1038/ncomms9355 | http://www.nature.com/naturecommunications

Web End =www.nature.com/naturecommunications 1

& 2015 Macmillan Publishers Limited. All rights reserved.

ARTICLE NATURE COMMUNICATIONS | DOI: 10.1038/ncomms9355

Adolescent idiopathic scoliosis (AIS) is a structural deformity of the spine that is estimated to affect millions of children, with a prevalence of 24% (refs 1,2).

Commonly understood as a complex trait (polygenic disease), the genetic aetiology of AIS has been widely investigated by linkage analysis and candidate gene association analysis35. During the era of candidate study, several genes were reported to be associated with AIS, such as ERa (ref. 6), MATN1 (ref. 7),

MTNR1B (ref. 8) and TGFB1 (ref. 9). However, few of them have been successfully replicated in different ethnicities1013. The rst genome-wide association study (GWAS) in a Caucasian population reported that rs1400180 of CHL1 was associated with AIS14, which however could not be replicated in a Chinese Han population15. Subsequently, a GWAS performed in a Japanese population revealed that LBX1 and GPR126 could play a role in the etiopathogenesis of AIS16,17. In a recent study that combined the GWAS data of AIS in Caucasian and Japanese populations, Sharma et al.18 reported an enhancer locus of PAX1 was associated with the susceptibility of female AIS. As estimated by Kou et al.17, variants of LBX1 and GPR126 can only explain B1% of the trait variance in AIS. Herein additional genetic risk factors need to be discovered for a better understanding of the genetic susceptibility of AIS.

To map novel risk loci for AIS, here we report the results of a four-stage GWAS in a Chinese population. Overall, we identify three new susceptibility loci: 1p36.32 near AJAP1, encoding adherens junction associated protein 1 (rs241215, Pcombined 2.95

10 9); 2q36.1 between PAX3 and EPHA4, encoding paired box 3 and EPH receptor A4, respectively (rs13398147,

Pcombined 7.59 10 13); and 18q21.33 near BCL-2, encoding

B-cell CLL/lymphoma 2 (rs4940576, Pcombined 2.22 10 12).

In addition, we rene a previously reported region associated with AIS at 10q24.32 (rs678741, Pcombined 9.68 10 37), which

suggests LBX1AS1, encoding an antisense transcript of LBX1 (lady-bird like homeobox), might be the functional variant of AIS. This is the rst GWAS to our knowledge investigating the genetic variants associated with AIS in a Chinese population, and the ndings provide new insight into the multiple aetiological mechanisms of AIS.

ResultsAssociation analyses in the discovery and replication. We performed a four-stage GWAS in cohorts of 4,317 patients and 6,016 controls, all of the Chinese Han population. In the discovery stage, 906,703 single nucleotide polymorphisms (SNPs) were genotyped in 1,001 patients presenting a main thoracic curve (major curve located in thoracic region) and 1,500 healthy controls. After quality control ltering of the genotyping data, 516,220 autosomal SNPs in 960 cases and 1,499 controls were retained for the subsequent CochranArmitage trend test. Analysis of population stratication with principal-component analysis (PCA) showed that all subjects of our study were clustered in Chinese population (Supplementary Fig. 1). The genomic ination factor (l) as shown by the quantilequantile plot was 1.031, indicating that there was only a minimal amount of population stratication (Supplementary Fig. 2). As shown by the Manhattan plot (Fig. 1), 43 SNPs were observed to surpass the threshold of suggestive whole-genome signicance (Po1.0 10 5) (Supplementary Table 2 and 3).

From each putative haplotype block dened by linkage disequilibrium (r24 0.8), we chose one representative SNP with the lowest P value for replication. Sixteen representative SNPs (Supplementary Table 2) were then included in replication stage 1, composed of 1,916 cases and 2,001 controls. Signicant associations were found in four SNPs, which were further genotyped in 1,400 cases and 2,515 controls. Combining the results from the discovery stage and three replication stages, four SNPs were found to reach a genome-wide signicance threshold of Po5 10 8, including rs241215 at 1p36, rs13398147 at 2q36,

rs678741 at 10q24 and rs4940576 at 18p21 (Table 1). We performed imputation analysis around the four signicant regions, and the logistic regression model showed similar association with those of the discovery stage (Supplementary Table 4). The strongest evidence of association was attained with rs678741 (Pcombined 9.68 10 37, odds ratio (OR) 1.44, 95%

condential interval (CI): 1.371.52, Phet 0.95) located in the

intron of the LBX1AS1 (Table 1, Fig. 2a and Supplementary Fig. 3). SNP rs4940576 (Pcombined 2.22 10 12, OR 1.23,

95% CI: 1.141.29, Phet 0.13) was located in the intronic region

Log10 (P)

15.014.514.013.513.012.512.011.511.010.510.09.59.08.58.07.57.06.56.05.55.04.54.03.53.02.52.01.51.00.5 0

Chr1

Chr2

Chr3

Chr4

Chr5

Chr6

Chr7

Chr8

Chr9

Chr10

Chr11

Chr12

Chr13

Chr14

Chr15

Chr16

Chr17

Chr18

Chr19

Chr20

Chr21

Chr22

Figure 1 | Genome-wide association results for AIS in Han Chinese girls. Here we present Manhattan plot of P values on the log10 scale in the

discovery stage (960 cases and 1,499 controls). The x axis shows Chromosomes as delineated by different colours. The y axis shows log10 P values of

each SNP. The red line represents P 5 10 8, and the blue dashed line represents P 1 10 5.

2 NATURE COMMUNICATIONS | 6:8355 | DOI: 10.1038/ncomms9355 | http://www.nature.com/naturecommunications

Web End =www.nature.com/naturecommunications

& 2015 Macmillan Publishers Limited. All rights reserved.

NATURE COMMUNICATIONS | DOI: 10.1038/ncomms9355 ARTICLE

Table 1 | Association results for four SNPs in the GWAS, replication and combined samples.

SNP Chr. Genes MA Stage MAF P value OR (95% CI) Phet

Cases Controlrs678741 10q24.32 LBX1AS1 A GWAS 0.532 0.444 1.79 10 9 1.42 (1.271.60)

Replication 1 0.546 0.455 5.96 10 14 1.44 (1.321.58)

Replication 2 0.549 0.452 3.83 10 12 1.48 (1.321.66)

Replication 3 0.538 0.453 1.23 10 4 1.41 (1.181.68)

Combined 0.543 0.451 9.68 10 37 1.44 (1.371.52) 0.95

rs4940576 18q21.33 BCL-2 T GWAS 0.471 0.393 7.01 10 8 1.38 (1.231.54)

Replication 1 0.419 0.384 2.06 10 3 1.16 (1.061.27)

Replication 2 0.430 0.386 1.68 10 3 1.20 (1.071.35)

Replication 3 0.431 0.382 2.37 10 2 1.23 (1.031.46)

Combined 0.435 0.387 2.22 10 12 1.23 (1.141.29) 0.13

rs13398147 2q36.1 PAX3/EPHA4 T GWAS 0.267 0.198 1.47 10 8 1.48 (1.291.70)

Replication 1 0.233 0.201 9.59 10 4 1.20 (1.081.34)

Replication 2 0.239 0.202 1.39 10 3 1.24 (1.091.42)

Replication 3 0.241 0.202 3.38 10 2 1.25 (1.021.54)

Combined 0.243 0.201 7.59 10 13 1.28 (1.201.37) 0.11

rs241215 1p36.32 AJAP1 A GWAS 0.261 0.321 7.55 10 6 0.75 (0.650.85)

Replication 1 0.299 0.329 5.49 10 3 0.87 (0.790.96)

Replication 2 0.291 0.323 7.09 10 3 0.85 (0.750.97)

Replication 3 0.284 0.327 3.71 10 2 0.82 (0.680.99)

Combined 0.287 0.325 2.95 10 9 0.83 (0.780.89) 0.32

Chr., chromosome; MA, minor allele; MAF, minor allele frequency; OR, odds ratio for the minor allele; P , P value from the heterogeneity test based on GWAS and the two replication study groups.

a b

r

0.8

0.60.40.2

100

Recombination rate (cM/Mb)

r

0.80.60.40.2

100

Recombination rate (cM/Mb)

Log10 (Pvalue)

Log10 (Pvalue)

15 10

8

6

4

2

Log10 (Pvalue)

10

5

0

80

60

40

20

0

100

Recombination rate (cM/Mb)

80

60

40

20

0

100

Recombination rate (cM/Mb)

0

102.7

102.8 102.9 103 103.1 103.2 103.3 58.6 58.8 59 59.2Position on chr10 (Mb)

Position on chr2 (Mb)

c d

10

8

6

4

2

10

r

0.80.60.40.2

r

0.80.60.40.2

8

6

4

2

0

80

60

40

20

0

Log10 (Pvalue)

80

60

40

20

0

0

222.2 222.4 222.6 222.8

Position on chr18 (Mb)

4.2 4.3 4.4 4.5 4.6 4.7 4.8 Position on chr1 (Mb)

Figure 2 | Regional association plots of four susceptibility loci for AIS. Each plot shows log 10 (P value) for SNPs in the specic region. Imputation was

performed by using 1000 Genomes Project CHB (Han Chinese in Beijing, China) data (November 2010 release) as a reference. The genotyped SNP with the strongest association signal in each locus is indicated by a purple diamond, and the other SNPs are coloured according to the r2 values with the proxy SNP. The genes within the regions of interest are annotated with the direction of transcription represented by arrows. (a) 10q24.32, (b) 18q21.33, (c) 2q36.1 and (d) 1p36.32.

NATURE COMMUNICATIONS | 6:8355 | DOI: 10.1038/ncomms9355 | http://www.nature.com/naturecommunications

Web End =www.nature.com/naturecommunications 3

& 2015 Macmillan Publishers Limited. All rights reserved.

ARTICLE NATURE COMMUNICATIONS | DOI: 10.1038/ncomms9355

of BCL-2 (Table 1, Fig. 2b and Supplementary Fig. 4). SNP rs13398147 (Pcombined 7.59 10 13, OR 1.28, 95% CI:

1.201.37, Phet 0.11) was located 323 kb upstream of EPHA4

and 304 kb downstream of PAX3 (Table 1, Fig. 2c and Supplementary Fig. 5). SNP rs241215 (Pcombined 2.95 10 9,

OR 0.83, 95% CI: 0.780.89, Phet 0.32) maps 111 kb upstream

of AJAP1 (Table 1, Fig. 2d and Supplementary Fig. 6).

Functional annotation. We evaluated the overlap of the novel associated SNPs with the Encyclopedia of DNA Elements (ENCODE)-annotated genomic elements19. Some signs of regulatory activity, including histone modications at promoter or enhancer, DNase hypersensitivity sites, binding proteins and motifs changed, were observed for the four associated SNPs and their surrogates (Supplementary Tables 58).

DiscussionIn this four-stage GWAS, we identied four SNPs that are signicantly associated with risk of AIS in the Chinese population. The linkage disequilibrium block of the top signal SNP (rs678741) at 10q24.32 contains a widely validated risk variant of AIS (rs11190870) mapping to the 30 region of LBX1 (refs 16,2022). Highly expressed in spinal cord and skeletal muscle, LBX1 was identied as an AIS susceptibility gene by a previous GWAS16. LBX1 may be implicated in the aetiology of scoliosis through abnormal somatosensory function16. In the current study, we further pinpointed a potentially regulatory SNP in this region. Data from ENCODE showed that rs678741 may be located in a strong enhancer region marked by peaks of several active histone methylation modications. In addition, a search of HaploReg indicated that rs678741 maps to a DNase hypersensitivity site in multiple cell lines (Supplementary Table 5)23. These ndings indicated that rs678741 has the potential to regulate the transcriptional activity of LBX1AS1. Previous studies have shown that antisense long non-coding RNA can regulate the expression of target genes by directly binding to the transcript24,25. Therefore, it is possible that LBX1AS1 can be a biologically functional element that contributes to the risk of AIS. Further studies are necessary to clarify the role of LBX1AS1 in the pathogenesis of AIS.

SNP rs4940576 is located in the intron of BCL-2 gene that encodes an integral outer mitochondrial membrane protein playing a key role in apoptosis26,27. Amling et al.28 observed that BCL-2 is expressed at highest levels in late proliferative and maturing chondrocytes. In BCL-2 knockout mice, endochondral ossication was accelerated owing to the premature loss of terminal hypertrophic chondrocytes in the growth plate29. Thus, the product of BCL-2 appears to have an important role in controlling osteoblast activity and regulating the endochondralossication process. Interestingly, overgrowth of anterior spinal endplate is a classic hypothesis in the etiopathogenesis of AIS30. Proliferative and hypertrophic chondrocytes in the anterior column of AIS patients seemed more active than that of the posterior column31. Moreover, AIS patients were found to have different growth kinetics in bilateral growth plate of the vertebrae as indicated by differences in cellular activity between the convex and concave side of the curve32.

The third SNP rs13398147 is located between EPHA4 and PAX3. EPHA4 belongs to the EPH receptor subfamily of the proteintyrosine kinase family33. EPH receptors are implicated in mediating developmental events, particularly in the nervous system34. The PAX3 subfamily regulates both myogenesis and neurogenesis in the neural tube3537. The neural tube/notochord complex has a critical function during the development of vertebral muscle38. In addition, PAX3 mutation can lead to

muscular and neural tube defects, as well as malformation of the vertebral column39,40. Abnormality of the paravertebral muscles has been proposed as the cause of AIS30. Many diseases that affect muscle function can present a secondary scoliosis41,42.The fourth SNP rs241215 maps near AJAP1, encoding a transmembrane protein that interacts with E-cadherin/b-catenin complexes and the adaptor protein complex AP-1B in polarized epithelial cells43. Recent ndings suggested a potential role for AJAP1 in cellcell and cell-extracellular matrix interactions that could be involved in cell adhesion, migration and invasion44,45. Regulation of cell adhesion has been found to be essential in the bone growth and osteoblast differentiation4648.

To determine if there was any association between the novel associated SNPs and curve severity, we investigated a subgroup of 632 patients with a mean Cobb angle of 37.29.4 (range 2749). Of these, 214 patients received fusion surgery, and all other 418 who did not require surgical correction had been observed until skeletal maturity. As shown in Supplementary Table 9, rs241215 might be associated with the curve severity. Patients with genotype AA had a higher curve magnitude (39.67.3) than those with genotype TA(36.97.7) or TT (36.78.1). As for the other three SNPs, no difference regarding the curve severity was found among the genotypes.

The primary limitation of our study lies in the lack of experiments supporting the functional role of the newly identied susceptibility genes. We believe that future functional analysis concerning the expression of these genes in tissues of AIS patients and related regulatory pathways can shed light on their relationships with the development of AIS. Moreover, future ne-mapping studies based on data from the resequencing of these regions may provide a clear understanding of the association at these loci.

In summary, our GWAS conducted in a Chinese Han population identies three new susceptibility loci of AIS at 1p36.32, 2q36.1 and 18q21.33. This study also renes a region previously associated with AIS at 10q24.32, suggesting that LBX1AS1, an antisense transcript of LBX1, might be the functional variant of AIS. These ndings provide additional insights into the genetic architecture of AIS. Further dissection of these loci is likely to expand our understanding of the aetiology of AIS.

Methods

Subjects. The current case-control analysis was composed of four stages, including the initial genome-wide discovery stage followed by three replication stages. The study has been approved by the following local Institutional Review Boards (IRB): IRB of Nanjing Drum Tower Hospital, IRB of The Chinese University of Hong Kong, IRB of The West China Hospital, IRB of The Second Afliated Hospital of Zhejiang University and IRB of China-Japan Union Hospital of Jilin University. Patients who came to our joint scoliosis research clinics between 2000 and 2015 were reviewed for the eligibility to be included in this study. The inclusion criteria were as follows: (1) female; (2) diagnosed as AIS through clinical and radiologic examinations; (3) with Cobb angle 420; and (4) having right major thoracic curve. Patients with scoliosis secondary to known aetiology, such as congenital scoliosis, neuromuscular scoliosis, scoliosis secondary to skeletal dysplasia and connective tissue abnormalities, were excluded from the study. Healthy subjects who came to our hospital for routine physical examinations were recruited as controls. Overall, there were 4,317 patients and 6,016 controls included in our study, all from ethnic Han Chinese population and collected in Nanjing of Jiangsu Province, Hangzhou of Zhejiang Province, Chengdu of Sichuan province, Changchun of Jilin Province and Hong Kong. The detailed prole of the subjects was as follows: in the discovery stage, 728 cases and 1,500 controls were collected in Nanjing, and the other 273 cases were collected in Hong Kong; in replication stage 1, all the 1,916 cases and 2,001 controls were collected in Nanjing; in replication stage 2, 307 cases and 1,735 controls were collected in Nanjing, and 600 cases and 180 controls were collected in Hong Kong, 15 cases were collected in Hangzhou, 10 cases were collected in Chengdu and 18 cases were collected in Changchun; in replications stage 3, all the 450 cases and 600 controls were collected in Nanjing. Under the approval of each participating centres Ethical Committee, informed

4 NATURE COMMUNICATIONS | 6:8355 | DOI: 10.1038/ncomms9355 | http://www.nature.com/naturecommunications

Web End =www.nature.com/naturecommunications

& 2015 Macmillan Publishers Limited. All rights reserved.

NATURE COMMUNICATIONS | DOI: 10.1038/ncomms9355 ARTICLE

consent was obtained from all the participants and from the guardians of the adolescents. Clinical characteristics of the participants in the discovery stage were summarized in Supplementary Table 1.

Genotyping and quality control for GWAS. Genomic DNA were extracted from peripheral blood leukocytes using a commercial DNA extraction kit (QIAGEN). Genotype analysis in the discovery stage was conducted using the Affymetrix Genome-Wide Human SNP Array 6.0 according to the manufacturers protocol. Systemic quality control on the raw genotyping data in the discovery stage was carried out using PLINK (v1.90). The exclusion of unqualied SNPs was performed on the basis of following criteria: (1) with call rates o95%; (2) with minor allele frequency o0.05; (3) showing deviation from HardyWeinberg equilibrium (Po1 10 5); and (4) mapping on non-autosomal Chromosomes. For sample

quality control, cryptic relatedness for each sample was evaluated with an identity-by-state method, and subjects showing second-degree relatedness or closer were removed (n 35). In addition, seven individuals with more than 10% missing

genotypes were removed. There was no individual showing gender discrepancies based on their X chromosome genotypes. To detect population outliers and stratication, we used PCA implemented in the software package EIGENSTRAT49. A total of 2,459 participants remaining after sample quality control were analysed together with 210 HapMap subjects (60 YRI, 60 CEU, 45 JPT and 45 CHB individuals). Using the top two associated principal components, we identied that all individuals were of Chinese ancestry. Subsequently, we performed PCA using only the genotype information of the case and control subjects, and the scatterplot indicated that all cases and controls were genetically matched with minimal evidence of population stratication (Supplementary Fig. 1). After the quality control process, a total of 960 cases and 1,499 controls within Chinese Han population and 516,220 SNPs were used in the nal analysis.

Selection of SNPs for the replication stage. SNPs surpassing a threshold of suggestive genome-wide signicance (Po1 10 5) in the discovery stage were

potential candidates for replication. Among these SNPs, those with ambiguous genotype scatter plots were excluded from the replication stage. A total of 43 SNPs remained after ltering, from which we chose 16 representative SNPs for replication stage 1 and removed the other 27 SNPs due to high linkage disequilibrium with at least 1 of the 16 representative SNPs (Supplementary Table 2 and 3). After replication stage 1, 4 SNPs were included in replication stage 2 and 3, including rs678741, rs241215, rs13398147 and rs4940576 (Table 1 and Fig. 1). The SNP genotyping in replication stage 1 was performed with Sequenom MassARRAY system (Sequenom Inc.). Sixteen SNPs were multiplexed in one assay, and the primers were designed by custom software to give unique mass ranges for each SNP. Samples of replication stage 2 and 3 were genotyped using TaqMan SNP Genotyping Assay, which was read with an ABI Step-One-Plus sequence detection system (Applied Biosystems, Foster City, CA). Approximately 20 ng of genomic DNA was used to genotype each sample. Several quality control measures were applied as listed below. Laboratory technicians performing the genotyping assay were kept blind of case or control samples. In all, 5% of the samples were randomly selected as blind duplicates, with a reproducibility rate of 100%. DNA samples with 410% missing genotyping were excluded.

Statistical analysis. In the discovery stage, a CochranArmitage trend test was used to calculate the association of each SNP with the disease. We calculated OR values and 95% CIs from a 2 2 allele frequency table. Data from the discovery

stage and three replication stages were combined using the MantelHaenszel method. The heterogeneity between studies was examined using the Q-statisticP value50. The ination factor l was calculated by dividing the mean of the lower 90% of the test statistics by the mean of the lower 90% of the expected values from a w2 distribution with 1 degree of freedom. We used PLINK 1.90 for general statistical analysis51. The Manhattan plot of log10 P was generated using

Haploview (v4.2)52. A quantilequantile plot generated with R (v2.6) was used to evaluate the overall signicance of the GWAS results and the potential impact of population stratication.

Imputation. We performed genotype imputation within 800 kb around the four signicant regions using MaCH-Admix software53. We used the linkage disequilibrium and haplotype information from the 1000 Genomes Project (phased CHB and CHS data; March 2012) as the reference set to impute ungenotyped markers (Supplementary Table 4). SNPs with minor allele frequency r0.05 or with a low imputation quality score (R2r0.30) were excluded to ensure good imputation quality. Logistic regression analysis was used to determine the association of the SNPs with the AIS risk. The linkage disequilibrium pattern around the risk-associated SNPs was analysed to identify susceptibility genes that were covered by the linkage disequilibrium blocks. Linkage disequilibrium structures and haplotype block plots were generated with Haploview (v4.2)52. The four signicant regions were plotted using the online tool LocusZoom54.

Functional annotation. To explore the regulatory properties of the association signals, we used chromatin state segmentation in LCL data generated by the

ENCODE project19. HaploReg was used to examine whether any of the SNPs or their proxies (that is, r24 0.8 in the 1000 Genomes CHB reference panel) annotate putative transcription factor binding or enhancer elements (Supplementary Table 5-8)23.

Subgroup analysis. To investigate the relationship between the novel associated SNPs and curve severity, two types of patients were included in the subgroup analysis: (1) being consistently observed until skeletal maturity; and (2) undergoing fusion surgery due to curve progression. Skeletal maturity was determined based on chronological age (Z16-years-old) or Risser sign (44). The curve severity was recorded at the latest visit for patients being consistently observed, or at the last preoperative visit for those undergoing surgery. One-way analysis of variance test was used to compare the curve severity among different genotypes. Statistical signicance was set at a P value of o0.05.

References

1. Miller, N. H. Cause and natural history of adolescent idiopathic scoliosis. Orthop. Clin. North Am. 30, 343352 (1999).

2. Weinstein, S. L. Adolescent idiopathic scoliosis: prevalence and natural history. Instr. Course Lect. 38, 115128 (1989).

3. Wise, C. A., Gao, X., Shoemaker, S., Gordon, D. & Herring, J. A. Understanding genetic factors in idiopathic scoliosis, a complex disease of childhood. Curr. Genomics 9, 5159 (2008).

4. Riseborough, E. J. & Wynne-Davies, R. A genetic survey of idiopathic scoliosis in Boston, Massachusetts. J. Bone Joint Surg. Am. 55, 974982 (1973).

5. Cheng, J. C., Tang, N. L., Yeung, H. Y. & Miller, N. Genetic association of complex traits: using idiopathic scoliosis as an example. Clin. Orthop. Relat. Res. 462, 3844 (2007).

6. Inoue, M. et al. Association between estrogen receptor gene polymorphisms and curve severity of idiopathic scoliosis. Spine (Phila Pa 1976) 27, 23572362 (2002).

7. Chen, Z. et al. Promoter polymorphism of matrilin-1 gene predisposes to adolescent idiopathic scoliosis in a Chinese population. Eur. J. Hum. Genet. 17, 525532 (2009).

8. Qiu, X. S. et al. The role of melatonin receptor 1B gene (MTNR1B) in adolescent idiopathic scoliosis--a genetic association study. Stud. Health Technol. Inform. 123, 38 (2006).

9. Ryzhkov, I. I. et al. Transforming growth factor beta 1 is a novel susceptibility gene for adolescent idiopathic scoliosis. Spine (Phila Pa 1976) 38, E699E704 (2013).

10. Bae, J. W., Cho, C. H., Min, W. K. & Kim, U. K. Associations between matrilin-1 gene polymorphisms and adolescent idiopathic scoliosis curve patterns in a Korean population. Mol. Biol. Rep. 39, 55615567 (2012).

11. Takahashi, Y. et al. Lack of association between adolescent idiopathic scoliosis and previously reported single nucleotide polymorphisms in MATN1, MTNR1B, TPH1, and IGF1 in a Japanese population. J. Orthop. Res. 29, 10551058 (2011).

12. Tang, N. L. et al. A relook into the association of the estrogen receptor [alpha] gene (PvuII, XbaI) and adolescent idiopathic scoliosis: a study of 540 Chinese cases. Spine (Phila Pa 1976) 31, 24632468 (2006).

13. Takahashi, Y. et al. Replication study of the association between adolescent idiopathic scoliosis and two estrogen receptor genes. J. Orthop. Res. 29, 834837 (2011).

14. Sharma, S. et al. Genome-wide association studies of adolescent idiopathic scoliosis suggest candidate susceptibility genes. Hum. Mol. Genet. 20, 14561466 (2011).

15. Qiu, X. S. et al. Lack of association between the CHL1 gene and adolescent idiopathic scoliosis susceptibility in Han Chinese: a case-control study. BMC Musculoskelet Disord. 15, 38 (2014).

16. Takahashi, Y. et al. A genome-wide association study identies common variants near LBX1 associated with adolescent idiopathic scoliosis. Nat. Genet. 43, 12371240 (2011).

17. Kou, I. et al. Genetic variants in GPR126 are associated with adolescent idiopathic scoliosis. Nat. Genet. 45, 676679 (2013).

18. Sharma, S. et al. A PAX1 enhancer locus is associated with susceptibility to idiopathic scoliosis in females. Nat. Commun. 6, 6452 (2015).

19. Birney, E. et al. Identication and analysis of functional elements in 1% of the human genome by the ENCODE pilot project. Nature 447, 799816 (2007).

20. Jiang, H. et al. Association of rs11190870 near LBX1 with adolescent idiopathic scoliosis susceptibility in a Han Chinese population. Eur. Spine J. 22, 282286 (2013).

21. Gao, W. et al. Association between common variants near LBX1 and adolescent idiopathic scoliosis replicated in the Chinese Han population. PLoS One 8, e53234 (2013).

22. Fan, Y. H. et al. SNP rs11190870 near LBX1 is associated with adolescent idiopathic scoliosis in southern Chinese. J. Hum. Genet. 57, 244246 (2012).

NATURE COMMUNICATIONS | 6:8355 | DOI: 10.1038/ncomms9355 | http://www.nature.com/naturecommunications

Web End =www.nature.com/naturecommunications 5

& 2015 Macmillan Publishers Limited. All rights reserved.

ARTICLE NATURE COMMUNICATIONS | DOI: 10.1038/ncomms9355

23. Ward, L. D. & Kellis, M. HaploReg: a resource for exploring chromatin states, conservation, and regulatory motif alterations within sets of genetically linked variants. Nucleic Acids Res. 40, D930D934 (2012).

24. Kerin, T. et al. A noncoding RNA antisense to moesin at 5p14.1 in autism. Sci. Transl. Med. 4, 128ra40 (2012).

25. Mizuta, I. et al. YY1 binds to alpha-synuclein 3-anking region SNP and stimulates antisense noncoding RNA expression. J. Hum. Genet. 58, 711719 (2013).

26. EI-Emshaty, H. M., Saad, E. A., Toson, E. A., Abdel, M. C. & Gadelhak, N. A. Apoptosis and cell proliferation: correlation with BCL-2 and P53 oncoprotein expression in human hepatocellular carcinoma. Hepatogastroenterology 61, 13931401 (2014).

27. Siddiqui, W. A., Ahad, A. & Ahsan, H. The mystery of BCL2 family: Bcl-2 proteins and apoptosis: an update. Arch. Toxicol. 89, 289317 (2015).

28. Amling, M., Posl, M., Hentz, M. W., Priemel, M. & Delling, G. PTHrP and Bcl-2: essential regulatory molecules in chondrocyte differentiation and chondrogenic tumors. Verh Dtsch Ges Pathol. 82, 160169 (1998).

29. Boot-Handford, R. P. et al. The bcl-2 knockout mouse exhibits marked changes in osteoblast phenotype and collagen deposition in bone as well as a mild growth plate phenotype. Int. J. Exp. Pathol. 79, 329335 (1998).

30. Kouwenhoven, J. W. & Castelein, R. M. The pathogenesis of adolescent idiopathic scoliosis: review of the literature. Spine (Phila Pa 1976) 33, 28982908 (2008).

31. Zhu, F., Qiu, Y., Yeung, H. Y., Lee, K. M. & Cheng, J. C. Histomorphometric study of the spinal growth plates in idiopathic scoliosis and congenital scoliosis. Pediatr. Int. 48, 591598 (2006).

32. Wang, S. et al. Histomorphological study of the spinal growth plates from the convex side and the concave side in adolescent idiopathic scoliosis. J. Orthop. Surg. Res. 2, 19 (2007).

33. Farenc, C., Celie, P. H., Tensen, C. P., de Esch, I. J. & Siegal, G. Crystal structure of the EphA4 protein tyrosine kinase domain in the apo- and dasatinib-bound state. FEBS Lett. 585, 35933599 (2011).

34. Frisen, J. & Barbacid, M. Genetic analysis of the role of Eph receptors in the development of the mammalian nervous system. Cell Tissue Res. 290, 209215 (1997).

35. Schubert, F. R. et al. Early mesodermal phenotypes in splotch suggest a role for Pax3 in the formation of epithelial somites. Dev. Dyn. 222, 506521 (2001).

36. Buckingham, M. & Relaix, F. The role of Pax genes in the development of tissues and organs: Pax3 and Pax7 regulate muscle progenitor cell functions. Annu. Rev. Cell Dev. Biol. 23, 645673 (2007).

37. Young, A. P. & Wagers, A. J. Pax3 induces differentiation of juvenile skeletal muscle stem cells without transcriptional upregulation of canonical myogenic regulatory factors. J. Cell Sci. 123, 26322639 (2010).

38. Rong, P. M., Teillet, M. A., Ziller, C. & Le Douarin, N. M. The neural tube/ notochord complex is necessary for vertebral but not limb and body wall striated muscle differentiation. Development 115, 657672 (1992).

39. Xiao, Y. et al. Characterization of a novel missense mutation on murine Pax3 through ENU mutagenesis. J. Genet. Genomics 38, 333339 (2011).

40. Farin, H. F., Mansouri, A., Petry, M. & Kispert, A. T-box protein Tbx18 interacts with the paired box protein Pax3 in the development of the paraxial mesoderm. J. Biol. Chem. 283, 2537225380 (2008).

41. Hsu, J. D. The natural history of spine curvature progression in the nonambulatory Duchenne muscular dystrophy patient. Spine (Phila Pa 1976) 8, 771775 (1983).

42. Allam, A. M. & Schwabe, A. L. Neuromuscular scoliosis. PM&R 5, 957963 (2013).

43. Bharti, S. et al. Novel membrane protein shrew-1 targets to cadherin-mediated junctions in polarized epithelial cells. Mol. Biol. Cell 15, 397406 (2004).

44. McDonald, J. M. et al. The SHREW1 gene, frequently deleted in oligodendrogliomas, functions to inhibit cell adhesion and migration. Cancer Biol. Ther. 5, 300304 (2006).

45. Lin, N. et al. Deletion or epigenetic silencing of AJAP1 on 1p36 in glioblastoma. Mol. Cancer Res. 10, 208217 (2012).

46. Sosa-Garcia, B. et al. A role for the retinoblastoma protein as a regulator of mouse osteoblast cell adhesion: implications for osteogenesis and osteosarcoma formation. PLoS ONE 5, e13954 (2010).

47. Hay, E. et al. N- and E-cadherin mediate early human calvaria osteoblast differentiation promoted by bone morphogenetic protein-2. J. Cell Physiol. 183, 117128 (2000).

48. Li, H. et al. Role of neural-cadherin in early osteoblastic differentiation of human bone marrow stromal cells cocultured with human umbilical vein endothelial cells. Am. J. Physiol. Cell Physiol. 299, C422C430 (2010).

49. Price, A. L. et al. Principal components analysis corrects for stratication in genome-wide association studies. Nat. Genet. 38, 904909 (2006).

50. Mantel, N. & Haenszel, W. Statistical aspects of the analysis of data from retrospective studies of disease. J. Natl Cancer Inst. 22, 719748 (1959).

51. Purcell, S. et al. PLINK: a tool set for whole-genome association and population-based linkage analyses. Am. J. Hum. Genet. 81, 559575 (2007).

52. Barrett, J. C., Fry, B., Maller, J. & Daly, M. J. Haploview: analysis and visualization of LD and haplotype maps. Bioinformatics 21, 263265 (2005).

53. Liu, E. Y., Li, M., Wang, W. & Li, Y. MaCH-admix: genotype imputation for admixed populations. Genet. Epidemiol. 37, 2537 (2013).

54. Pruim, R. J. et al. LocusZoom: regional visualization of genome-wide association scan results. Bioinformatics 26, 23362337 (2010).

Acknowledgements

We thank all the research staff and students who participated in this research and devoted their time and effort. We sincerely thank Dr Gorman Kristen of CHU Saint Justine Research Center (Montreal, Canada) for advice on language and style of the manuscript. We especially thank Linlin Zhang, Le Zhang, Jian Wang and Lili Sha for their great help in solving the administration issues. Besides, we thank Wei Liu for his coordination of the sample transportation. This work was funded by the Science and Technology Project of Jiangsu Provinces Clinical Medicine (BL2012002), Jiangsu Provinces Key Medical Center (ZX201107) and the Natural Science Foundation of China (81171767, and 81171672). J.C.-Y.C. and N.L.-S.T. received grant support from Fondation Yves Cotrel, Institut de France.

Author contributions

Z.Z., N.L.-S.T., Y.Q., and J.C.-Y.C. directed the study and obtained nancial support.L.X., and J.C. performed the statistical analysis. L.X. and N.L.-S.T. processed the data and drafted the manuscript. B.-p.Q., S.M., J.Q., B.S., F.W., X.S., Z.L. and F.Z. were responsible for the collection of blood samples in Nanjing. Z.W., D.X., X.Q., Z.C., L.C., Z.H., M.J.,X.H., Z.Z., and L.J. were responsible for the extraction of DNA of Nanjing samples. N.L.-S.T., K.M.L., B.K.-W.N., W.Y.W.L. and T.P.L. were responsible for subject recruitment and DNA preparation for the Hong Kong samples. Y.S., L.L., P.L., and F.L. were responsible for subject recruitment and DNA preparation in Chengdu, Changchun and HangZhou, respectively. All authors approved the nal version of the manuscript.

Additional information

Supplementary Information accompanies this paper at http://www.nature.com/naturecommunications

Web End =http://www.nature.com/ http://www.nature.com/naturecommunications

Web End =naturecommunications

Competing nancial interests: The authors declare no competing nancial interests.

Reprints and permission information is available online at http://npg.nature.com/reprintsandpermissions/

Web End =http://npg.nature.com/ http://npg.nature.com/reprintsandpermissions/

Web End =reprintsandpermissions/

How to cite this article: Zhu, Z. et al. Genome-wide association study identies new susceptibility loci for adolescent idiopathic scoliosis in Chinese girls. Nat. Commun. 6:8355 doi: 10.1038/ncomms9355 (2015).

This work is licensed under a Creative Commons Attribution 4.0 International License. The images or other third party material in this article are included in the articles Creative Commons license, unless indicated otherwise in the credit line; if the material is not included under the Creative Commons license, users will need to obtain permission from the license holder to reproduce the material. To view a copy of this license, visit http://creativecommons.org/licenses/by/4.0/

Web End =http://creativecommons.org/licenses/by/4.0/

6 NATURE COMMUNICATIONS | 6:8355 | DOI: 10.1038/ncomms9355 | http://www.nature.com/naturecommunications

Web End =www.nature.com/naturecommunications

& 2015 Macmillan Publishers Limited. All rights reserved.