A multi-ethnic meta-analysis confirms the association of rs6570507 with adolescent idiopathic scoliosis

Adolescent idiopathic scoliosis (AIS) is the most common type of spinal deformity and has a significant genetic background. Genome-wide association studies (GWASs) identified several susceptibility loci associated with AIS. Among them is a locus on chromosome 6q24.1 that we identified by a GWAS in a Japanese cohort. The locus is represented by rs6570507 located within GPR126. To ensure the association of rs6570507 with AIS, we conducted a meta-analysis using eight cohorts from East Asia, Northern Europe and USA. The analysis included a total of 6,873 cases and 38,916 controls and yielded significant association (combined P = 2.95 × 10−20; odds ratio = 1.22), providing convincing evidence of the worldwide association between rs6570507 and AIS susceptibility. In silico analyses strongly suggested that GPR126 is a susceptibility gene at this locus.

Adolescent idiopathic scoliosis (AIS) is defined as a lateral spinal curvature that occurs without obvious cause between age 10 and skeletal maturity. The prevalence of AIS in the adolescent population is approximately 2-3% 1 and AIS occurs predominantly in females 2, 3 . AIS has been regarded as a multifactorial disease and a number of population, family and twin studies strongly suggest the importance of genetic factors in its etiology and pathogenesis [4][5][6][7] .
Genome-wide association study (GWAS) is a powerful tool to detect genetic variants associated with common diseases. Recently, several GWASs have identified several loci associated with AIS such as chromosome 10q24.31, 6q24.1, 20p11.22, and 1p36.32 loci [8][9][10][11] . Our initial GWAS identified rs11190870, a common variant on chromosome 10q24.31 that showed significant association with AIS in a Japanese cohort 8 and subsequent multi-ethnic meta-analysis provided convincing evidence of the association 12 . In addition, recent validation studies in different populations have also reported its significant associations [13][14][15] . Thus, the association of the 10q24.31 locus with AIS has been confirmed in multiple studies; however, other AIS loci have not been fully investigated.
In our previous GWAS, we identified an additional genetic variant, rs6570507 on chromosome 6q24.1, associated with AIS in a Japanese cohort 16 . The association was replicated in independent cohorts of Chinese and Caucasian in the USA 16 ; however, additional studies would be necessary to confirm the association and to identify the susceptibility gene at the locus.
In this study, we performed a meta-analysis using multi-ethnic cohorts of ~46,000 subjects. The result provided convincing evidence of the association of rs6570507 with AIS susceptibility, suggesting that the chromosome 6q24.1 locus is related to the global risk of AIS. In silico analyses strongly suggested that GPR126 is a susceptibility gene at this locus.

