Association Between Sleep Apnea Syndrome and Osteoarthritis: Insights from Bidirectional Mendelian Randomization and Bioinformatics Analysis

Introduction

Sleep apnea syndrome (SAS) is characterized by the partial or complete collapse of the throat during sleep, leading to apnea and hypoventilation.¹ This intermittent hypoventilation causes the body to experience intermittent hypoxia, exposing all body cells to a cycle of deoxygenation-reoxygenation. This process leads to oxidative stress and systemic inflammation, thereby damaging other organs.² If this type of injury occurs in the articular cartilage, it may lead to the development of osteoarthritis(OA).

SAS prevalence is notably elevated in individuals with OA. Shannon Stark Taylor and others have noted a 66% incidence rate of SAS among OA patients.³ The two conditions share many similarities, including aging and metabolic disorders as common risk factors, as well as cardiovascular diseases and neurodegenerative diseases being common comorbidities.⁴ However, there’s no consensus on the mutual relationship between SAS and OA. Kanbay Asiye et al found a strong correlation between SAS and OA.⁵ Conversely, Diaz et al found in a retrospective cohort study that OA is not related to chronic SAS.⁶ This highlights that observational studies are easily influenced by confounding factors, making it difficult to establish a precise causal relationship between SAS and OA.

To further elucidate the causal relationship between them, our study employed a bidirectional two-sample Mendelian Randomization (MR) approach. MR analysis, leveraging the principles of Mendel’s Second Law, effectively manages confounding bias and reverses causality dilemmas in observational epidemiology, thus presenting a novel method for ascertaining causality. In this study, exposure-related single nucleotide polymorphisms (SNPs) were chosen as instrumental variables (IVs), thereby utilizing genetic variations to deduce the causal link between exposure and outcomes.^7–9

The results of the Mendelian randomization analysis suggest a potentially bidirectional causal relationship between SAS and OA. To further explore the connection between the two, this study employs bioinformatics to investigate the mechanisms underlying this relationship. Bioinformatics combines multiple disciplines, such as statistics and computer science, to manage, analyze, and interpret vast and complex biological data. With the continuous development of high-throughput sequencing and chip-based sequencing technologies, it has become increasingly easy to obtain transcriptome data on SAS and OA from public databases (such as the GEO database) and subsequently use bioinformatics to analyze the data to find the comorbidity mechanisms of the two. This approach aims to provide some guidance for the prevention and treatment of both diseases.

Material

Data Sources

