Family-based integrative genomics identifies neural–metabolic–skeletal pathways underlying scoliosis susceptibility

WES Data from Scoliosis Families

The WES data of scoliosis families included in this study were obtained from 34 scoliosis families previously collected by our research group. The study protocol was approved by the ethics committee (approval number: 20181106003-GZ2021). After obtaining informed consent, peripheral blood samples were collected, and genomic DNA was extracted. WES was performed for sequencing, and following quality control and variant filtering, candidate SNP data were generated for subsequent intersection analysis with the SMR results.Detailed sequencing data, sampling details, and systematic sample information for the 34 families are provided in Supplementary Table 1.

cis-eQTL Data

The cis-eQTL summary statistics used in this study were obtained from the eQTLGen Consortium, which provides standardized cis-eQTL summary files suitable for SMR analysis via its official website. This dataset integrates whole-blood samples from 31,684 individuals of European ancestry across 37 independent studies, offering comprehensive and high-quality cis-eQTL data

⁸

. The data are available at:https://www.eqtlgen.org/cis-eqtls.html.

GWAS Data

The scoliosis GWAS summary statistics were obtained from the FinnGen R12 database, which integrates genomic and health data from 500,000 Finnish biobank participants. FinnGen is a large-scale genomics initiative that has analyzed over 500,000 Finnish biobank samples and linked genetic variants to health data to understand disease mechanisms and susceptibility

⁹

. Case-control summary statistics for scoliosis were retrieved using the keyword “scoliosis”. The data are available at: (1) FinnGen R12 summary statistics: https://www.finngen.fi; (2) Scoliosis data:https://storage.googleapis.com/finngen-public-data-r12/summary_stats/ release/finngen_R12_M13_SCOLIOSIS.gz.

Data Preprocessing

The raw summary statistics from the FinnGen R12 scoliosis GWAS and the eQTLGen cis-eQTL data were subjected to standardized preprocessing and format harmonization. For the GWAS data, core fields including SNP identifier, effect allele, reference allele, allele frequency, effect size, standard error, and P-value were uniformly extracted, and field naming conventions were standardized. Duplicate and blank SNP records were removed to ensure unique variant identification. Total sample size information was added, and the data were converted into the outcome data format required for SMR analysis. Subsequently, the cis-eQTL and GWAS data underwent allele orientation alignment, genomic coordinate matching, and gene symbol standardization to eliminate format discrepancies and naming inconsistencies across different data sources, ensuring complete matching of variant information between the two datasets. Finally, standardized input files directly applicable to SMR colocalization analysis were generated.

SMR Analysis

Based on the principle of Mendelian Randomization (MR), this study used SMR (Windows 1.3.1) software to perform colocalization and causal inference analysis between gene expression and scoliosis. By integrating cis-eQTL and GWAS summary data, this method overcomes the limitation of traditional GWAS that can only identify genetic associations but cannot pinpoint functional causal genes. Single nucleotide polymorphisms (SNPs) significantly associated with gene expression in the cis-eQTL data were used as instrumental variables, with gene expression level (cis-eQTL) as the exposure and scoliosis (GWAS) as the outcome. The effect estimates of the same variants from the GWAS summary data were integrated to quantitatively infer the potential causal effect of gene expression on the risk of scoliosis.

The significance threshold was set at P_SMR< 0.05. Since two different variants located in linkage disequilibrium (LD) may show significant mutual associations, relying solely on the significance of SMR results cannot directly distinguish true causal associations from false positives due to LD, i.e., it is difficult to exclude pleiotropic effects

¹⁰

. To effectively differentiate LD effects from true pleiotropic causal effects, Zhu et al.

¹⁰

proposed the HEIDI test based on instrumental variables. In this study,P_HEIDI > 0.05 was used as a filtering criterion to remove false positive associations caused by LD interference. Ultimately, significantly associated genes passing both filters were obtained, along with their association information, generating the SMR.filter.csv dataset, which serves as a core candidate set for subsequent functional enrichment and mechanistic analysis.

Family-Based Gene Intersection Analysis

To integrate the results from family-based WES and public omics data, and to improve the reliability and specificity of candidate genes, an intersection analysis was performed between the two independent gene sets. The significant genes identified by SMR analysis were intersected with those carrying key SNPs from the WES analysis of scoliosis families. Venn diagrams were plotted using R software (version 4.4.2) with the ggvenn package

¹¹

to visualize the overlapping genes. Finally, 31 common genes present in both datasets were identified as core candidates for subsequent functional annotation and mechanistic exploration of scoliosis pathogenesis.

Data Visualization Analysis

All genome-wide data visualization analyses were performed using R software (version 4.4.2). The dplyr and tidyr packages

^12,

¹³

were used for data manipulation, and the CMplot package

¹⁴

was used to generate Manhattan plots. A Manhattan plot was drawn with chromosomes on the x-axis and -log₁₀ P-value on the y-axis for the 31 common genes, providing a visual representation of the genome-wide distribution of scoliosis-associated SNPs across human chromosomes. The Manhattan plot clearly shows the chromosomal localization of significant association loci and preliminarily reveals the distribution pattern of scoliosis susceptibility regions, offering a visual basis for subsequent core gene screening and functional interpretation.

To further validate the colocalization association between gene expression and scoliosis risk, scatter plots of eQTL effect size versus GWAS effect size were generated for each of the 31 common genes using the magick and TeachingDemos packages

^15,

¹⁶

in R 4.4.2. In these plots, the x-axis represents the eQTL effect size of the variant on gene expression, and the y-axis represents the GWAS effect size of the variant on scoliosis risk, with the top cis-eQTL and common cis-eQTL sites annotated for each gene. These scatter plots effectively illustrate the correlation trend between eQTL effects and GWAS effects, assist in evaluating the strength of the causal association between gene expression and disease phenotype, and provide visual support for the colocalization validation of core candidate genes.