Results
Association of rs6570507 and AIS susceptibility. We conducted the meta-analysis of rs6570507 using eight cohorts (Table 1, Fig. 1). Three cohorts were previously reported 16 , and the other five were recruited for this study that included cohorts from Guangzhou (case 647 and control 1,048), Hong Kong (case 300 and control 788), Beijing (case 482 and control 861), USA (case 1,360 and control 7,267), and Scandinavia (case 1,522 and control 1,804). Finally, 6,873 cases and 38,916 controls were included in the meta-analysis. The analysis showed a convincing association between rs6570507 and AIS: combined P = 2.95 × 10 −20 ; odds ratio (OR) = 1.22; 95% confidence interval (CI) = 1.17-1.27 (Table 1, Supplementary Fig. S1). ORs were >1 in all eight cohorts, with little difference between ethnic groups. The analysis did not show any significant heterogeneity, suggesting no statistical difference between studies. Because AIS has clinical evidence of sexual dimorphism 17 , we performed gender-stratified analyses to determine whether a genetic difference existed between males and females ( Table 2, Supplementary Fig. S1). However, no gender difference was observed in the association.
Fine mapping. The sentinel SNP rs6570507 is present in an intron of GPR126 (encoding G protein-coupled receptor 126), and GPR126 is the only gene contained within the linkage disequilibrium (LD) block (r 2 > 0.8) represented by rs6570507 (ref. 16 ). To identify the candidate susceptibility gene in the locus, we evaluated the topologically associated domains (TADs) around the associated SNPs. Hi-C data 18 (http://promoter.bx.psu.edu/hi-c/view. php) revealed that GPR126 and VTA1 (encoding vesicle trafficking 1) were included in the TAD that contained the LD block (Fig. 2). Using expression quantitative trait loci (eQTL) data from the Genotype-Tissue Expression (GTEx) project 19 , we investigated the target genes regulated by rs6570507. We found that the expression level of GPR126 in subcutaneous adipose tissue and sun-exposed skin was significantly associated with rs6570507 and its risk allele rs6570507-A decreased the expression (Supplementary Fig. S2). These data strongly suggested that GPR126 is the most plausible AIS susceptibility gene at this locus.  [20][21][22] . To assess the potential contribution of rare variants at the chromosome 6q24.1 locus to AIS, we carried out targeted sequencing. We screened the coding regions and exon-intron boundaries of VTA1 (chr.6: 142468421-142539785) and GPR126 (chr.6: 142623463-142764657) across 10,121 individuals (2,721 cases and 7,400 controls). The coverage rate of the target region, presented as an average in all individuals covered with ≥20 reads at each position, was 99.96%. We identified 183 single nucleotide variant (SNVs), 2 insertions and 9 deletions. The single-variant association analysis for each nonsynonymous variant (including missense, splice site, nonsense, and frameshift variants) with MAF <0.05 demonstrated no rare variant with significant associations (Supplementary Table S1). Subsequently, we performed the gene-based association analyses by the cohort allelic sums test (CAST) 23 and the sequence kernel association test (SKAT) 24 , using the nonsynonymous variants.
Only GPR126 showed a nominal significant association (P = 4.65 × 10 −2 ) with AIS by CAST; however, it did not reach the threshold for significance after Bonferroni correction (Supplementary Table S2).

Discussion
In our previous study, we identified the significant association between rs6570507 and AIS in Japanese and the association is replicated in the independent cohorts of Chinese and Caucasian in the USA. However, as the numbers of Chinese and Caucasian subjects were both limited, the association between rs6570507 and AIS in these populations does not provide convincing evidence. In the current study, we have conducted a meta-analysis using eight independent multi-ethnic-cohorts and showed significant and consistent association of rs6570507 with AIS across the different ethnic groups. Recently, evidence for association of rs6570507 with AIS was reported in a Chinese cohort of Zhejiang 25 ; however, we did not use the data for our meta-analysis because the allele frequency of rs6570507 in the control group of the study was quite different from the frequency for the Chinese in the public database (for example, the 1000 Genomes Project) and the minor allele was opposite. Nevertheless, our meta-analysis demonstrates convincing association between rs6570507 and AIS susceptibility. For some common complex diseases, rare large-effect variants have recently been identified by targeted exon sequencing of genes within GWAS associated loci 26,27 . This approach is helpful to identify causal genes and to explain a further portion of disease variance. In the present study, we assessed the contribution of rare variants in the candidate genes at the locus; however, we could not find any statistically significant association in individual variants. GPR126 showed a nominally significant association by a gene-based study. The frequency of rare and low-frequency nonsynonymous variants was increased in cases, suggesting that AIS is caused by functional impairment of the GPR126 gene. However, in the current study, the evidence for association did not surpass Bonferroni correction. This negative result may be due to the insufficient sample size in the present study and further larger studies are required to evaluate the effect of rare and low-frequency variants in GPR126.
To identify the susceptibility gene at the chromosome 6q24.1 locus, we evaluated the TADs around the associated SNPs using the Hi-C database. The data showed that GPR126 and VTA1 were involved in the TAD that contained the LD block of the AIS associated SNPs. All SNPs strongly correlated with rs6570507 (r 2 > 0.8) were present in the GPR126 gene region. The eQTL data indicated that these SNPs were significantly associated with the expression level of GPR126 in some tissues, such as subcutaneous adipose tissue and sun-exposed skin. It is well known that AIS patients are associated with low body weight and low body mass index (BMI) 28,29 , and are also reported to be significantly associated with lower body fat 30 . These observations suggest that GPR126 has a potential to affect AIS susceptibility through the adipose tissues. However, tissues used in GTEx project are limited, and the sample size is different by tissue. Therefore, it may fail to identify tissues contributing to the disease Figure 1. The flow diagram of the meta-analysis using eight cohorts. Three cohorts (Japanese GWAS, replication study and Nanjing study) were previously reported 16 , and the other five were recruited for this analysis.
pathology. In order to more optimally determine the effect of the associated SNPs on scoliosis, eQTL analysis in AIS-related tissues such as intervertebral disc, cartilage and bone should be performed. Moreover, our previous experiments indicated that GPR126 is highly expressed in the cartilage of human and proliferating chondrocytes of the vertebral body in mouse embryo 16

Materials and Methods
Subjects and genotyping. Informed consent was obtained from all subjects participating in this study.
This study was approved by the ethics committee of the RIKEN center for Integrative Medical Sciences and all experiments were performed in accordance with relevant guidelines and regulations. All AIS subjects met clinical criteria including scoliosis with Cobb angle of 10° or more on standing spinal postero-anterior (P-A) radiographs and excluding other non-idiopathic forms of scoliosis. The subjects in the Japanese and Nanjing-Chinese cohorts were recruited and genotyped as previously described 16 . The details of additional studies: i.e., Guangzhou, Hong Kong, Beijing, USA, and Scandinavia studies were described as below.
Guangzhou study. AIS subjects were recruited from the First Affiliated Hospital and Sun Yat-sen Memorial Hospital of Sun Yat-sen University. They provided detailed histories, accepted physical examinations, underwent standard up-standing P-A radiography of the whole spine, and other testing such as magnetic resonance imaging (MRI), computed tomography (CT) and nuclear scintigraphy, when necessary. All AIS subjects were diagnosed at ages 10-16 years. Control subjects were recruited from individuals who received scoliosis screening at middle and primary schools in Guangzhou and fracture patients selected from the participating hospitals. All subjects were confirmed for not having AIS by x-ray scans of the spine. Routine history and physical examinations were also conducted to exclude other relevant diseases. Genomic DNA was extracted from blood using DNA Blood Mini-kit (Tiangen Biotech, Beijing, China). Primer extension sequencing (SNaPshot) assay (Applied Biosystems, CA, USA) was used for genotyping and the results were analyzed by GeneMarker software (SoftGenetics LLC, PA, USA) at Beijing Genomics Institute (Shenzhen, China) and checked by visual inspection of I.K. and H.D.
Hong Kong study. AIS subjects were recruited from the Duchess of Kent Children's Hospital in Hong Kong.
The inclusion criteria were as previously described 32 . Control subjects were randomly selected from the subjects recruited for the Genetic Study of Degenerative Disc Disease project 33 . All were confirmed not to have AIS by MRI examination of the spine. Genomic DNA was extracted from peripheral blood lymphocytes using standard procedures. The PCR-based invader assay (Third Wave Technologies, WI, USA) was used for genotyping. ). The dbGaP controls were previously genotyped on the same microarray platform used for cases. Only subjects of self-reported Non-Hispanic White were included in the present study. Phenotypes of all controls were reviewed to exclude any with musculoskeletal or neurological disorders. We applied initial per sample quality control measures and excluded sex inconsistencies and any with missing genotype rate per person more than 0.03. Remaining samples were merged using the default mode in PLINK.1.9 (ref. 34 ). Duplicated or related individuals were removed as previously described 35 . We used principal component analysis (PCA) 36 on the merged data projected onto HapMap3 samples to correct possible stratification 37 . After quality controls, 8,647 subjects (1,360 AIS subjects and 7,287 controls) were included for the current study. We applied initial per-SNPs quality control measures using PLINK including genotyping call-rate per marker (>95%), minor allele frequency (>0.01) and deviation from Hardy-Weinberg equilibrium (cutoff p-value = 10 −4 ). We noted that the SNP success rate was 99.99% for the rs6570507 and there was no significant difference in the missingness rate between cases and controls after quality controls.
Scandinavia study. AIS subjects were recruited from six hospitals in Sweden and one in Denmark as described previously to the Scoliosis and Genetics in Scandinavia (ScoliGeneS) study 15,[38][39][40] . Individuals with a history or clinical sign of a non-idiopathic scoliosis and with neural abnormalities in a MRI of the spine were excluded. All control subjects were females and recruited from the Osteoporosis Prospective Risk Assessment cohort and PEAK-25 cohort 41,42 . Dual-energy X-ray absorptiometry (DXA) scan was performed in both cohorts and subjects showing any sign of a curved spine on DXA were excluded. Genomic DNA was extracted from blood or saliva using the QIAamp 96 DNA Blood Kit and the Autopure LS system (Qiagen). iPLEX Gold chemistry and the MassARRAY system (Agena Bioscience, CA, USA) were used for genotyping. Genotype calls were checked by two persons individually using the MassARRAY Typer v4.0 Software (Agena Bioscience).
Statistical analysis. The association between rs6570507 and AIS in each study was evaluated by the Cochrane-Armitage trend test. Data from the eight studies were combined using the inverse-variance method assuming a fixed-effects model in the METAL software package (http://csg.sph.umich.edu//abecasis/Metal/) 43 .
The heterogeneity among studies was tested using Cochran's Q test based upon inverse variance weights using METAL. and the PharmaSNP Consortium 44,45 . Multiplex PCR-based target sequencing was carried out as previously described 45 . Primers were designed using the Primer 3 software (ver. 2.3.4) to obtain 180-300 bp PCR products. We amplified 9,847 bp consisting of 47 genomic fractions using three multiplex-PCR products with dual barcodes for each individual. After purifying of each library using Agencourt AMPure XP (Beckman Coulter, CA, USA), the library was applied to the bioanalyzer (Agilent Technologies, CA, USA) to check the size distribution and quantified using the KAPA library quantification kit (Kapa Biosystems, MA, USA) on an ABI Prism 7700 sequence detection system (Life Technologies, CA, USA). Sequencing was carried out using the HiSeq. 2500 instrument (Illumina) with 2 × 150-bp paired-end reads and dual 8-bp barcode sequences. Sequence reads were aligned to the human reference sequence (hg19) by Burrows-Wheeler Aligner (ver. 0.7.9a) and then applied to the RealignerTargetCreator and IndelRealigner tools using GATK (ver. 3.2.2) for each bam file. The quality control was performed as previously described 45 . We selected variants that clearly showed three peaks corresponding to three genotypes and two peaks if they were considered rare variants by visual inspection.