Data for the Mendelian randomization study were sourced from extensive genome-wide association studies (GWAS). The GWAS data for SAS, OA and its subtypes (hip and knee OA), and colorectal cancer were obtained from the Finnish database (https://www.finngen.fi/en/access_results). The dataset for SAS includes 375,657 samples from the European population, with 38,998 in the disease group and 336,659 in the control group. OA samples, also from Europeans, included 22,254 in the disease group and 240,862 in the control group. The knee OA group included 44,688 European samples, with 240,862 in the control group. The hip OA group included 22,254 European samples, with 240,862 in the control group. Colorectal cancer samples, also from Europeans, included 6509 cases and 287,137 controls. The GWAS data for immune cells are from the GWAS Catalog (https://www.ebi.ac.uk/gwas/downloads/summary-statistics), with IDs from GCST0001391 to GCST0002121,¹⁰ including a total of 731 immune phenotypes. These comprised Absolute Cell Count (AC) (118 instances), Median Fluorescence Intensity for surface antigen levels (MFI) (389 instances), Morphological Parameters (MP) (32 instances), and Relative Cell Count (RC) (192 instances). MFI, AC, and RC features included combinations of B cells, CDCs, mature T cells, monocytes, myeloid cells, TBNK (T cells, B cells, natural killer cells), and Treg, while MP features included combinations of CDC and TBNK.

Transcriptome data for OA and SAS were obtained from the GEO database (https://www.ncbi.nlm.nih.gov/geo/). Datasets GSE12021, GSE55235, and GSE55457 all originated from platform GPL96 and included samples of OA and normal samples, OA is diagnosed according to the criteria for OA;¹¹ GSE32317 was from platform GPL570 and did not include normal samples but included early-stage OA samples (cartilage degeneration without the presence of full-thickness cartilage loss, and a Kellgren-Lawrence score of ≤2) and late-stage OA samples (full-thickness cartilage loss, with a Kellgren-Lawrence score of >2). GSE135917 originates from platform GPL6244 and includes both normal samples and samples of sleep apnea syndrome. GSE75097 (GPL10904) focuses exclusively on sleep apnea syndrome, presenting samples at different stages (based on the Apnea-Hypopnea Index (AHI)): primary snoring (PS, AHI<5), moderate to severe sleep apnea (MSO, 15<AHI≦50), very severe sleep apnea (VSO, AHI>50), along with very severe sleep apnea patients undergoing long-term continuous positive airway pressure therapy (CPAP) (VSOC, regular CPAP use: >4 hours/night, >1 year), this project has selected the first three categories of GSE75097 for analysis. The datasets are freely available to the public, and each GWAS and GEO dataset has received ethical approval from its respective institution. Table S1 shows detailed information about the GWAS and GEO datasets analyzed in this study.

Selection of Instrumental Variables

During the initial screening of instrumental variables, a threshold of p<5×10−8 was set to identify statistically significant single nucleotide polymorphisms (SNPs). However, when selecting instrumental variables for immune cells, the number of SNPs identified at a threshold of p<5×10−8 was limited. To obtain more instrumental variables, the threshold was set to p<1×10−5. This study utilized local Plink software to test for linkage disequilibrium between SNPs, setting the allele imbalance coefficient r² to 0.001 and the region width to 10,000kb to minimize the impact of genetic pleiotropy on the results. The explanatory variance R² of the genetic instruments used in the study was calculated using Equation (1) ¹²:

(1)

where EAF is the allele frequency and beta is the effect size of a specific SNP. Next, the F-statistic for each SNP was calculated using Equation (2) to assess its strength¹³:

(2)

where R² is the variance in the exposure factor explained by a specific SNP, N is the sample size of the exposure data in the GWAS study, and an F-value greater than 10 is sufficient to eliminate bias.

Mendelian Randomization Analysis and Sensitivity Test

Inverse variance weighting (IVW)¹⁴ is used as a key MR analysis method to estimate the causal effect of SAS on OA and its subtypes. Its basic principle can be simply understood as follows: first, a weight is assigned to the effect of each exposure variant on the outcome variation, which is inversely proportional to the variance of the study results, meaning that results with smaller variance receive greater weight. Therefore, inverse variance is used to represent the reciprocal of variance, with the larger the inverse variance, the greater the weight. Then, in estimating the overall effect, a weighted average based on the weight of each result is used to enhance the influence of more reliable (less error-prone) results, thus making the estimates more reflective of actual conditions.

Additionally, this study employs four other methods as supplements to IVW. The MR Egger method treats the intercept term as part of the analysis, allowing it to vary freely during the estimation process. By estimating this non-zero intercept, the MR-Egger method can identify and correct for bias, thus providing a more robust estimate of causal effects;¹⁴

The Weighted Median method involves ranking each SNP’s effect estimate and calculating the weighted median as the estimate of causal effect. This method can provide a robust estimate of the causal effect, even if up to 50% of the instrumental variables are invalid.¹⁵

The Simple mode method selects one or a few highly correlated SNPs as instrumental variables, using the effect values of these instruments to estimate the causal effect of exposure on outcomes. The Weighted mode method is based on providing a weighted average estimate using valid instrumental variables, with each instrument weighted according to its predictive power for exposure or other relevance. This weighting strategy can improve the accuracy of the overall estimate.

This study uses the Egger intercept for pleiotropy testing to ensure the effect of horizontal pleiotropy. If the intercept shows statistical significance, the MR-PRESSO (Mendelian Randomization-Pleiotropy Residual Sum and Outlier) method is used to remove outlier SNPs. The MR-PRESSO method involves three steps:¹⁶ i) detecting horizontal pleiotropy; ii) correcting for pleiotropy by excluding genetic variants exhibiting horizontal pleiotropy; iii) comparing differences in causal association before and after correction, and re-testing for pleiotropy. Additionally, this study calculates Cochran’s Q value to quantify the heterogeneity of individual causal effects. If the p-value of Cochran’s Q is less than 0.05, heterogeneity is considered to exist, and a random-effects IVW model is adopted; otherwise, a fixed-effect IVW model is used. All statistical analyses were performed using the “TwoSampleMR” and “MR-PRESSO” packages in R software.

