Genome-wide association study identifies GAK and KLF12 associated with curve severity of adolescent idiopathic scoliosis

Subjects

This study was conducted under the approval of the Institutional Review Board of Nanjing Drum Tower Hospital (Approval No. 2019-066-01). All participants were informed of the study’s purpose and provided written consent. The current association analysis was composed of the initial genome-wide discovery stage followed by an independent replication stage. Patients who received treatment for scoliosis in our center between 2000 and 2017 were reviewed. The following inclusion criteria were used: (1) female patients diagnosed as AIS through clinical and radiologic examinations; (2) with major thoracic curve; (3) with Cobb angle less than 30 degrees or more than 40 degrees at skeletal maturity. Overall, 620 female AIS patients were enrolled in the discovery stage of GWAS. 323 patients with curve more than 40 degrees were assigned to the severe group and the other 297 patients with curve less than 30 degrees were assigned to the mild group. For the replication stage, an independent cohort of 634 patients with Cobb angle more than 40 degrees and 546 patients with Cobb angle less than 30 degrees were then included. All participants included in this study were from ethnic Han Chinese populations. Informed consent was obtained from the participants or from their guardians.

Genotyping and quality control

Genomic DNA was isolated from peripheral blood leukocytes using a commercial DNA extraction kit (Qiagen, Hilden, Germany). In the discovery phase, genotyping was conducted with the Affymetrix Genome-Wide Human SNP Array 6.0 platform, adhering to the manufacturer’s protocol. Quality control of the raw genotyping data in the discovery phase was performed using PLINK (v1.90), as we previously described (Zhu et al., 2015). Single nucleotide polymorphisms (SNPs) that did not meet the following criteria were excluded: (1) call rates <95%; (2) minor allele frequency <0.05; (3) significant deviation from Hardy-Weinberg equilibrium (P <  1  × 10⁻⁵); and (4) located on non-autosomal chromosomes. For sample quality control, cryptic relatedness was assessed using an identity-by-state method, and individuals with second-degree or closer relationships were excluded. Additionally, individuals with more than 10% missing genotype data were also removed. To detect population outliers and stratification, we used principal component analysis (PCA) implemented in the software package EIGENSTRAT.

Selection of SNPs for the replication

SNPs surpassing a threshold of P <  1.0  × 10⁻⁴ in the discovery stage were selected for further replication. The SNP genotyping in the replication stage was performed using TaqMan SNP Genotyping was conducted using an ABI Step-One-Plus sequence detection system (Applied Biosystems, Foster City, CA). To assess reproducibility, 5% of the samples were randomly chosen as blind duplicates, which demonstrated a 100% concordance. DNA samples with more than 10% missing data were excluded from further analysis.

Genotype imputation

Genotype imputation within a 400kb window surrounding the two replicated loci (rs2061846 and rs7330031) was conducted using MaCH-Admix software (Liu et al., 2012). The reference data for linkage disequilibrium (LD) and haplotypes were sourced from the 1,000 Genomes Project, specifically from the phased CHB and CHS datasets (March 2012 release). Variants with imputation quality R² < 0.30 or minor allele frequency <0.10 were excluded. Logistic regression models assuming an additive genetic effect were used to evaluate the association between imputed variants and AIS curve severity. Regional association plots were generated using LocusZoom (Pruim et al., 2010).

Functional annotation

The regulatory properties of the novel susceptible signals were analyzed using the Chromatin state segmentation in LCL data generated by the ENCODE project. HaploReg Version 4.2 (Ward & Kellis, 2016) and RegulomeDB Version 2.2 (Boyle et al., 2012) were used to explore the annotations of the susceptible loci on haplotype blocks and to determine whether the variants are located in putative transcription factor binding sites or enhancer elements.

Tissue collection and gene expression analysis

Paraspinal muscle samples from the proximal vertebra region were obtained during corrective surgery from 24 AIS patients for analysis of susceptible gene expression. The collected tissue was immediately flash-frozen in liquid nitrogen and stored at −80 °C for subsequent analysis. Total RNA was extracted with Trizol (Thermo Fisher Scientific, Waltham, MA, USA) according to the manufacturer’s protocol. A total of one µg RNA from each sample was reverse-transcripted to cDNA with the HiScript III All-in-one RT SuperMix Perfect for qPCR (Vazyme Biotech, Shanghai, China), followed by qPCR with the HiScript III All-in-one RT SuperMix Perfect for qPCR (Vazyme Biotech, Shanghai, China). Target gene expression values were calculated using the 2^−ΔCt method and normalized to GAPDH expression (Liu et al., 2012). Specific primers were listed in Table S1.

Statistical analysis

PLINK 1.90 was used for general statistical analysis. The Cochran-Armitage trend test was used to calculate the association of each SNP with curve severity in the discovery stage. For the replication stage, the Chi-square test was used to compare the frequencies of genotype and allele between the mild and the severe group. OR values and 95% CIs were calculated from a 2 ×2 allele frequency table. The Manhattan plot was generated using Haploview version 4.2 (Barrett et al., 2004). A quantile–quantile plot generated with R (v2.6) was used to evaluate the potential impact of population stratification (Turner, 2018). Gene expression correlation analysis was performed with Pearson’s correlation test, or Spearman’s correlation test. Statistical significance was set at a P value of 0.05.