Functional Enrichment

The Gene Ontology (GO) database provides a standardized annotation system for gene functions and serves as one of the core bioinformatics tools for interpreting gene expression data. GO enrichment analysis systematically annotates gene functions from three dimensions: Biological Process (BP), Molecular Function (MF), and Cellular Component (CC), effectively identifying statistically enriched functional categories and biological characteristics under disease conditions

¹⁷

. In this study, GO enrichment analysis was performed using the clusterProfiler package

¹⁸

in R software (version 4.4.2), and bubble plots were generated. A P-value < 0.05 was considered statistically significant for enrichment results

¹⁹

.

The Kyoto Encyclopedia of Genes and Genomes (KEGG) is a comprehensive molecular function database constructed by the Bioinformatics Center of Kyoto University, covering 17 sub-databases and integrating information on cellular processes, genetic information, human diseases, and more

^20,

²¹

. Enrichment analysis based on KEGG data is a common approach to identify biological pathways and functional modules involved in target genes. In this study, KEGG enrichment analysis was performed using the clusterProfiler package in R 4.4.2, and bubble plots were visualized. A P-value < 0.05 was considered statistically significant for KEGG enrichment results

¹⁹

.

Family-based analysis of scoliosis cohorts

Based on the integrated analysis of WES data from 34 scoliosis families and SMR results, a total of 31 common genes present in both datasets were identified. These genes passed stringent filtering criteria (P_SMR < 0.05 and P_HEIDI > 0.05), indicating that they carry key genetic variants at the family level while also showing causal associations with scoliosis risk at the population level. The P_SMR and P_HEIDI values for these 31 common genes are detailed in Table 1.

Colocalization validation at the LINC00299 locus

The Manhattan plot illustrates the distribution of the 31 common genes across all chromosomes (Figure 1A), revealing that multiple genomic loci are significantly associated with scoliosis susceptibility. Colocalization analysis indicated that the long non-coding RNA LINC00299 exhibits a significant cis-eQTL effect in the context of scoliosis genetic susceptibility.

The single‑gene scatter plot for LINC00299 (Figure 1B) showed that most data points fall along a negative correlation trend line, indicating that alleles increasing its expression tend to decrease the risk of scoliosis. That is, higher expression of LINC00299 is associated with lower disease susceptibility, suggesting that down-regulation of this gene may promote abnormal spinal development. The top cis-eQTL for this locus (marked by a red triangle) also perfectly aligns with the overall fitted trend. Statistical analysis showed a significant SMR causal association test (P_SMR = 0.01113481, P < 0.05) and a HEIDI heterogeneity test (P_HEIDI = 0.06945593, P > 0.05) indicating no horizontal pleiotropy, thus excluding false colocalization bias.

Further analysis using the regional colocalization plot (Figure 1C) revealed that within the 8.5-9.0 Mb region on chromosome 2, the significant association peak from the scoliosis GWAS and the strong cis-eQTL association peak for LINC00299 are highly overlapping. The two genetic signals are consistently distributed across the genomic interval, exhibiting clear colocalization characteristics.

eQTL-GWAS effect characteristics of core candidate genes

Through integrative colocalization analysis of scoliosis GWAS signals and cis-eQTL signals, this study identified four candidate genes with potential pathogenic relevance:NCAM1, AGPAT3, SMIM12, and CTSH.

For NCAM1(Figure 2A), AGPAT3(Figure 2B), and SMIM12(Figure 2C), the cis-eQTL effect sizes showed a stable positive correlation with the GWAS disease effect sizes, indicating that the direction of genetic variant effects on gene expression was consistent with the direction of variant effects on scoliosis risk. In contrast, CTSH(Figure 2D) exhibited a clear negative correlation between eQTL effect sizes and GWAS effect sizes. Moreover, the top cis-eQTL (marked by red triangles) for each gene fell exactly on the overall fitted trend line, suggesting that altered expression levels of these four genes provide reliable genetic evidence for an association with the risk of scoliosis.

Functional annotation of key genes associated with scoliosis

Functional enrichment analysis was conducted on the 31 consensus genes associated with scoliosis. GO analysis identified significant enrichment across three standard functional dimensions(Figure 3A-B).

In the BP category, these genes were mainly clustered in pathways related to energy coenzyme biosynthesis and neural signal transduction. The most significantly enriched terms included NAD biosynthetic process, nicotinamide nucleotide biosynthetic process, pyridine nucleotide biosynthetic process, pyridine-containing compound biosynthetic process, Rho protein signal transduction, and regulation of cold-induced thermogenesis.

For the CC category, prominent enrichment was found in synaptic complexes, neuronal structures and specialized cell membrane domains, such as basal part of cell, postsynaptic density, asymmetric synapse, juxtaparanode region of axon, postsynaptic specialization, and neuron-to-neuron synapse.

In the MF group, enriched functions were primarily involved in catalytic activity, ion channel regulation and molecular binding. Key significant terms covered exopeptidase activity, dipeptidyl-peptidase activity, bile acid binding, myristoyltransferase activity, nuclear export signal receptor activity, and calcium-activated potassium channel activity. Full detailed statistics for all GO enrichment results are listed in Table 2.

Further KEGG pathway enrichment analysis revealed that these genes were significantly enriched in nicotinate and nicotinamide metabolism, insulin secretion, glucagon signaling pathway, Hypoxia-Inducible Factor-1 (HIF-1) signaling pathway, cofactor biosynthesis, and the renin-angiotensin system (Figure 3C).