Expression Difference Analysis

Datasets GSE12021, GSE55235, and GSE55457, originating from platform GPL96 and comprising synovial tissue samples, were combined, followed by batch correction and normalization of the merged data. After preparing the data for each disease, the “Limma” package in R software was used to compare the SAS and OA datasets, obtaining differentially expressed genes (DEGs) between disease and control groups with a selection threshold set at p-value <0.05 and |log2FC| >0.5. Volcano plots were utilized to display the differential analysis results for each group, where green represents downregulated differential genes and red represents upregulated differential genes.

Weighted Gene Co-Expression Network Analysis (WGCNA)

WGCNA is utilized to discover relationships among genes and between genes and clinical traits. Genes highly correlated in expression are often involved in specific biological processes or pathways and thus are grouped into a module. Additionally, the WGCNA method can associate gene sets with clinical traits, thus identifying modules related to clinical features, thereby helping to unveil potential biological mechanisms. When performing WGCNA, standardized mRNA expression data (using the “Limma” package) are used as input files to enable the exploration of correlations between genes, and between gene modules and clinical traits (such as SAS and OA).

To conduct WGCNA effectively, (1) perform hierarchical clustering analysis using the R package “gplots” to identify potential outlier samples; (2) then use the “pickSoftThreshold” function to choose an appropriate soft-thresholding “power” parameter, typically ranging from 1 to 20; (3) subsequently, convert the correlation matrix into an adjacency matrix, and further into a topological overlap matrix (TOM), to reflect the topological similarity among genes; (4) following that, construct a hierarchical clustering tree of genes based on average linkage, and identify different gene modules using the dynamic tree cutting algorithm (minModuleSize=30). Finally, merge similar modules within each group by adjusting the cut height; (5) Lastly, conduct correlation analysis between gene modules and clinical phenotypes (such as SAS and OA) using the Pearson correlation coefficient.

Identification of Shared Genes and Functional Enrichment Analysis

By integrating DEGs and module genes identified by WGCNA, shared genes for SAS and OA are obtained. This integration facilitates the identification of common molecular mechanisms between these diseases. The process then uses the Venn package to draw Venn diagrams for visualization, providing a clear depiction of the overlap between different gene sets. The “clusterProfiler” package is used for enrichment analysis of shared genes in Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathways, enabling the identification of where significant enrichment in biological functions and signaling pathways occurs. Finally, bar graphs show the results of significant enrichment when the p-value is less than 0.05, thus offering a quantitative understanding of the enrichment levels.

Colocalization Analysis

Colocalization analysis is used to identify genetic variant loci associated with related traits or phenotypes and to determine whether there are shared causal variants in a given genomic region. This study conducted a colocalization analysis between SAS GWAS and OA GWAS. Colocalization analysis involves five hypotheses: H0: There is no significant association between phenotype 1 and phenotype 2 across all SNP loci in a genomic region; H1/H2: Either phenotype 1 or phenotype 2 shows significant association across all SNP loci in a genomic region; H3: There is a significant association between phenotype 1 and phenotype 2 across all SNP loci in a genomic region, but it is driven by different causal variant loci; H4: There is a significant association between phenotype 1 and phenotype 2 across all SNP loci in a genomic region, and it is driven by the same causal variant locus. These hypotheses are crucial as they fundamentally test the posterior probabilities of these five hypotheses. The posterior probability values of hypothesis H4 range from 0% to 100%, and typically, loci with an H4 posterior probability greater than 0.75 are considered colocalized loci. To operationalize this analysis, this study employs the “coloc” package for colocalization analysis, which is utilized to calculate SNPs within approximately 1MB of the shared genetic region, and subsequently, visualizes the results using the “locuscomparer” package.