Demographic characteristics of the subjects

Clinical characteristics of the subjects in the discovery and replication stage are summarized in Table 1. The mean age of the patients was 17.8 ± 1.9 years, ranging from 16.5 years to 22 years. All patients were followed up until skeletal maturity with a Risser grade of 5. The average curve magnitude was 41.5 ± 16.7 degrees, with a range of 20 to 72 degrees. For subjects included in the gene expression analysis, the mean Cobb angle were 55.5 ± 6.5 degrees, ranging from 45 to 72 degrees.

Association analysis and identification of severity-related loci

As shown in the Q-Q plot (Fig. S1A), the genome inflation factor (λ) was 1.02, indicating a lower possibility of population stratification. PCA analysis indicated that all subjects were genetically well-matched with minimal evidence of population stratification (Fig. S1B). After quality control, 606,254 SNPs were retained in the severe versus mild AIS association analysis.

Fifteen SNPs showed a p-value of less than 1.0  × 10⁻⁴ (Table 2). These SNPs located within or near several genes, including GAK, DST, SMC2/ NIPSNAP3A, KLF12, ST8SIA5-DT and SEZ6L. We also evaluated previously reported AIS progression–associated variants identified in earlier GWASs (Miyake et al., 2013; Nada et al., 2019; Ogura et al., 2015; Ogura et al., 2017), but none reached suggestive significance in our cohort (Table S2). The top-associated SNPs from each locus (rs2061846, rs12200301, rs10820637, rs7330031, rs2469472, and rs738650) were selected for replication in an independent cohort (Fig. 1).

Replication of severity-associated loci and locus-specific imputation

Replication analysis (Table 3) validated the association of two SNPs: rs2061846 in GAK and rs7330031 in KLF12. For rs2061846, the frequency of allele G was significantly higher in the severe AIS group compared to the mild group (35.3% vs. 29.3%, P = 0.001), with an odds ratio (OR) of 1.32 and a 95% confidence interval (CI) [1.11–1.57]. For rs7330031, allele C was more frequent in severe cases than in mild cases (16.9% vs. 12.2%, p = 0.001), with an OR of 1.46 and a 95% CI [1.16–1.85]. The remaining four SNPs did not show significant differences in allele or genotype distribution between groups.

To further refine the regional association signal, locus-specific imputation was performed across ±400 kb of rs2061846 and rs7330031 (Fig. 2). The association patterns of imputed variants closely mirrored those of the directly genotyped SNPs, and no imputed variant exceeded the statistical significance of the two lead SNPs. These results indicate that rs2061846 and rs7330031 represent the primary association signals within their respective LD blocks.

Relationship between expression levels of susceptibility genes and phenotypes in AIS

The two novel susceptibility loci, rs2061846 and rs7330031, are located within intronic regions of the GAK and KLF12 genes, respectively. To investigate their potential roles in AIS pathogenesis, we examined the mRNA expression levels of GAK and KLF12 in paraspinal muscle samples from AIS patients undergoing corrective surgery (n = 24).

As shown in Fig. 3A, KLF12 mRNA levels were significantly negatively correlated with the Cobb angle (r = −0.451, P = 0.027), suggesting a potential association between lower KLF12 expression and increased curve severity. In contrast, GAK expression showed no significant correlation with curve magnitude (r = −0.203, P = 0.342). Considering the critical role of muscle fiber composition in the pathogenesis of AIS, we further evaluated the relationship between KLF12 expression, markers of slow and fast-twitch muscle fibers, and curve severity. As shown in Fig. 3B, KLF12 expression was positively correlated with slow-twitch myofiber related genes, including MYL3 (r = 0.566, P = 0.004), TNNC1 (r = 0.505, P = 0.012) and MYH7 (r = 0.621, P = 0.001). In contrast, Fig. 3C demonstrated that KLF12 expression was inversely correlated with fast-twitch myofiber–associated genes, such as MYL1 (r = −0.557, P = 0.005), TNNC2 (r = −0.578, P = 0.003), and TNNT3 (r = −0.516, P = 0.010).

Functional annotation of the susceptibility loci

Functional annotation of rs2061846 and rs7330031 revealed that both variants lie within intronic regions of GAK and KLF12, respectively, and map to genomic segments enriched for enhancer-associated histone modifications, including H3K36me3 and H3K4me1 (Fig. S2). Motif analysis further indicated that the allelic substitutions at these loci overlap predicted NKX2 and HMX2 transcription factor binding motifs. Consistent with these observations, genome browser inspection shows that both SNPs are positioned within putative regulatory elements characterized by accessible chromatin. As summarized in Table 4, SNPs in high linkage disequilibrium (r² ≥ 0.9) within the same haplotype blocks were also predicted to have regulatory potential. These findings further support the involvement of these loci in transcriptional regulation processes potentially related to AIS development and progression.