An international meta-analysis confirms the association of BNC2 with adolescent idiopathic scoliosis

Adolescent idiopathic scoliosis (AIS) is a common spinal deformity with the prevalence of approximately 3%. We previously conducted a genome-wide association study (GWAS) using a Japanese cohort and identified a novel locus on chromosome 9p22.2. However, a replication study using multi-population cohorts has not been conducted. To confirm the association of 9p22.2 locus with AIS in multi-ethnic populations, we conducted international meta-analysis using eight cohorts. In total, we analyzed 8,756 cases and 27,822 controls. The analysis showed a convincing evidence of association between rs3904778 and AIS. Seven out of eight cohorts had significant P value, and remaining one cohort also had the same trend as the seven. The combined P was 3.28 × 10−18 (odds ratio = 1.19, 95% confidence interval = 1.14–1.24). In silico analyses suggested that BNC2 is the AIS susceptibility gene in this locus.

Sex-stratified association. AIS has an ample clinical evidence of sexual dimorphism 13 . In our previous study, we investigated BNC2 expression in a variety of human tissues and found that BNC2 expression is highest in uterus, suggesting its sex-related biological function 6 . Therefore, we performed sex-stratified analyses to determine whether a genetic difference existed between male and female. We could obtain sex information for both cases and controls in five cohorts. We could obtain 6,266 cases and 15,292 controls in the female-only analysis, and 485 cases and 10,490 controls in the male-only analysis (Supplementary Tables 1 and 2). In both sexes, we could not find genome-wide level significant association (P = 5 × 10 −8 ); particularly in male, the P value did not even reach to the nominal association level (P = 5 × 10 −2 ) (Tables 2 and 3). However, the ORs were similar between male and female, which were similar to that in the analysis disregarding the sex (Table 1).
Fine mapping. The landmark SNP rs3904778 is located in intron 3 of BNC2, and BNC2 is the only gene contained within the linkage disequilibrium (LD) block (r 2 > 0.8) represented by rs3904778. The topologically associated domain (TAD) is the partition of the genome that represents a regulatory unit within which enhancers and promoters can interact 14 . To identify the candidate susceptibility gene in the locus, we evaluated the TAD around the associated SNPs using H1-mesenchymal stem cell. Hi-C data 15 (http://promoter.bx.psu.edu/hi-c/view. php) revealed that BNC2 was the only gene included in the TAD that contained the LD block of the associated SNPs (Fig. 2). The data strongly suggested that BNC2 is the most plausible AIS susceptibility gene at this locus.

Discussion
In the present study, we have performed a meta-analysis for the genetic association of rs3904778 with AIS using more than 36,000 subjects from eight independent multi-ethnic cohorts. To date, no large-scale replication study for the association of the AIS locus has been conducted. Previously, we demonstrated that rs3904778 had significant association with AIS in Japanese and Chinese 6 ; however, no evidence has been reported regarding its association in non-East Asian populations. The present study not only gave solid evidence of association of the locus in additional Chinese cohorts, but also revealed that it had significant association in Caucasian, suggesting the global significance of this AIS locus. Previous lack of association in Caucasian may be due to lack of power because the OR of this locus is about 1.2, suggesting relatively large sample size is optimal for identification. The most significantly associated SNPs are clustered in intron 3 of BNC2. BNC2 is the only gene contained in the LD block of the associated SNPs. TAD containing the LD block only contained BNC2 (Fig. 2). These genome data strongly suggest that BNC2 is the AIS susceptibility gene in the locus. rs10738445 in the locus is in high LD (r 2 = 0.9) with rs3904778. Genevar (Gene Expression Variation) data revealed that the risk allele of the functional SNP in this locus, rs10738445, increased BNC2 expression (p = 0.048) 6 . Our previous in vitro analyses Forest plots for the association of rs3904778 with AIS susceptibility. The odds ratios and 95% confidence intervals were estimated based on the fixed-effect model. The contributing effect from each study is shown by a square with its confidence interval indicated by a horizontal line. Summary: the combined metaanalysis estimate.  revealed that the risk allele of rs10738445 functioned as an enhancer element and caused increased BNC2 expression through the increased binding of a transcription factor, YY1 (Ying-Yang 1) 6 . BNC2 was highly expressed in musculoskeletal tissues such as spinal cord, bone and cartilage 6 . GTEx database also showed similar expression pattern; BNC2 expression was the highest in uterus followed by ovary and nerve. We hypothesized that increased BNC2 expression in these tissues lead to susceptibility of AIS. Actually, the over-expression of Bnc2 in zebrafish caused scoliosis-like deformity 6 .
To gain insight into the sex difference in AIS susceptibility, we examined sex-stratified association of rs3904778. While the association was almost genome-wide significant level in the female-only analysis (6,266 cases and 15,292 controls), no significant association was obtained in the male-only analysis (485 cases and 10,490 controls) (Tables 2, 3). This is most probably due to be lack of power in the male analysis; in the analysis, sample size was small, especially in the case group, which reflected the female prevalence in all ethnic populations 6,7,16 . It is of note that the ORs were similar in both sex-stratified analysis. Further analysis with a sufficient sample size will be necessary for the male AIS study, which would inevitably be an international, mutli-center study.

Methods
Subjects and genotyping. We obtained informed consent from all subjects and/or their parents. The ethics committee of RIKEN approved this study. All experiments were performed in accordance with relevant guidelines and regulations. The datasets generated during the current study are available from the corresponding authors on reasonable request. AIS subjects were diagnosed through clinical and radiological examinations according to the previously described criteria 4,6,9 . The subjects in the Japanese and Nanjing-Chinese cohorts were recruited and genotyped as previously described 4,6,9 . The detail of beadchip information, quality control and statistical analysis were also previously described 4,6,9 . The details of additional studies (Guangzhou, Hong Kong, Beijing, USA, and Scandinavia studies) were described as below.
Guangzhou study. We recruited AIS subjects from the First Affiliated Hospital and Sun Yat-sen Memorial Hospital of Sun Yat-sen University as previously described 12 . We recruited control subjects from individuals who received scoliosis screening at middle and primary schools in Guangzhou and fracture patients selected from the First Affiliated Hospital and Sun Yat-sen Memorial Hospital of Sun Yat-sen University. Orthopedic surgeons evaluated these subjects with Adam's forward bending test and scoliometers to screen scoliosis. We extracted genomic DNA from blood using DNA Blood Mini-kit (Tiangen Biotech, Beijing, China). The 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. We recruited AIS subjects from the Duchess of Kent Children's Hospital in Hong Kong with previously described inclusion criteria 11 . We randomly selected control subjects from the subjects recruited for the Genetic Study of Degenerative Disc Disease project 17 . We confirmed control subjects did not have scoliosis by MRI examination of the spine. We extracted genomic DNA from peripheral blood lymphocytes using standard procedures. We used the PCR-based invader assay (Third Wave Technologies, WI, USA) for genotyping.
Beijing Study. We recruited AIS subjects from Peking Union Medical College Hospital. All subjects underwent clinical and radiologic examination and expert spinal surgeons evaluated scoliosis. We extracted genomic DNA from peripheral blood using QIAamp DNA Blood Mini Kit (Qiagen, Hilden, Germany). We used the MassARRAY system (Agena Bioscience, San Diego, CA, USA) for genotyping. USA study. We recruited AIS subjects at Texas Scottish Rite Hospital for Children as previously described 7 and used the Illumina HumanCoreExome Beadchip array for genotyping. For controls, we utilized a single dataset of individuals downloaded from the dbGaP web site (http://www.ncbi.nlm.nih.gov/sites/entrez?db = gap) from Geisinger Health System-MyCode, eMERGE III Exome Chip Study under phs000957.v1.p1 (https://www.ncbi. nlm.nih.gov/projects/gap/cgi-bin/study.cgi?study_id = phs000957.v1.p1). 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 subjects having 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. 15 ). Duplicated or related individuals were removed as previously described 18 . We used principal component analysis (PCA) 19 on the merged data projected onto HapMap3 samples to correct possible stratification 20 . After quality controls, 9,312 subjects (1,360 AIS patients and 7,952 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 imputed genotypes for the region around rs3904778 using minimac3 21 with the 1000G-Phase3.V.5 reference panel accoding to the instructions of the software (http://genome.sph.umich. edu/wiki/Minimac3_Imputation_Cookbook).
Scandinavia study. We recruited AIS subjects from six hospitals in Sweden and one in Denmark as with previously described inclusion criteria [22][23][24][25] . We recruited control subjects from the Osteoporosis Prospective Risk Assessment cohort and PEAK-25 cohort 26,27 . Dual-energy X-ray absorptiometry (DXA) scan was performed in both cohorts and subjects with any sign of scoliosis on DXA were excluded. We extracted genomic DNA from blood or saliva using the QIAamp 96 DNA Blood Kit and Autopure LS system (Qiagen, Hilden, Germany). We used iPLEX Gold chemistry and MassARRAY system (Agena Bioscience, CA, USA) for genotyping. Two persons checked genotype calls using the MassARRAY Typer v4.0 Software (Agena Bioscience).

Statistical analysis. The association between rs3904778 and AIS in each study was evaluated by the
Cochrane-Armitage trend test aside from the Japanese 1 and USA studies since rs3904778 was an imputated SNP in the two studies. The Japanese 1 study was analyzed as previously described 6 . For the USA study, Mach2dat software 28 was used to test the imputed allele dosages of rs3904778 by logistic regression with gender and principal components as covariates. 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/) 29 . The heterogeneity among studies was tested using Cochran's Q test based upon inverse variance weights using METAL.