Immune Microenvironment Immune Score

The “estimate” package was used to analyze the status of immune cells in OA and SAS samples, assessing the proportion of immune cells. By conducting this analysis, the “estimate” package can be utilized to acquire the immune score. The immune score reflects the number of immune cells in each sample, a higher score denotes greater quantities of immune cells. Furthermore, the “ggpubr” package is used to draw violin plots to visualize the score differences between the disease group and the normal group, as well as among different stages of the disease group. A p-value less than 0.05 is considered to indicate a statistically significant difference between the two groups.

Results

Bidirectional Mendelian Randomization of Two Samples of Sleep Apnea and OA

We employ MR analysis to examine the association between SAS and OA, as well as OA subtypes (hip OA, knee OA). Subsequently, we reversed the direction of the MR analysis. Through a rigorous selection process, we identified SNPs strongly associated with both SAS and OA, with each instrumental variable having a P-value< 5×10-8.

The results primarily relied on the IVW model, while other methods were used for supplementary analysis. Firstly, the MR-Egger intercept test did not find significant horizontal pleiotropy (Table S2). Then combined with the heterogeneous analysis results (Table S3). If heterogeneity existed, the random-effects IVW model was adopted; Otherwise, the fixed-effects IVW model was employed.

The causal effects of SAS on OA were reported in the form of odds ratios (OR) and forest plots in the MR analysis (Figure 1). SAS was a risk factor for OA (IVW OR = 1.33, 95% CI = 1.18–1.50, p = 0.000004), knee OA (IVW OR = 1.48, 95% CI = 1.28–1.72, p = 0.0000001), and hip OA (IVW OR = 1.35, 95% CI = 1.13–1.60, p = 0.001) (Table S4). Consistent results were also observed using MR-Egger and Weighted Median methods.

Figure 1 Results of Mendelian randomization of OA and SAS.

Subsequently, we performed a reverse analysis by swapping the exposure and outcome variables. OA (IVW OR = 1.22, 95% CI = 1.10–1.36, p = 0.000195) and knee OA (IVW OR = 1.15, 95% CI = 1.05–1.24, p = 0.001) were identified as significant risk factors for SAS (Figure 1, Table S4). MR-Egger and Weighted Median analyses also yielded consistent results. However, no significant association was found between SAS and hip OA (IVW OR = 0.99, 95% CI = 0.91–1.08, p = 0.892) (Figure 1, Table S4). The results of the five methods were not statistically significant (p>0.05).

Identification of Differentially Expressed Genes and Screening of Key Modules by WGCNA

Firstly, we performed differential expression analysis on the OA combined dataset. A total of 2170 differentially expressed genes were identified using the filtering criteria of LogFC > 0.5 and P < 0.05. Among them, 1081 genes were downregulated and 1089 genes were upregulated. The results were visually represented in a volcano plot (Figure 2A). Similarly, differential expression analysis was also conducted on the SAS dataset using the same filtering criteria. A total of 2069 genes met the criteria, with 264 genes upregulated and 1805 genes downregulated. The results were visualized in a volcano plot as well (Figure 2B).

Figure 2 Volcano plot and WGCNA for SAS and OA: (A) Volcano plot of OA (B) Volcano plot of SAS (C) Determination of soft threshold power for OA (D) Determination of soft threshold power for SAS (E) Cluster tree of highly connected genes in key modules of OA.(F) Cluster tree of highly connected genes in key modules of SAS (G) The relationship between modules and traits of OA, Each cell contains a correlation and a p-value (H) Modular-feature relationship in SAS. Each cell contains a correlation and a p-value.

