Introduction

Scoliosis, characterized by a lateral curvature of the spine, represents a multifactorial disorder highly influenced by genetic factors. It may manifest as part of the phenotypic spectrum of genetic syndromes such as Marfan syndrome, spinal muscular atrophy and osteogenesis imperfecta. However, in approximately 80% of patients, scoliosis appears to be idiopathic scoliosis (IS), presenting without an identifiable underlying aetiology, with adolescent idiopathic scoliosis (AIS) forming the largest subset when classified by age of onset.1 2 Despite the absence of vertebral malformations, IS can develop into severe spinal deformities resulting in dysmorphic feature, low back pain, incapacity to work or even paralysis if not treated.3

To date, the molecular aetiology of AIS remains largely unclear. Prior twin studies have estimated AIS heritability to exceed 85%,4 5 underscoring a pronounced genetic predisposition. Genome-wide association studies (GWAS) have discovered numerous significant AIS susceptibility loci near several matrisome-related genes including LBX1,6 GPR126,7 BNC2,8 PAX1,9 SLC39A8 10 and SOX9,11 along with other novel candidate genes11 12 for subsequent functional validation. However, these SNPs account for only approximately 5% of AIS heritability, and the application of the polygenic risk score (PRS) merely elevates this figure to 11%,13 suggesting the inevitable contribution of rare variants to AIS genetic architecture. Individual large familial studies have pinpointed rare variants in HSPG2,14 POC5 15 and CELSR2R 16 with co-segregation patterns, as putative pathogenic variants for AIS. Rare-variant burden analyses have demonstrated mutational enrichment of FBN1, FBN2 17 and musculoskeletal collagen genes18 in severe AIS cases, although these results were restricted to European ancestry.

In this study, we performed next-generation sequencing (NGS) for a Han Chinese cohort with severe AIS to analyse whole exome-wide rare variants, with the objective of unveiling new elements of the ‘AIS etiology puzzle’. Gene-based and pathway-based burden analyses were undertaken to identify significant genes and biological processes implicated in AIS pathogenesis.

MethodsStudy subject

