Genetic polymorphisms of GPR126 are functionally associated with PUMC classifications of adolescent idiopathic scoliosis in a Northern Han population

Abstract GPR126 has been identified to be associated with AIS (Adolescent Idiopathic Scoliosis) in different populations, but data on the northern Chinese population are unavailable. Additionally, it is important to know the exact clinical phenotypes associated with specific genetic polymorphisms. Fourteen SNP (single nucleotide polymorphism) loci in GPR126 were genotyped in 480 northern Chinese Han AIS patients and 841 controls. These patients were classified into three types based on the PUMC classification system. Luciferase assays were used to investigate their regulation of GPR126 transcription activity. Combined and stratified genotype–phenotype association analyses were conducted. The alleles rs225694, rs7774095 and rs2294773 were significantly associated with AIS (P = 0.021, 0.048 and 0.023, respectively). rs225694 and rs7774095 potentially have regulatory functions for the GRP126 gene. Correlation analysis revealed that allele A of rs225694 was a risk allele only for PUMC type II AIS (P = 0.036) and allele G of rs2294773 was a risk allele only for PUMC type I AIS (P = 0.018). In summary, rs225694, rs7774095 and rs2294773 are significantly associated with disease in northern Chinese Han AIS patients. The SNPs rs225694 and rs2294773 are associated with different AIS PUMC classifications.

Background AIS is the most common spinal deformity, affecting approximately 1-3% of children throughout the world [1,2]. There are many classification systems, including the King classification [3], Lenke classification [4] and the Peking Union Medical College (PUMC) classification systems [5]. According to the PUMC classification, AIS can be classified into single (I), double (II), or triple (III) curve types. The PUMC classification is a three-dimensional classification system that is relatively simple to use and has corresponding recommended surgical approaches [5]. Despite decades of research on AIS, its exact aetiological and pathological causes remain unclear.
Previous studies have suggested that genetic polymorphisms play a pivotal role in the pathogenesis of AIS [6,7]. Although a recent genomewide association study (GWAS) and other studies have identified several susceptible gene variants, such as GPR126 [8], LBX1 [9], PAX1 [10], BCN2 [11] and PTK7 [12], few of these variants have been replicated in different populations, and even fewer have been analysed with regard to different AIS subtypes.
GPR126 (encoding G protein-coupled receptor 126), also known as ADGRG6, located at 6q24.1, is 144,348 bp in length and consists of 25 exons. According to Monk et al. [13], gpr126 is essential for mouse survival due to Schwann cell myelination. Courtney et al. also found [14] that a gpr126 deletion in cartilage could lead to idiopathic scoliosis and pectus excavatum in mice. It has been shown that GPR126 gene SNPs are associated with AIS. Kou et al. [8] identified rs6570507 as a susceptibility locus with AIS in 1,819 Japanese cases and 25,939 controls and found that the developing spine of a gpr126 knockdown zebrafish model had delayed ossification. Subsequently, Xu et al. [15] found this locus and two other loci (rs7774095 and rs7755109) in 352 southern Chinese cases and 149 controls. Recently, Qin et al. [16] found another functional locus (rs9403380) that regulates GPR126 expression in the paraspinal muscles of southern Chinese AIS patients.
As far as we know, there has not been an association analysis of the GPR126 gene with AIS in the northern Chinese Han population. It will be of great significance to identify the specific subgroups of AIS that are associated with GPR126 genetic polymorphisms. We suggested that GPR126 risk SNPs are potentially associated with the AIS PUMC classification system. Therefore, we enrolled a cohort from the northern Chinese Han population and corresponding controls. SNPs in GRP126 were genotyped with subsequent phenotypic association analysis based on their PUMC classification.

Materials and methods Participants
A total of 480 patients diagnosed with AIS and 841 healthy in-house controls from a northern Chinese Han population at Peking Union Medical College Hospital were enrolled between July 2011 and August 2016. The inclusion and exclusion criteria were as follows: Inclusion criteria (1) Diagnosed as AIS with a Cobb angle > 20 degrees.
(3) Origin was from the northern Chinese Han population as identified by a Native Place questionnaire. (4) Having complete imaging data, including X-ray, three-dimensional imaging of the spine CT or spinal MRI.
Exclusion criteria (1) Scoliosis secondary to a known aetiology, including congenital scoliosis, neuromuscular scoliosis, scoliosis secondary to skeletal dysplasia or connective tissue abnormalities.
(2) Origin was not from the northern Chinese Han population.
All patients diagnosed as AIS were classified into three types based on the PUMC classification system [5] by at least two experienced orthopaedic surgeons. Written informed consent was obtained from all participants or their parents. The Ethics Committee of Peking Union Medical College Hospital, Chinese Academy of Medical Sciences approved this study.