In addition, we also performed WGCNA. Using the soft thresholding approach, a co-expression network was constructed in this study. The parameter “power” is crucial for maintaining scale-free topology in a co-expression network, as gene expression data is likely to exhibit scale-free characteristics. Therefore, a fitting index greater than 0.85 was considered indicative of scale-free topology in the OA group, and a power value of 8 was set accordingly (Figure 2C). Similarly, WGCNA was also applied to the SAS group, with a power value of 13 identified as the optimal soft threshold (Figure 2D). An adjacency matrix was generated using the adjacency function, and a hierarchical clustering was performed based on the TOM dissimilarity measurement, as shown in Figure 2E and F. In OA group, a total of 8 co-expression modules were identified, with modules showing P < 0.05 considered as key modules. As shown in Figure 2G, the MEbrown and MEturquoise modules exhibited significant correlation and contained 1185 genes. In SAS group, a total of 13 modules were identified, with the MEblue, MEgreenyellow, MEtan, and MEturquoise modules showing significant correlation with the disease (Figure 2H), including 1636 genes.

Enrichment Analysis of Shared Genes and Pathway Mendelian Randomization

We intersected the genes screened by DEGs and WGCNA to obtain shared genes between the two diseases, and performed enrichment analysis on the shared genes to explore the functional pathways of the two diseases. Figure 3A illustrates the overlap between the DEGs and WGCNA genes related to SAS and OA, consisting of a total of 78 genes. We speculate that these 78 genes may be associated with the pathogenesis of both SAS and OA, indicating a shared relationship between the two conditions. Through functional annotation and enrichment analysis of these genes (Figure 3B and C), we aimed to investigate potential biological changes between the two diseases.

Figure 3 Analysis of shared gene enrichment and pathway MR in SAS and OA: (A) Venn diagram of SAS, OA and shared genes (B) GO enrichment analysis of shared genes. (C) KEGG enrichment analysis of shared genes. (D) Results of MR of pathways in colorectal cancer and SAS and OA.

Gene Ontology (GO) analysis revealed that these shared genes are significantly enriched in biological processes such as leukocyte migration, leukocyte chemotaxis, response to corticosteroids, and response to glucocorticoids. This suggests that the migration and chemotaxis of immune cells may play an important role in the common mechanisms of SAS and OA.

Additionally, KEGG enrichment analysis showed significant involvement of TNF signaling pathway, IL-17 signaling pathway, colorectal cancer, lipid and atherosclerosis. We noticed that the association of inflammatory factors and lipid metabolism with both diseases is in line with previous research. However, the relationship between colorectal cancer and both SAS and OA has not been reported before. Therefore, we conducted Mendelian randomization analysis on the pathway of colorectal cancer in relation to SAS and OA (Figure 3D, Table S5). The colorectal cancer pathway was found to be a protective factor for SAS (IVW OR = 0.96, 95% CI = 0.93–0.99, p = 0.025). However, for OA (IVW OR = 1.04, 95% CI = 1.00–1.08, p = 0.036) and hip OA (IVW OR = 1.07, 95% CI = 1.01–1.14, p = 0.028), it was identified as a risk factor. There was no evidence supporting the relationship between this pathway and knee OA (IVW OR = 1.04, 95% CI = 0.99–1.10, p = 0.13)(Table S5). The MR-Egger intercept test did not detect significant levels of pleiotropy (Table S6). The results of heterogeneity tests are shown in Table S7

Colocalization Analysis of Shared Genes

To further screen the key genes of the two diseases from 78 shared genes, we conduct colocalization analysis. The results indicated that SAS and OA may share a causative variant within the JUNB (PP.H4=0.812), COL8A1 (PP.H4=0.973), FOSB (PP.H4=0.935), and IER2 (PP.H4=0.811) loci (Figure 4A–D, Table S8). Based on the co-localization analysis, we identified four potential genes that could potentially play a role in the relationship between SAS and OA, with evidence suggesting a common genetic effect between SAS and OA.

Figure 4 Results of JUNB(A), COL8A1(B), FOSB(C), and IER2(D) colocalization analysis.

Immune Microenvironment Immune Score

Enrichment analysis suggests that SAS and OA are associated with immunity and immune cells, which we further confirm with immune microenvironment scoring. We utilized an immunomicroenvironment analysis method to score the immune responses in OA and normal individuals, as well as different stages of OA and SAS. Statistically significant differences in immune scores were observed between the normal and OA groups, as well as between early-stage and moderate-severe stage SAS, and between early-stage and severe stage SAS. However, there were no significant differences in immune scores between early and late stages of OA, or between SAS and normal individuals. (Figure 5).