Blood samples of 368 unrelated Han Chinese individuals with severe AIS (Cobb angle ≥40°, age of onset between 10 and 18 years) who underwent spinal surgery at Peking Union Medical College Hospital (PUMCH) under the Deciphering Disorders Involving Scoliosis & COmorbidities (DISCO, http://discostudy.org/) study group, from November 2016 to November 2021 for correction of scoliosis or kyphosis were consecutively recruited. The diagnosis of AIS was strictly confirmed by both an orthopaedic surgeon and a radiologist using series of clinical and radiological techniques (online supplemental methods). More clinical features of the cases are listed in online supplemental table 1.

The control individuals were aggregated from individuals who underwent exome sequencing at the Medical Research Center, PUMCH for clinical or research purposes without IS or any other skeletal malformation. The resultant control dataset constituted 3794 unrelated individuals, maintaining a balanced sex distribution.

Next-generation sequencing

NGS, including whole exome sequencing (WES) or whole genome sequencing (WGS), was performed on peripheral blood DNA samples from the study participants. The number of case and control samples that underwent WES or WGS are provided in online supplemental table 2.

WES was performed in paired-end mode on the Illumina HiSeq 2000 platform, with an average depth of coverage of 70X, while PCR-free WGS was performed on the Illumina HiSeq X-Ten or NovaSeq platform with an average depth of coverage of 35X (online supplemental methods).

Variant calling and quality control

According to Peking Union Medical College Hospital pipeline,19 20 the raw sequencing reads were mapped to the GRCh37 human reference genome using Picard (broadinstitute.github.io/picard/) and Burrow-Wheeler Aligner (BWA).21 To optimally preserve the genomic information, the mixed WES and WGS data were initially processed as a unified WGS dataset. Using the Sentieon’s Haplotyper, the resulting variant call files (VCFs) for each sample were then converted to genomic VCFs (gVCFs) and the individual gVCFs were subsequently combined and jointly called, utilizing a population-based approach to take into account information from all samples concurrently to detect low-frequency variants that may be missed when calling independently. Single-nucleotide variants and small insertions/duplications (indels) were called using the DNA sequencing module of the Sentieon software.22 Variant-level and sample-level quality control (QC) were performed on the joint call-set (online supplemental methods). To note, variants with a missing rate >10% across the whole cohort were filtered out, mitigating the possible biases between WES and WGS data, and ensuring a more reliable and consistent dataset for analysis.

After QC filtering, a total of 1.44 million genetic variants and 362 case samples and 3740 control samples were included for subsequent analysis.

Gene-based rare variant association analysis

We classified each variant into different mask levels (from 1 to 6) and assigned weight values (ranging from 1 to 0) to each mask level according to the variant type and bioinformatic prediction results (online supplemental methods, online supplemental table 3). The mutational burden of a given gene in each individual was assigned as the maximum weight value among all ultra-rare variants carried by the individual (if any). The ACAT package23 was used to perform the weighted mutational burden test for each gene, with the minor allelic count threshold for each gene being included set to 5 (weight≠0). P values were adjusted for multiple testing using the Benjamini-Hochberg false discovery rate (FDR) (online supplemental methods).

Gene-set enrichment analysis for lactate transporters

Based on the gene-based analysis result, we performed enrichment analysis for six established lactate transporter genes set including proton-coupled and sodium-dependent monocarboxylates transporters (MCTs and SMCTs) based on the function and structure according to the literatures24 25 (online supplemental table 7). We used the statistic derived from the gene-based association result and matching background genes to the specified gene set via a k-dimensional tree. A one-sided Wilcoxon test was implemented to evaluate whether the given gene set had a lower p value than expected.

Pathway-based rare variant burden analysis

We tested the aggregated mutational burden of gene sets related to biological processes, molecular pathways or disease phenotypes. The associated gene sets were retrieved from The Molecular Signatures Database (MSigDB, https://www.gsea-msigdb.org/gsea/msigdb/, accession date: 15 August 2022), collecting curated canonical pathways and ontology gene sets from diverse resources including Kyoto Encyclopedia of Genes and Genomes (KEGG, https://www.genome.jp/kegg/), Gene Ontology (GO, http://geneontology.org/) biological process, Reactome (https://reactome.org/), The Pathway Interaction Database (PID, http://pid.nci.nih.gov), BioCarta (http://www.biocarta.com/), WikiPathways (https://www.wikipathways.org/) and hierarchical terms under ‘Abnormality of the skeletal system’ (HP:0000924) in The Human Phenotype Ontology (HPO, https://hpo.jax.org/). All of these resources were accessed on 15 August 2022. Moreover, random gene sets were generated for permutation.

The algorithm for calculating the gene-set burden was referred to a latest study.26 For a candidate gene set, the mutational weight (same as the gene-based test) of each gene was summed (per individual) into a burden score, which was then used to fit a logistic regression model to predict the case-control status. ORs and p values are presented as measures of enrichment in the tested gene sets. The Benjamini-Hochberg FDR procedure was used to adjust for multiple tests (online supplemental methods).

RNA sequencing and preprocessing

RNA data were obtained from muscle tissues of 14 patients with AIS from the initial cohort and another 15 patients diagnosed as congenital vertebral malformation (CVM) who also underwent corrective surgery for scoliosis at PUMCH. The gender distribution aligns with the primary cohort and prior studies, showing a female-to-male ratio of 13:1 for AIS and 8:7 for the control group (online supplemental table 8). The patients were selected based on their diagnosis and the availability of the tissue samples.

The RNA sequencing raw reads were processed through in-house QC procedures and aligned to the human reference genome GRCh37 using the Hisat227 with default parameters (online supplemental methods). The gene expression levels were then quantified using the htseq-count28 programme with Ensembl gene annotations (release 106), generating a count matrix with raw expression values.

Differential expression analysis and pathway enrichment

Differential gene expression analysis was performed on the raw count matrix of 14 patients with AIS and 15 CVM controls using the R package DESeq2.29 We imported the count matrix data and generate dataset for differential expression analysis using ‘DESeqDataSetFromHTSeqCount’ function. To mitigate the impact of noise in the data, genes with <10 counts were filtered out as these genes are more likely to have unreliable expression estimates. Genes that showed more than a 1.5-fold change in expression when compared with controls, either upregulated or downregulated, were retained for further analysis. Pathway enrichment analysis was performed using the clusterProfiler (V.4.0)30 with the predefined gene sets in GO and KEGG database. All the original p values were adjusted for multiple testing using the Benjamini-Hochberg FDR procedure. Similar analyses were performed between the AIS subgroups based on their variant carrying status to explore the integrated consequences of theses variants.

ResultsGene-based rare variant analysis identifies SLC16A8 as a novel candidate gene for severe AIS

We recruited a total of 368 patients with severe AIS, defined as a Cobb angle of ≥40° or a history of surgical treatment. Genome sequencing or exome sequencing was performed on all subjects, as well as on 3794 unrelated control samples. To increase the power of the association and decrease the possibility of latent disease in the control population, the study was limited to patients with severe AIS with an average Cobb angle of 51.6°. Cases in this cohort had a female bias of 84.5%, which aligns with previous research suggesting that the female-to-male ratio of affected individuals with severe malformation may be as high as 10:1. The study population consisted of individuals of Han Chinese ancestry with distant relatedness, as confirmed by genetic estimation (see ‘Methods’ section).

To investigate the role of rare genetic variants in AIS susceptibility, a weighted gene-based burden test was implemented on rare variants. The weight value assigned to each variant was determined based on its influence (see ‘Methods’ section).

As a result, MOK showed the strongest association with severe AIS compared with controls (p=8.86E-07, FDR=0.005, OR=0.90) (table 1, figure 1A). MOK encodes a ubiquitously expressed protein kinase of the MAPK superfamily, involved in cell growth, neuronal development and immune response. Homologous mok-knockout mice exhibited neurological phenotypes such as abnormal behaviour and decreased grip strength, according to the Mouse Genome Informatics database. The identification of MOK aligns with recent studies,31 32 suggesting a potential role for neuronal development and function in the pathogenesis of AIS. However, all likely gene-disrupting variants identified in AIS cases fall outside the critical protein kinase domain and no MOK rare variant carrier exhibited neurological phenotype, attenuating the likelihood of MOK as a monogenic causal gene (online supplemental figure 1).

Table 1

Gene-based rare variant association analysis results

Figure 1

Gene-based rare variant association analysis identifies SLC16A8 as a novel candidate gene for adolescent idiopathic scoliosis (AIS). (A) Quantile-quantile (Q-Q) plots of gene-based rare variant association analysis. The y-axis represents the observed data, while the x-axis represents the expected values of a normal distribution. Q-Q plots were generated for all rare non-synonymous variants (left) to identify disease-associated genes, and for synonymous variants (right) as a background control. The plotted points show the deviation of the observed distribution from the expected normal distribution. (B) Mutational spectrum of deleterious missense (D-mis) and likely gene-disrupting (LGD) variants in SLC16A8 identified from AIS cases. MFS, major facilitator superfamily domain.

Second to the top signal, we identified SLC16A8 as a novel candidate gene associated with severe AIS (p=3.08E-06, FDR=0.009, OR=7.42) (table 1, figure 1A). SLC16A8 encodes a proton-coupled lactate transporter that is preferentially expressed in the retinal pigment and choroid plexus of the eye,33 and is mildly expressed in various tissue types according to GTEx data (online supplemental figure 2). As a member of the MCT family, SLC16A8 displays high affinity for lactate and shares common substrates with other MCTs, including pyruvate, ketone bodies acetoacetate, β-hydroxybutyrate and acetate, which can potentially regulate carbohydrate metabolism. Additionally, SLC16A8 was implicated to be the major MCT isoform responsible for efflux of glycolytically derived lactic acid from white skeletal muscle.34

All the deleterious variants identified in our patients with AIS, including one nonsense variant (p.Glu231Ter), one frameshift variant (p.Arg235GlyfsTer28), one splicing variant (c.214+1G>C) and two predicted damaging missense variants (p.Gly296Asp, p.Leu134Pro), were located in the major facilitator superfamily domain, which is critical for the symport of lactate and protons (figure 1B). The bi-allelic variant in the same location (c.214+1G>C) has been reported to result in the absence of the SLC16A8 protein and cause a deficit in transepithelial lactate transport in retinal pigment epithelial cells, thereby increasing the susceptibility to macular degeneration,35 confirming the pathogenicity of the splicing variant. Notably, within our cohort, a significant proportion (80%, 4/5) of patients harbouring deleterious SLC16A8 variants developed myopia from an early age prior to the onset of scoliosis, with 75% of these cases classified as high myopia (exceeding 600 degrees) (table 2), indicating a putative association between SLC16A8 and the development of both severe scoliosis and high myopia.

Table 2

Clinical information of patients with rare deleterious variants in SLC16A8 and corresponding variant annotation

Furthermore, to evaluate the general contribution of lactate transport to the pathogenesis of AIS, we performed one-sided Wilcoxon test based on the gene-based burden result rank for established lactate transporter genes (see ‘Methods’ section). Our analysis revealed a noteworthy enrichment of rare variants in lactate transporter genes compared with the background level, as indicated by a significant p value of 0.01, provide further evidence for the potential involvement of lactate transporters in the genetic susceptibility to AIS.

Rare variants in carbohydrate metabolism pathway are enriched in severe AIS

To assess the cumulative effect of genetic variants on different genes functioning within the relevant biological pathways, we conducted a pathway-based gene-set burden analysis.

The result showed a prominent enrichment of rare variants in several pathway involved in carbohydrate metabolism, including the galactose metabolic process as the top signal (Gene Ontology Biological Process: GO:0006012, p=2.93E-05), followed by amino sugar and nucleotide sugar metabolism (KEGG pathway: map00520, p=4.04E-05), carbohydrate-derived catabolic process (GO:1901136, p=6.46E-05) and deoxyribose phosphate catabolic process (GO:0046386, p=7.06E-05), galactose catabolic process (GO:0019388, p=2.34E-04) (figure 2A, figure 2B, online supplemental table 4).

Figure 2

Pathway-based gene-set burden analysis implicates carbohydrate metabolism in the pathogenesis of adolescent idiopathic scoliosis (AIS). (A) Manhattan plot of pathway-based rare variant association. The y-axis represents the −log10 p values for pathways, while the x-axis represents their resource. The data sources selected for the analyses were BioCarta, Gene Ontology Biological Process (GOBP), Kyoto Encyclopedia of Genes and Genomes (KEGG), The Pathway Interaction Database (PID), Reactome and WikiPathways (WP). The red line indicates the threshold for statistical significance. The plot shows a peak for carbohydrate metabolism, indicating significant association of rare variants in this pathway with AIS. (B) Forest plot of the top 10 pathways of pathway-based rare variant association result. The squares represent the ORs for each pathway, with the horizontal lines representing the corresponding 95% CIs. (C) Histogram of the distribution of non-synonymous rare variants in cases (n=362) and controls (n=3740) in genes involved in carbohydrate metabolism. The proportions of variants for each variant mask level are shown above the bar. (D) Bar chart showing the distribution of patients with different variant numbers in the cases. The number of patients harbouring one, two, three or more variants are shown.

In our cohort, 48.5% of the patients with AIS (176/362) exhibited non-synonymous rare variants in 200 genes involved in these candidate carbohydrate metabolism pathways, a significantly higher proportion than that observed in the in-house controls (39.1%; 1466/3740) (p=1.11E-05, Fisher’s exact test), with higher proportion of deleterious variants according to the mask level (figure 2C, online supplemental table 5). However, unlike the observations in SLC16A8, we found that variants in these carbohydrate metabolism genes mostly exhibited variants with moderate or low impact and did not show substantial loss-of-function variants (figure 2C). Meanwhile, we noted that 80.7% of the patients (142/176) carried more than one rare variant per individual (figure 2D). Additionally, the ORs of each carbohydrate metabolism pathway were also not as high as those of other significant pathways (figure 2B).

These significant pathways, which are derived from glucose metabolism, are interdependent and play crucial roles in energy supply, synthesis of glycosaminoglycans, glycosylation of proteins and lipids and generation of derivatives for cellular signalling (figure 3). Of particular significance is the galactose metabolic process pathway, which mainly converts galactose into glucose-6-phosphate, a form of glucose that can be used for energy production or stored as glycogen for later use. Among the different components of the galactose metabolic process, the galactose catabolic processes exhibited an independent significance and a comparative effect size (OR=2.97) (figure 2B), highlighting the possible critical role of energy supplement in the association between galactose metabolic process and AIS.

Figure 3

Illustration of the metabolic pathway of carbohydrates and pinpointed candidate genes. The figure depicts part of the carbohydrate metabolic pathway and the biological processes derived from this pathway, highlighting the enzymes and reactions involved in the breakdown of glucose and other sugars, as well as the generation of extracellular matrix (ECM). Candidate genes of top 20 nominal significance, identified through gene-based burden analysis, are emphasised by underlining and a red colour scheme. TCA, tricarboxylic acid.

Additionally, genes critical for carbohydrate pathway function, such as GALE, DERA and LDHD, were also identified as being among the highest-ranking results in the gene-based burden analysis (table 1, figure 3). Dysfunction of these genes could cause imbalances in energy storage, cell growth and proliferation, developmental signalling and extracellular matrix (ECM) production and maintenance, which could potentially contribute to scoliosis development.

Differential expression analysis of muscle tissue highlights the participation of metabolism in AIS development

To validate the preliminary genomic findings and further explore the potential role of metabolism in the pathogenesis of AIS, we performed a differential expression analysis on muscles tissues obtained from 14 patients with AIS. Our choice of muscle tissue for subsequent analyses hinged on the accessibility of the tissue and the tissue-specific expression data gleaned from the Genotype-Tissue Expression project (online supplemental figure 2).

According to the genomic variant result, we divide the AIS cases into two subgroups. Five of these patients carried likely gene-disrupting variants in different carbohydrate metabolism-related genes, a trait absent in the rest of the cohort (online supplemental table 9). As a result, differentially expressed genes (DEGs) of the AIS subgroups demonstrated a consistent enrichment in glucose metabolic process (p=1.93E-04, FDR=0.12), complemented by additional enrichment in long-chain fatty acid import process (p=2.00E-04, FDR=0.12) (online supplemental figure 3 and table 10), indicating the potential of the candidate variants to disrupt carbohydrate metabolism and lipid metabolism.

As a general control, we collected muscle tissues from 15 patients diagnosed with CVM, characterised by vertebral dysplasia and/or vertebral fusion. Given the distinct vertebral anomalies observed in CVM, we hypothesised that the underlying genetic mechanisms may differ from those involved in AIS. Our analysis highlighted numerous enriched metabolic pathways among the DEGs between the AIS cases and CVM controls. Notably, there was a significant enrichment in the positive regulation of catabolic processes (p=2.95E-06, FDR=0.01) as the top signal. Aside from catabolic processes, we also discerned a noteworthy disparity in the response to nutrient levels between patients with AIS and controls (p=2.37E-05, FDR=0.02) (online supplemental figure 4 and online supplemental table 11), suggesting that patients with AIS may have heightened sensitivity to changes in nutrient availability, leading to a more vulnerable pattern in energy regulation. In terms of enriched cellular components, the mitochondrial matrix emerged with the highest significance (p=8.12E-05, FDR=0.01), which is a critical site for metabolism, including the tricarboxylic acid cycle (figure 3, online supplemental table 11).

Distinct clinical manifestation of AIS patients with rare variants in carbohydrate metabolism pathways

To determine whether patients with AIS carrying rare variants in carbohydrate metabolism pathways show distinct clinical features, we conducted a comparative analysis between patients with AIS with rare variants and without rare variants in any relevant gene. The clinical information of certain patients was missing due to lost follow-up or lack of documentation, especially the body mass index (BMI), which was omitted in each independent comparison.

Our results did not reveal any significant different clinical features between patients with AIS with rare non-synonymous variants in carbohydrate metabolism pathways and those without rare variant in relevant genes (online supplemental table 6). However, when we restricted our comparison analysis to deleterious rare variants (mask level 1–3, online supplemental table 3) in the pathways, we observed a significantly larger Cobb angle in patients with AIS with deleterious rare variants compared with those without any rare variant in the carbohydrate metabolism pathway genes (p=0.015, one-sided two-sample t-test). Similar larger curve angle was observed in patients with rare deleterious variants in the galactose metabolism pathways (p=0.0379, one-sided two-sample t-test) (online supplemental table 6). These finding suggest that abnormal energy metabolism is associated with a high risk of progression into more severe scoliosis.

Despite the general low BMI in patients with AIS (mean: 18.8 kg/m2, upper quartile: 16.9 kg/m2, lower quartile: 20.04 kg/m2), consistent to previous studies,36 there was no significant difference in BMI between the two patient groups, except a tendency towards lower BMI in patients with AIS who carried deleterious variants in galactose metabolism pathways (p=0.12) (online supplemental table 6).

Interplay of GWAS loci and carbohydrate metabolism genes revealed the genetic complexity underlying AIS

To elucidate the interplay between carbohydrate metabolism mechanism and common variants, we further investigate the impact of previously reported AIS-associated GWAS loci.

Initially, we implemented Functional Mapping and Annotation of GWAS37 to perform SNP-to-gene analysis including position mapping, expression quantitative trait loci (eQTL) mapping and three-dimensional chromatin interaction mapping strategies, locating 282 genes underlying all the reported significant GWAS loci (online supplemental table 12). On the 282 prioritised genes, we identified five genes, DERA, GALT, HPSE, FUCA2 and PGM3, that are involved in distinct carbohydrate metabolism pathways, Notably, DERA and GALT are in our gene-based burden results, ranking within the top 20 (table 1), suggesting the substantial and potentially inevitable role of carbohydrate metabolism in AIS pathogenesis, mediated through both rare and common genetic variants.

Subsequently, we conducted a PRS analysis on 1663 WGS samples of our cohort using PRScise-238 (online supplemental methods, online supplemental table 13). Our result did not reveal a significant difference between the PRS distributions of patients with AIS and the control group (online supplemental figure 5A), consistent with the understanding that the best-performing PRS model for AIS only accounts for approximately 10% of the heritability,13 which was anticipated an even lesser contribution in our severe AIS cases. However, a notable leftward deviation was observed in the PRS distribution of patients with AIS carrying the SLC16A8 variants (online supplemental figure 5B), indicating the significant effect of these SLC16A8 variants.

We additionally investigated the inheritance pattern of a patient carrying SLC16A8 deleterious variants (p.Gly296Asp) who also exhibited an extremely high PRS score (1.6133), falling within the top 15% of the PRS distribution. The SLC16A8 variant was validated to inherit from her father, who was unaffected by scoliosis (online supplemental figure 6). Interestingly, the PRS scores for the parents were 0.7572 and 0.6979, placing them in the bottom 15% of PRS distribution. This significant disparity in PRS scores between the affected proband and the unaffected parents provides insights into the mechanism of incomplete penetrance mediated by polygenic factors, suggesting a protective effect against AIS in individuals with lower PRS scores, even when carrying rare deleterious variants.

Discussion

The genetic architecture of AIS is multifaceted and remains enigmatic, characterised by a considerable degree of locus heterogeneity and complexity. In this study, we first performed a comprehensive genetic analysis of a Han Chinese cohort with severe AIS to identify potential genetic factors underlying the development of AIS and uncovered a substantial enrichment of rare variants in carbohydrate metabolism pathways and related genes, providing novel insights into the aetiology of AIS.

The identification of SLC16A8, a lactate transporter gene, as a potential candidate gene associated with severe AIS is noteworthy. Remarkably, all three patients harbouring loss-of-function SLC16A8 variants demonstrated pronounced AIS severity and exhibited high myopia preceding scoliosis onset. As SLC16A8 is dominantly expressed in the retinal pigment epithelium, which is an extremely energy-demanding tissue,39 a dysfunction can lead to lactate accumulation, in the consequence, a suppressed consumption of glucose,40 which may strongly affect the regulation of ocular growth and contribute to the earlier development of myopia compared with scoliosis. Also, it has been recently elucidated that lactate can curb cartilage matrix degradation gene expression in osteoarthritic chondrocytes and stimulate matrix synthesis via HIF1A upregulation.41 Thus, the plausible association of SLC16A8 with both severe scoliosis and myopia underscores the potential involvement of glucose metabolism and retinal pigment and choroid in ocular growth regulation, thereby potentially accelerating myopia onset relative to scoliosis.

Furthermore, our gene-set burden analysis disclosed a significant enrichment of rare variants in lactate transporter genes and carbohydrate metabolism pathways, thereby underlining a probable correlation between carbohydrate metabolism-related biological processes and severe AIS development. Although the impact of individual variants in the carbohydrate metabolism pathway in patients with AIS is mild, the disruption is widespread across relevant genes, resulting in a generalised perturbation and a polygenic burden in broad aspects of carbohydrate-derived metabolism. Of the most significance is the galactose metabolism, which is important for energy production and the biosynthesis of glycoproteins and glycolipids. Recent research has indicated dysregulation of galactose metabolism and glycolysis pathways in progressive patients with AIS compared with non-progressive ones, with aberrant LBX1 expression associating with AIS paraspinal muscle energy metabolic markers.42 The severity of curvature was significantly higher in patients carrying deleterious variants in carbohydrate-derived metabolism or galactose metabolism pathways, suggesting early stage impairments in carbohydrate, especially galactose, metabolism as a possible risk factor for severe curvature development.

Supporting the involvement of metabolism in scoliosis pathogenesis, our independent differential expression analysis of muscles revealed notable differential expression in glucose metabolism as well as lipid import process between patients with AIS harbouring deleterious variants in carbohydrate metabolic pathway genes and those who do not, suggesting that these variants can potentially modulate both carbohydrate and lipid metabolism. Furthermore, the positive regulation of catabolic processes and response to nutrient levels exhibited a significant divergence in patients with AIS compared with CVM controls. Given the central role carbohydrates play as a primary nutrient, dysfunction in carbohydrate metabolism, coupled with a high sensitivity to nutrient availability changes, could induce direct or indirect serious outcomes, such as the emergence of scoliosis without vertebrate anomalies.

Previous studies have demonstrated diverse concepts and hypotheses for the pathogenesis of AIS, including axial rotational instability,43 relative anterior spinal overgrowth and asynchronous spinal neuro-osseous growth,44 hormonal and metabolic factors,45 biomechanical and neuromuscular factor46 and more recently, genetics factors influencing maintenance or composition of the cartilage matrisome.47 Compared with other aspects, the understanding of the interplay between metabolic processes and AIS pathogenesis is lacking. Emerging studies have reported abnormal metabolism phenomena in patients with AIS, such as low BMI, low circulating leptin level and osteopenia.36 48 Moreover, the inhibition of glycolysis in chondrocytic cell lines was reported to lead to hypertrophy-like changes and abnormal ECM production,49 50 suggesting a potential role for energy metabolism in AIS development. Our findings provide a novel and comprehensive elucidation of the association between carbohydrate metabolism and AIS, suggesting that the dysregulation of metabolic processes may contribute to an elevated AIS risk, potentially via aberrant ECM production or other yet-to-be-discovered mechanisms.

In concordance with the observation that patients with AIS manifest altered body composition prior to the clinically detected scoliosis,48 the genetic predisposition impacting carbohydrate metabolism may induce an increased spinal instability in late childhood before AIS onset. During the adolescent growth spurt, the human body experience a surge in nutrient demand, energy utilisation and biosynthesis, necessitating meticulous regulation of glucose-dependent processes. Given the exceptional flexibility of spine and the gravitational stress, an improper growth trajectory may amplify spinal instability, thereby facilitating postural instability and the onset and progression of scoliosis (figure 4). Our comparative analysis, however, did not reveal a significant metabolic manifestation, with only a slight trend towards lower BMI in patients carrying non-synonymous variants in the galactose metabolism pathway. Given the generally low BMI of patients with AIS, a substantial metabolic disruption would be required to further decrease BMI against this already low baseline, which is not expected in our findings. Nevertheless, BMI is not an optimal or sensitive indicator of metabolic status. Future studies could employ more refined and information-rich phenotypic data, such as metabolomics, proteomics or transcriptomics of relevant cell types. An important caveat in interpreting these results is the necessity to discern whether the observed metabolic changes are primary alterations driving the disease or compensatory responses aimed at maintaining homeostasis in the face of the disease.

Figure 4

Illustration of carbohydrate metabolism dysfunction on severe adolescent idiopathic scoliosis (AIS) across life stages. This figure illustrates the evolution and progression of severe AIS, triggered by a generalised dysfunction in carbohydrate metabolism. The ‘Initiation Phase’ marks the early onset of AIS, reflecting the impact of metabolic disturbances. The ‘Progression Phase’ underlines the swift progression of scoliosis during periods of typical growth spurts in preteens and teens. The ‘Challenge Phase’ in adulthood encapsulates the enduring struggle with the disease, marked by increased risks of spinal degeneration, chronic pain, cardiopulmonary complications and diminished quality of life. Environmental factors exert influence throughout all these stages. The figure was created by BioRender.com.