SNP selection and genotyping
Fourteen candidate SNPs around GPR126 were selected by the following criteria from a literature review and the NCBI SNP database (www. ncbi.nlm.nih.gov/SNP; Table S1): (i) SNP with a minor allele frequency (MAF) above 5%; (ii) Tag SNPs were preferred; (iii) reported SNPs associated with idiopathic scoliosis; and (iv) potential functional SNPs predicted by HaploReg (v4.1) [17].
Genomic DNA was extracted from the peripheral blood of each participant using DNeasy Blood & Tissue Kits (QIAGEN, Eastwin Scientific, Inc. Beijing, China). SNP genotyping was performed using the MassArray system from Sequenom. All of the details were the same as in our previous studies [18][19][20].

Luciferase assay
We used luciferase assays to test the influence of significantly associated SNPs on transcription. We cloned oligonucleotides around each SNP (Table S2) and cloned them into a pGL3 promoter luciferase reporter vector (Promega, Madison, WI, USA). Human embryonic kidney 293a cells (HEK293a) and HeLa cells were transfected with the reporter vector by lip3000 (Invitrogen, Carlsbad, CA, USA). After a 48-hr incubation, cells were harvested and luciferase activity was measured using a Dual-Luciferase Reporter Assay System (Promega, Madison, WI, USA).

Data analysis
Hardy-Weinberg equilibrium (HWE) tests and primary analyses, including allelic, genotypic, haplotypic analyses and linear regression controlling for sex were conducted using PLINK software (v1.07, 2009 Shaun Purcell) and the Haploview program (version 4.2, Broad Institute of MIT and Harvard, Cambridge, MA, USA). The baseline characteristics of AIS participants presenting different PUMC types were compared through SPSS software (16.0 version, SPSS Inc., Chicago, IL, USA). Intergroup differences were assessed with the chi-square test for categorical variables and one-way ANOVA for continuous normally distributed variables.

Baseline characteristics of participants
Seventy-eight males and 402 females were enrolled in the case group, and 483 males and 358 females were enrolled in the control group (Table 1). All patients included in our study were classified into three different PUMC AIS types based on their radiology data (Fig. S1). There were no significant differences in the main Cobb angle, age or gender between each of the PUMC groups (Table 1).

Genotyping and Hardy-Weinberg equilibrium test
Among all participants, 14 SNPs in or around GPR126 were genotyped. rs6929442 was excluded from subsequent analyses because of the low call rate (Table S1). The other 13 SNPs were successfully genotyped with a minimum call rate of 98%. None of the 13 SNPs deviated from HWE and, thus, were subjected to subsequent analyses [21].