Figure 5 Violin plot of immune score (A) OA group and control group (B) early OA group and end OA group (C) SAS group and control group (D) violin plot of SAS at different stages (“**”:P<0.01, “*”: P<0.05, “ns”: P>0.05, PS: primary snoring, MSO: moderate to severe obstructive sleep apnea, VSO: very severe obstructive sleep apnea).

Mendelian Randomization Analysis of Immune Cells and Sleep Apnea and Osteoarthritis

Analysis of the immune microenvironment suggests that the content of immune cells differs between the OA group and the control group, as well as across different stages of SAS, indicating that immune cells may be related to both diseases. To further screen immune cell phenotypes, we conducted Mendelian randomization analysis of 731 immune cell phenotypes with both. Using 731 immune cell types as exposure factors and SAS and OA as outcomes, a Mendelian randomization analysis was conducted with IVW as the primary analytical method. In the presence of heterogeneity, a random-effects IVW model was used; otherwise, a fixed-effect IVW model was applied. In total, 57 immune phenotypes were related to SAS, of which 31 immune phenotypes were protective factors for SAS, 26 immune phenotypes were risk factors for SAS (Table S9, Figure 6A and B), and the MR Egger intercept indicated no pleiotropy (Table S10). There were 95 immune phenotypes associated with OA, of which 43 served as protective factors for OA, 52 immune phenotypes were detrimental to OA (Table S11, Figure 6A and B), with the MR-Egger intercept indicating no evidence of horizontal pleiotropy (Table S12). There are 6 immune phenotypes influencing both OA and SAS, of which 2 are protective factors and 4 are risk factors (Figure 6A–C).

Figure 6 Results of Mendelian randomization of immune cell phenotypes and SAS and OA. (A) Venn diagram of immunophenotypes that protect against OA and immunophenotypes that protect against SAS. (B)Venn diagram of immunophenotypes that promote OA and immunophenotypes that promote SAS. (C)Forest plots of immunophenotypes that can affect both OA and SAS.

CD4 on CD39+ resting Treg cells offer protection against OA (IVW OR=0.946, 95% CI: 0.913–0.980, p=0.0024) and SAS (IVW OR=0.937, 95% CI: 0.901–0.975, p=0.0012). CD4 on HLA DR+ CD4+ cells also confer protection against OA (IVW OR=0.97, 95% CI: 0.941–0.999, p=0.041) and SAS (IVW OR=0.972, 95% CI: 0.946–0.9985, p=0.039).(Table S13)

However, IgD on IgD+ CD24+ has an adverse effect on both OA (IVW OR=1.024, 95% CI:1.0015–1.048, p=0.0366) and SAS (IVW OR=1.024, 95% CI:1.004–1.043, p=0.0158). CD33 on CD33br HLA DR+ is a common risk factor for both OA (IVW OR=1.011, 95% CI:1.002–1.01996, p=0.014) and SAS (IVW OR=1.013, 95% CI:1.0022–1.0235, p=0.018). CD33 on CD33br HLA DR+ CD14- is also detrimental to both OA (IVW OR=1.012, 95% CI:1.003–1.021, p=0.01003) and SAS (IVW OR=1.012, 95% CI:1.001–1.023, p=0.026). CD64 on CD14- CD16- is a common risk factor for OA (IVW OR=1.048, 95% CI:1.0199–1.078, p=0.0008) and SAS (IVW OR=1.033, 95% CI:1.0003–1.067, p=0.048).(Table S13)

Discussion

The MR method uses SNPs strongly associated with SAS and OA as instrumental variables, excluding potential confounding factors that might affect experimental outcomes, hence the conclusions are more reliable than observational studies. Additionally, by employing bioinformatics methods to deeply investigate the mechanisms of both diseases, it offers strategies for future treatment and prevention of these diseases.

SAS is a causative factor in OA, especially in knee OA. Silva Andressa et al found a strong correlation between the severity of SAS and the grading of OA.¹⁷