Association analysis
Of the 13 SNPs, rs225694, rs7774095 and rs2294773 presented significant allelic differences between the case and control groups. This result was demonstrated with LocusZoom, which showed that the alleles of rs225694, rs2294773 and rs7774095 had the highest association signal with AIS (http://csg.sph.umich.edu/locuszoom/) (Fig. 1). Controlling for sex, allele A of rs225694, allele A of rs774095 and allele G of rs2294773 had risk effects based on the different test models (P = 0.021, 0.048 and 0.023, respectively; Table 2).

Haplotypic analysis
We constructed haplotypes based on genotype data from thirteen SNPs using Haploview software (version 4.2). The pairwise linkage disequilibrium r-square values between SNPs and the LD plots are presented in Figure S2. According to the Four Gamete Rule, we identified two haplotype blocks: rs7774095 and rs6570507 were in one block and rs35699755, rs17280293, rs11155242, rs2143390, rs4896582, rs7755109, rs2294773, rs2294775 and rs3748069 were in another block. There was no significant association between these two haplotype blocks and AIS when controlling for sex (Table 3).

Luciferase assay
Among the three SNPs with positive results in association analyses, rs225694 and rs7774095 had potential enhancer activity. The luciferase activity of the risk allele and no-risk allele was compared. The construct containing risk allele A of rs225694 had a 1.97-fold higher level of luciferase activity than non-risk allele C (P = 0.007). The construct containing risk allele A of rs7774095 had a 0.86-fold lower level of luciferase activity than the non-risk allele C (P = 0.015; shown in Fig. 2). Experiments using the HeLa cell line showed similar results (Fig. S3).

PUMC subgroup analysis
We conducted PUMC subgroup analysis in the three SNPs with positive association results. The analysis indicated that the genotypic association of rs225694 had positive results for the PUMC type II subgroups (P = 0.036). rs2294773 presented a positive result for the PUMC type I group (P = 0.018), and the other subgroups showed no significant difference (Table 4).

Discussion
In this study, we investigated the association of fourteen SNP loci in the GPR126 gene with AIS in a northern Chinese Han population. We found a significant association signal of rs7774095 with AIS. Moreover, we identified two novel SNP loci (rs225694 and rs2294773) that not only contributed to predisposition but were also associated with different types of AIS.
rs225694 was first identified to be associated with height [22]. It is located in an intronic area of VTA1, but is in the same topologically associated domain (TAD) as GPR126 (Fig. 3). TAD is the compartment of the 3D chromatin space in which genomes of many bilaterian animals are organized. It is universally recognized that TAD plays an important role in gene regulation [23] and has potential gene-enhancer interactions [24]. Thus, rs225694 is supposed to have enhancer activity for the gene GPR126. Association analysis showed that this allele was most significantly associated with AIS, controlling for sex. Luciferase assays confirmed its enhancer function, including up-regulation of transcription activity. Considering the PUMC classification, rs225694 had a powerful association only with type II AIS. We therefore suggested that the risk allele of rs225694 up-regulated transcription of GPR126 and increased the risk of PUMC type II AIS. rs7774095, an intronic SNP of GPR126, was first reported to be associated with AIS in a Japanese population in 2013 [8]. Other studies in southern Chinese populations replicated the association of this locus [15,16]. However, there was no clear functional evidence provided for this conclusion. We also observed this association and showed that it down-regulates transcription activity. It has been shown that deletion of gpr126 in cartilage leads to an idiopathic scoliosis deformity in a mouse model [14]. However, detailed clinical phenotype analyses found no significant association with any PUMC AIS subgroup, which may be due to the limited sample size. Additional studies with larger sample sizes are needed to identify the relationship between rs7774095 and PUMC AIS subgroups.
Association analyses of rs2294773 also showed significant results with AIS. rs2294773 is located in the 3 0 UTR of GPR126. It is postulated that the 3 0 UTR is a regulatory region that has potential activity in regulating mRNA localization, translation and stability [25]. Different alleles at one locus can allow or eliminate miRNA binding. Song et al. found that mRNA containing the susceptibility allele rs4148941 had an enhanced interaction with miR-513a-5p, which may contribute to the aetiology of Lumbar disc degeneration (LDD) [26]. Zhu et al. found a functional SNP in the 3 0 UTR region associated with bone mineral density (BMD) [27]. As for rs2294773, in silico predictions by three online tools, including TargetScan (http://www.ta rgetscan.org/) [25,28,29], miRanda (http://www.microrna.org/) [30] and mirSNP (http://cmbi.bjmu.edu.cn/mirSNP) [31], showed that there are three microRNAs, hsa-miR-1343, hsa-miR-198 and hsa-miR-30b-3p that can bind to this mRNA. We also performed eQTL analysis in different tissues from the database and found that rs2294773 was associated with GPR126 gene expression in several tissues, including transformed fibroblasts, sun-exposed skin and so on (Table S3, Fig. S4). Although there was no eQTL data for cartilage or bone, this locus may influence GRP126 gene expression and potentially contribute to AIS formation by these above miRNAs. Regarding the PUMC classification, rs2294773 also had a powerful association with only type I AIS. Additional replication studies with larger sample sizes are needed to provide further evidence for this association.
In our study, we identified two haplotype blocks in the candidate GPR126 gene with AIS. In the first block, we found that the A-A Fig. 2 Relative Luciferase Activity for different allele of rs225694 and rs7774095. (A) Transcriptional enhancer activities of rs225694 constructs. The construct containing the risk allele (R) had 1.97-fold higher level relative luciferase activity than the non-risk allele. *P value < 0.05. Error bars indicate stand error. The assay was repeated three times. (B) Transcriptional enhancer activities of rs7774095 constructs. The construct containing the risk allele (R) had 0.86-fold lower level of relative luciferase activity than the non-risk allele. *P value < 0.05. Error bars indicate stand error. The assay was repeated three times. haplotype of rs7774095-rs6570507 was not significantly associated with AIS (P = 0.067), although previous studies had identified an association of rs6570507 with AIS. In contrast to previous studies, the lack of evidence of an association between rs6570507 and AIS may be due to different allele frequencies between geographic populations or the limitation of the sample size studied [8,15,32]. In the second block, we also did not find significant association.
In conclusion, our study found that three SNPs near or in GPR126 were associated with AIS in the northern Chinese Han population. Risk allele A of rs225694 had higher transcription activity, and risk allele A of rs7774095 had lower transcription activity in vitro.
We suggested that both overexpression and underexpression of GPR126 increased AIS susceptibility for different types of AIS and that SNPs at microRNA binding sites potentially contribute to AIS. GPR126 risk SNPs had potential associations with the type of AIS, as defined by the PUMC classification system. Rs225694 was functionally associated with PUMC type II AIS, and rs2294773 was mainly functionally associated with type I AIS, which contribute to our understanding of AIS development. SNPs of GRP126 were phenotypically associated with distinct PUMC types of AIS and might be useful as potential biomarkers for AIS classification and further surgical intervention.       Table S1. Fourteen target SNPs of GRP126 from literature review. Table S2. Oligonucleotide sequence for luciferase assay. Table S3. EQTL analysis result of rs2294773 in different tissues.