SAS can lead to intermittent hypoxia and inflammation, thereby inducing cellular damage.¹⁸ Patients experience oxidative stress due to intermittent hypoxia during sleep, leading to mitochondrial DNA damage and dysfunction.^19,20 Intermittent hypoxia leads to the production of reactive oxygen species, TNF, inflammatory cytokines (IL-2, IL-4, IL-6), lipid peroxidation, and extracellular DNA damage, further promoting the generation and progression of inflammation.²¹ Inflammation leads to the degradation of components of the extracellular matrix of cartilage cells, metabolic disorder of chondrocytes, and promotes oxidative stress and cartilage damage.²² The inflammatory environment can disrupt the synthetic metabolic balance of cartilage cells, promote the production of tissue degrading enzymes (especially MMP-13 and ADAMTS-5), and mediate the degradation of the cartilage matrix.²³ Resulting in damage to the structure of joints and the emergence of joint pain.

Body weight and metabolic disorder may play a role in both diseases. Obesity is a fundamental factor in the occurrence of SAS.²⁴ Obesity and metabolic abnormalities can damage joint tissues.²⁵ Adipose tissue releases fatty factors, increasing pro-inflammatory mediators, a phenomenon known as “metabolic inflammation”, which is a mechanism by which obesity damages joint tissues.²⁶

SAS is a detrimental factor for OA, possibly due to factors such as hypoxia, inflammation, and metabolic disorder caused by SAS. Further research on these mechanisms will help to better prevent and treat SAS and OA.

This study also found that knee OA is a risk factor for SAS, a finding not yet reported in the literature, with most studies typically considering SAS as the cause and OA as the outcome. Knee OA can exacerbate obesity, with Liao, J and others finding a higher incidence of obesity among patients with bilateral knee OA than among those with unilateral knee OA.²⁷ Obesity is a risk factor for SAS, with literature reporting that an increase in Body Mass Index (BMI) and neck circumference raises the incidence of SAS.²⁸ Thus, knee OA may lead to SAS by promoting obesity. However, there is insufficient evidence on the impact of hip arthritis on SAS, and large-scale clinical randomized controlled trials are needed in the future to study the relationship between hip arthritis and SAS.

Shared genes between the two diseases are significantly enriched in the TNF pathway, IL-17 pathway, lipid metabolism pathway, and colorectal cancer pathway. Fiedorczuk et al found that the plasma TNF concentrations in patients with SAS positively correlate with the severity of the disease.²⁹ Meta-analysis shows that continuous positive airway pressure (CPAP) treatment for SAS can reduce TNF-α levels in patients.³⁰ TNF-α is implicated in the early development of OA, with research showing elevated expression of TNF-α in the joint cartilage of mild OA cases.³¹ Stannus et al found that serum TNF-α levels in the elderly are associated with the loss of knee joint cartilage.³² Hussein M R et al reported that levels of IL-17 and TNF-α in the serum and synovial fluid of arthritis patients are higher than those in the healthy control group.³³

IL-17 also plays a significant role in both diseases. Plasma IL-17 levels are elevated in patients with SAS.³⁴ IL-17 is a cytokine associated with Th17, potentially regulating the function of helper T cells.³⁵ Th17 is involved in the occurrence and development of SAS, with studies finding an imbalance in the Th17/Treg cell ratio in patients with SAS.³⁶ This indicates the significant role of IL-17 and Th17 in SAS. IL-17 can also cause OA by stimulating human OA synovial fibroblasts to increase production of VCAM-1, thereby enhancing monocyte adhesion and leading to the onset of OA.³⁷

There is a close relationship between lipid metabolism and SAS. Research indicates that all non-traditional lipid indicators positively correlate with the severity of SAS.³⁸ Chronic hypoxia in SAS leads to the production of a large amount of reactive oxygen species (ROS), activating the JNK/SREBP/ACC pathway in neurons and glial cells, ACC catalyzes the carboxylation of acetyl-CoA to malonyl-CoA, thereby leading to an increase in fatty acid synthesis and abnormal accumulation.³⁹ In OA, changes in lipid metabolism may be a cause of the disease, studies have observed lipid deposition in joints before histological changes occur, proteomic analysis has revealed significant connections between OA and lipid metabolism, with adipokines identified as key regulatory factors in the pathogenesis of OA.⁴⁰

In-depth future research into TNF, IL-17, and lipid metabolism pathways could offer a valuable theoretical foundation for treating the comorbidity of OA and SAS.

In GO enrichment analysis, shared genes are significantly enriched in biological processes related to leukocyte migration, chemotaxis, and immune response. The possible association of both conditions with immune cells is confirmed by the analysis of the immune microenvironment. Intermittent hypoxia due to SAS leads to significant deposition and adhesion of leukocytes on endothelial cells, resulting in vascular damage.⁴¹ Damage to blood vessels can aggravate injuries to organs, including the respiratory tract. Research has identified a correlation between imbalances in T lymphocyte immunity and SAS.⁴² T cells contribute to airway inflammation and increased reactivity by recruiting neutrophils and inducing contractions in the airway smooth muscles.⁴³ Airway constriction leading to impaired ventilation is one of the pathogenic mechanisms of SAS. Immune cells also play an important role in OA. These cells can secrete various cytokines, among which Th17 cells produce IL-17, which in turn activates other immune cells (such as fibroblasts and macrophages), promoting the release of IL-6, TNF-α, and other chemokines.⁴⁴ These inflammatory factors can further amplify the inflammation process of OA, ultimately exacerbating the condition.⁴⁵ This study conducted Mendelian randomization analysis using 731 immune cell subtypes as exposures, with OA and SAS as outcomes, finding that 57 immune cell phenotypes are associated with SAS, and 95 with OA. In the future, immune cells may represent a direction for the treatment of both diseases.

SAS and OA share genetic variations at the JUNB, COL8A1, FOSB, and IER2 gene loci, and these four genes may play a significant role in the relationship between SAS and OA. Previous studies have shown that JUNB is highly expressed in Th17 cells and regulates the expression of Th17 hallmark genes.⁴⁶ JUNB is also crucial for the development of Th17 cells.⁴⁶ Th17 cells play a significant role in both OA and SAS. Lin et al proposed a mechanism involving the JUNB/FBXO21/ERK axis, which regulates autophagy and promotes the progression of cartilage degeneration in OA.⁴⁷ COL8A1 may be involved in stabilizing the extracellular matrix of chondrocytes, and silencing COL8A1 inhibits the expression of the COL2A1 gene, exacerbating OA.⁴⁸ COL8A1 may also be involved in the onset and progression of SAS by affecting smooth muscle function. As part of the extracellular matrix, COL8A1 can facilitate the proliferation of smooth muscle cells.⁴⁹ Smooth muscle plays a crucial role in maintaining airway patency, and when smooth muscle is damaged, it leads to airway constriction.

This study has several limitations. Firstly, the lack of individual information in GWAS data means differences among individuals could lead to biased results. Secondly, our choice of GWAS data from European populations may introduce certain errors if the research results are generalized to other populations. Thirdly, the limited number of GEO datasets for SAS could result in bias due to limited sample sizes. Finally, this study lacks fundamental experiments in exploring the impact of immune cells on SAS and OA. Future studies may need to conduct specific experiments to explore the effects of the aforementioned immune cell subtypes on SAS and OA.

This research establishes SAS as a novel risk factor for OA and identifies knee OA as a risk factor for SAS. Patients suffering from both OA and SAS require prompt SAS intervention, which holds significant importance for treating OA. Additionally, JUNB, COL8A1, FOSB, and IER2 might serve as future research directions for understanding the comorbidity between these two diseases. Finally, this study highlights the important roles of 6 immune cell phenotypes in these two diseases, offering new insights for the development of related therapeutic approaches.

Conclusion

A bidirectional causal link exists between SAS and OA, necessitating the proactive management of SAS simultaneously with OA treatment. Notably, the 4 genes (JUNB, COL8A1, FOSB, IER2) and 6 immune phenotypes are crucial for both diseases, these provide hopeful targets for future interventions against these two diseases.