Genetic and Molecular Interconnections Between Chronic Obstructive Pulmonary Disease and Osteoporosis: Insights from Single-Cell and Mendelian Randomization Analyses

Introduction

Chronic obstructive pulmonary disease (COPD) is a widespread chronic lung disorder, clinically defined by long-term and typically worsening airflow obstruction. This condition is marked by an ongoing inflammatory reaction within the airways and lung tissue, primarily triggered by the inhalation of harmful particles or gases.¹ According to the Global Burden of Disease Study 2021, COPD affected approximately 213.4 million people worldwide in 2021 and was responsible for 3.7 million deaths, ranking as the third leading cause of death globally.² The burden is particularly severe in China, where the prevalence of COPD is 13.7% among adults aged 40 years and older.³ Osteoporosis (OP) is a systemic skeletal disorder characterized by a reduction in bone mineral density, compromised bone microarchitecture, and the degradation of bone matrix quality.^4,5 Osteoporosis affects approximately 500 million people worldwide. The global prevalence of osteoporosis is 18.3% in the general population and 21.7% among the elderly.⁶ In 2021, osteopenia and osteoporosis led to approximately 459,661 deaths and 17.3 million disability-adjusted life years (DALYs) globally.⁷ In recent years, the prevalence of OP has risen markedly, largely due to the accelerated aging of the population, thereby placing a substantial burden on global healthcare systems.⁸

Numerous studies have confirmed the high prevalence of OP among individuals with COPD.^9,10 A 2019 meta‑analysis reported a pooled global OP prevalence of 38% in COPD patients, with COPD increasing the risk of OP nearly threefold (OR = 2.83).¹¹ The pathogenesis of COPD‑associated OP is multifactorial, involving shared risk factors (eg., aging, smoking, low BMI, physical inactivity),¹² as well as COPD‑specific drivers such as systemic inflammation, oxidative stress, hypoxemia, and long‑term corticosteroid use. Increasing evidence also suggests cross‑talk between muscle and bone via “myokines” and “osteokines”, contributing to adverse outcomes including falls, fractures, and mortality. This “lung‑bone axis” concept highlights how chronic lung inflammation may affect skeletal health.¹³ Nevertheless, the direct molecular mechanisms linking COPD and OP remain poorly understood. Notably, a recent Mendelian randomization study found no direct causal relationship, suggesting that the observed association may be due to shared genetic effects or environmental confounders rather than direct causality.¹⁴ The integration of single-cell data analysis with MR analysis facilitates the systematic identification of comorbid genes and elucidates their critical roles within disease-associated pathophysiological networks.

Multi‑modal data integration has proven effective in dissecting complex clinical phenotypes.¹⁵ This study integrates MR with single‑cell transcriptomics to identify shared genetic determinants between COPD and OP—a comorbidity whose molecular basis remains largely unexplored. Beyond gene discovery, we systematically characterize the shared biological pathways, gene regulatory networks, disease associations, and chemical‑disease predictions, further validated by molecular docking and molecular dynamics simulations. This integrative approach not only prioritizes candidate genes for functional validation but also identifies potential chemical factors and offers a framework for dissecting shared pathogenesis in other comorbid conditions.

Materials and Methods

Study Design

Based on data from the GEO database, this study analyzed single‑cell transcriptome profiles from COPD and OP to identify common DEGs. MR was then used to prioritize causal genes among the shared DEGs, with cis‑eQTLs employed as instrumental variables. Subsequently, functional annotation, cell‑type mapping, and pseudo-time trajectory analysis were performed. Immune cell infiltration was assessed using CIBERSORT, followed by molecular docking and molecular dynamics simulations for the proteins encoded by the key genes and the associated chemicals. The workflow of the analysis is shown in Figure S1.

Data Source

Single-cell datasets pertaining to COPD and OP were obtained from the GEO database, utilizing the accession numbers GSE173896 and GSE147287, respectively. GSE173896 includes data from 5 COPD patients and 3 healthy controls, while GSE147287 contains data from 1 OP patient and 1 osteoarthritis control. To evaluate the impact of candidate genes on COPD and OP, MR analysis was performed using these genes as exposure factors and COPD and OP as outcomes. GWAS data for COPD-related genes were obtained from the ebi-a-GCST90018807 database, and for OP-related genes from finn-b-M13_OSTEOPOROSIS, both involving European populations.

Single-Cell Data Processing

This analytical pipeline is grounded in the established frameworks for scRNA‑seq data processing.¹⁶ Single-cell datasets were processed utilizing the Seurat software package.¹⁷ A Seurat object was created for each dataset, and a quality control filtering process was applied to ensure that each gene was expressed in at least three cells and that each cell expressed a minimum of 200 genes. The quality control thresholds were specifically defined for datasets GSE173896 and GSE147287, where cells were retained if the number of detected genes (nFeature_RNA) ranged from 200 to 5000, and mitochondrial gene expression (percent.mt) was below 15%. The NormalizeData function was used to perform data normalization. The FindVariableFeatures function, using the variance stabilization transformation (vst) method, identified highly variable genes by selecting the top 2000 variable genes.¹⁸ Expression values were standardized utilizing the ScaleData function prior to undergoing dimensionality reduction. Principal component analysis (PCA) was subsequently performed, with the significance of each principal component evaluated through the JackStraw and ScoreJackStraw functions. Clustering was executed at a resolution of 0.1, followed by the application of UMAP to visualize the distribution of cells within the reduced dimensional space. Cell type annotation was conducted using established marker genes in conjunction with the CellMarker 2.0 database, facilitating the classification of cell clusters into distinct immune populations.¹⁹ The expression patterns of marker genes within these cell populations were illustrated using heatmaps, and variations in cell type abundance between disease and control cohorts were analyzed.

DEGs and Enrichment Analysis

To identify DEGs across distinct cell populations, the FindMarkers function in Seurat was utilized. Genes that satisfied the criteria of |logFC| > 1 and an adjusted p-value < 0.05 were designated as DEGs. To explore common molecular characteristics between COPD and OP, co-regulated DEGs, encompassing both upregulated and downregulated genes, were selected for further investigation. Subsequently, functional enrichment analysis of these candidate genes was conducted using the clusterProfiler package.²⁰ This analysis incorporated Gene Ontology (GO) assessments—encompassing biological processes (BP), cellular components (CC), and molecular functions (MF)—alongside KEGG pathway analysis.

MR Analysis Identifies Key Genes

A two-sample MR analysis was performed utilizing the R software package, drawing upon publicly accessible summary statistics derived from GWAS available in the IEU OpenGWAS database. To ensure the validity and strength of the genetic instruments, for each candidate gene we selected cis‑eQTL SNPs with stringent criteria: a strong association with the exposure trait (P < 5e-8), independence after linkage disequilibrium clumping (r² < 0.1 within a 10,000‑kb window), and an F‑statistic > 10 (calculated as (β/SE)²) to exclude weak instruments.²¹ To deduce causal relationships between exposure and outcome traits, five MR analytical methods were employed: IVW, MR-Egger regression, weighted median, simple mode, and weighted mode approaches.^22–24 The IVW method was the main statistical approach, using a significance level of p < 0.05. Odds ratios (OR) and 95% confidence intervals were determined. Significant causal genes were those with consistent OR direction across five methods (all >1 or all <1), no heterogeneity (Cochran Q test P > 0.05), and no pleiotropy (MR-Egger intercept P > 0.05). For the significant genes identified through the aforementioned screening process, a series of analyses were conducted to ensure the robustness of the findings. A leave-one-out analysis was performed using the `mr_leaveoneout` function, which involved sequentially removing individual SNPs to evaluate the stability of their effects. Additionally, a single SNP Wald odds ratio analysis was carried out using `mr_singlesnp` to validate the independent causal effects of individual instrumental variables. Radial MR analysis was also employed using the RadialMR package to detect potential outliers. Further diagnostic plots were generated to support these analyses: scatter plots were used to illustrate the relationships between SNP effects and regression lines across various MR methods; forest plots were employed to visualize the effects of individual SNPs; funnel plots were utilized to assess potential publication bias; and leave-one-out sensitivity plots were created to depict the fluctuations in effects following the removal of individual SNPs. Collectively, these analyses and visualizations provided a comprehensive validation of the reliability of the causal associations.

Construction of Key Gene Regulatory Networks, Chemicals Prediction, and Disease Correlation

To elucidate the functional interactions among key genes, a gene-gene interaction (GGI) network was constructed utilizing the GeneMANIA platform (https://genemania.org/). This analysis identified the top 20 genes most closely associated with the key genes, as well as the seven most significantly enriched pathways pertinent to these genes. To further investigate post-transcriptional regulatory mechanisms, miRNA prediction was conducted using the miRDB database (https://www.mirdb.org)²⁵ This approach enabled the identification of common miRNAs targeting the key genes, thereby offering insights into potential regulatory interactions that may influence disease pathogenesis. Additionally, an analysis of gene associations with diseases and chemicals was performed by consulting the Comparative Toxicogenomics Database (https://ctdbase.org/). The top five diseases and top ten chemicals associated with the key genes were selected for visualization.

Acquisition of Putative Relevant Cells Types and Pseudo-Time Analysis

To evaluate the differential expression of key genes across various cell populations, a Wilcoxon test was conducted on all samples obtained from the single-cell datasets GSE173896 and GSE147287. For each identified cell type, the expression levels of key genes were compared between disease and control conditions. Cell types exhibiting a significance level of P < 0.05 were considered pertinent for further analysis. Subgroup analysis was conducted on the selected relevant cell clusters, involving data normalization, identification of highly variable genes, and re-execution of principal component analysis after isolating target cells. Nonlinear dimensionality reduction was achieved using the UMAP technique. Clustering analysis was performed using the FindClusters function (resolution=0.1) to elucidate cellular heterogeneity along disease trajectories. The FindAllMarkers function (min.pct=0.2, only.pos=TRUE) was utilized to identify marker genes specific to each subpopulation. To model the temporal evolution of cellular states, pseudo-time trajectory analysis was conducted using the Monocle package.²⁶ A set of high-variability genes was selected for analysis. The DDRTree algorithm was employed to perform dimensionality reduction and to reconstruct cellular developmental trajectories, while the orderCells function was utilized to compute pseudo-time values. The plot pseudo-time heatmap function was used to visualize dynamic gene expression patterns, thereby offering insights into the temporal regulation of genes during disease progression.

Correlation Analysis Between Key Genes and Immune Cell Infiltration

We investigated immune cell dynamics in both COPD and OP using a deconvolution approach based on CIBERSORT, which applies linear support vector regression to systematically assess cellular heterogeneity in expression matrices. Through CIBERSORTx (https://cibersortx.stanford.edu/), we obtained expression profiles for 22 immune cell types. Additionally, lollipop charts were constructed to illustrate the correlation patterns between key genes and key immune cell types.

Molecular Docking

The SDF files of the main active ingredients of the associated chemicals were obtained from the PubChem database (https://pubchem.ncbi.nlm.nih.gov/), and the structures of key target proteins were retrieved from the PDB database (https://www.rcsb.org/). Target proteins were optimized using PyMOL 2.1.0 by removing water molecules and small‑molecule ligands. AutoDock Tools 1.5.6 was then used to add hydrogen atoms and assign charges, and the results were saved in pdbqt format. The key target served as the receptor and its corresponding active ingredient as the ligand. Molecular docking was performed using vina 2.0 within the PyRx software to calculate binding energies and output the result files. PyMOL was used for visualization of the docking results. The binding energy (kcal/mol) reflects the stability of the ligand‑receptor interaction; a more negative binding energy indicates stronger and more stable binding. Visual analysis was further conducted using Discovery Studio 2020 client (https://discover.3ds.com/discovery-studio-visualizer-download).²⁷

Molecular Dynamics Simulation

To further evaluate the stability and dynamic behavior of the TUBB2A‑estradiol complex identified by molecular docking, molecular dynamics simulations were performed using Amber 24. The ff14SB force field was used for proteins, and GAFF for compounds. Each complex was solvated in a TIP3P water box with a buffer distance of 10.0 Å from the protein surface and neutralized with Na⁺/Cl⁻. Two‑stage energy minimization was performed, followed by gradual heating to 300 K and equilibration at 300 K and 1 bar. A 100 ns production run was conducted without restraints. Trajectory analysis (RMSD, Rg, SASA, RMSF, hydrogen bonds) was carried out using CPPTRAJ. Binding free energies were calculated by MM/GBSA from 100 frames extracted from the last 1 ns of the trajectory.

Statistical Analysis

R software (version 4.5.2) was used for statistical analyses, and Adobe Illustrator (version 2025) was employed for image editing. A p-value below 0.05 was considered statistically significant.

Results

Analysis of Single-Cell Data for COPD

Following quality control filtering, the GSE173896 dataset retained a total of 24,019 cells, comprising 11,225 control cells and 12,794 cells from COPD patients, and representing 23,342 genes (Figure S2A). Utilizing the FindVariableFeatures method, the top 2,000 highly variable genes were identified based on the correlation between mean expression and variance. From this subset, the ten most variable genes—such as SCGB1A1, IGKC, BPIFB1, and IGHG1—were selected for visualization (Figure 1A). Principal component analysis (PCA) demonstrated effective integration between control and COPD samples across the principal components (Figure S2B). Dimensionality reduction and clustering analyses indicated that the first 15 principal components captured substantial biological variation and were thus retained for further investigation (Figure S2C). Following dimensionality reduction, 18 distinct cell clusters were identified (Figure 1B and C). These clusters were subsequently classified into nine specific cell types: B cells, dendritic cells, endothelial cells, fibroblasts, macrophages, monocytes, NK cells, plasma cells, and T cells (Figure 1D and E). Analysis of cell type abundance indicated a higher proportion of NK cells in the control group, whereas T cells were prevalent in both groups (Figure 1F). A total of 2,623 genes were selected for further comparative analysis between the COPD and control groups (Figure 1G).

Figure 1 Analysis of single-cell data for COPD. (A) Identification of highly variable genes. (B) UMAP-based dimension reduction and cell clustering. (C) Grouped dimension-reduced clustering results. (D) Cell type annotation. (E) Cell annotation distribution across samples. (F) Cell type proportion statistics across samples. (G) Manhattan plot of DEGs by cell type.

Analysis of Single-Cell Data for OP

The GSE147287 dataset encompasses a total of 13,753 cells, consisting of 6,549 control cells and 7,204 cells from patients with OP, and includes 18,832 genes (Figure S3A). Utilizing the FindVariableFeatures method, the top 2,000 highly variable genes were identified based on the correlation between mean expression and variance. From this subset, the ten most variable genes—such as MPO, SPP1, DEFA3, and IGHG1—were selected for visualization (Figure 2A). Principal component analysis (PCA) demonstrated effective integration between control and OP samples across the principal components (Figure S3B). Dimensionality reduction and clustering analysis revealed that the first 15 principal components encapsulated significant biological information, leading to their selection for subsequent analyses (Figure S3C). After performing dimensionality reduction, the cell populations were delineated into 19 distinct clusters (Figure 2B and C). These clusters were subsequently categorized into seven primary cell types: B cells, dendritic cells, fibroblasts, macrophages, monocytes, plasma cells, and T cells (Figure 2D and E). Analysis of cell type abundance indicated a higher proportion of fibroblasts in the OP group samples, whereas monocytes were more prevalent in both groups (Figure 2F). A total of 2,454 genes were identified for further comparative analysis between the OP group and the control group (Figure 2G).

Figure 2 Analysis of single-cell data for OP. (A) Identification of highly variable genes. (B) UMAP-based dimension reduction and cell clustering. (C) Grouped dimension-reduced clustering results. (D) Cell type annotation. (E) Cell annotation distribution across samples. (F) Cell type proportion statistics across samples. (G) Manhattan plot of DEGs by cell type.

Identification of Candidate Genes

A comprehensive analysis identified 161 genes exhibiting co-upregulation in both COPD and OP, alongside 106 genes demonstrating co-downregulation in these conditions. Collectively, these 267 genes were earmarked as candidate genes for further investigation (Figure 3A and B). Gene Ontology (GO) enrichment analysis highlighted the most significantly enriched biological processes, including response to mechanical stimulus, extracellular matrix binding, and related functions (Figure 3C–E and Tables S1–S3). Furthermore, KEGG pathway enrichment analysis identified significant pathways, notably the MAPK signaling pathway, transcriptional dysregulation in cancer, malaria, cytoskeletal dynamics in muscle cells, among others (Figure 3F and Table S4).

Figure 3 Screening and functional analysis of candidate genes. (A) Co-upregulated genes. (B) Co-downregulated genes. (C) GO enrichment analysis of biological process. (D) GO enrichment analysis of cellular component. (E) GO enrichment analysis of molecular function. (F) KEGG enrichment analysis.

MR Analysis Identifies Key Genes

In the MR analysis, which employed candidate genes as exposure factors with COPD as the outcome event, 16 genes demonstrated significant causal effects on the development of COPD. Among these, SLC29A1, DAPK2, C2orf40, ECE1, PADI4, DSE, TUBB2A, and LYZ were determined to be risk factors (odds ratio [OR] > 1), whereas the remaining genes served as protective factors (OR < 1) (Figure 4A). Similarly, when examining OP as the outcome event, 23 genes exhibited significant causal relationships. Among these, RHOB, EMP1, FAM129B, RAB3IP, IER3, PADI4, PID1, TUBB2A, and EBPL were identified as risk factors (OR > 1), with the rest being protective (OR < 1) (Figure 4B). The directionality of OR values was consistently maintained across all five analytical methods employed (Table S5). No evidence of heterogeneity was observed (Cochran’s Q test, P > 0.05) (Table S6), and no horizontal pleiotropy was detected (MR-Egger intercept, P > 0.05) (Table S7). Notably, two genes, PADI4 and TUBB2A, were identified as risk factors for both COPD and OP. No protective factors common to both conditions were identified (Figure 4C and D). Further independent and detailed MR analyses were conducted for the significant genes PADI4 and TUBB2A. A series of diagnostic plots were generated to evaluate the data: scatter plots depicted the linear relationship between SNP exposure effects and outcome effects; forest plots presented independent causal estimates and aggregated effects for each SNP; funnel plots were utilized to assess potential publication bias; and a leave-one-out sensitivity analysis indicated no significant alterations in the overall estimates upon the exclusion of any single SNP, thereby confirming the robustness of the findings (Figure S4).

Figure 4 Causal effects of candidate genes in MR analysis. (A) Candidate genes on COPD outcomes. (B) Candidate genes on OP outcomes. (C) Genes with shared positive effects. (D) Genes with shared negative effects.

Construction of Key Gene Regulatory Networks, Chemical Prediction, and Disease Correlation

A Gene-Gene Interaction (GGI) network analysis was conducted on two pivotal genes, PADI4 and TUBB2A, to investigate their functional associations. This analysis identified several genes closely related to these key genes, including SIRT2, PADI3, and TUBA1A, among others. The pathways enriched in association with these genes encompassed hydrolase activity and purine ribonucleoside binding (Figure 5A). Gene regulatory network analysis revealed that PADI4 is predicted to interact with 34 microRNAs and TUBB2A with 19 microRNAs in the miRDB database (Figure 5B). These genes are also linked to potential Chemical interactions with BPA, TCDD, and estradiol (Figure 5C). Furthermore, disease prediction analysis connects these genes to conditions like weight loss, chemical and drug-induced liver injury, and necrosis (Figure 5D).

Figure 5 Key analysis of gene regulatory networks and target prediction. (A) Gene-gene interaction network. (B) miRNA regulatory networks. (C) Predicting potential associated chemicals. (D) Identifying associated diseases.

Putative Relevant Cell Types Acquisition and Pseudo-Time Analysis

To identify cell types linked to COPD, we analyzed the GSE173896 dataset, finding a significant enrichment of NK cells in both genes, highlighting them as the key cell population (Figure 6A). Further clustering revealed five NK cell subpopulations (Figure 6B), with Subpopulation 0 as the progenitor, differentiating into other types over time (Figure 6D). We examined prognostic gene expression during NK cell differentiation, revealing two distinct expression clusters across stages (Figure 6C).

Figure 6 Identification of COPD-related differential cells and pseudo-time analysis. (A) Expression profiles of PADI4 and TUBB2A across samples. (B) Clustering analysis of NK cells. (C) Heatmap showing gene expression dynamics over pseudo-time. (D) Pseudo-time trajectory analysis of NK cells. *** represents P < 0.001, ** represents P < 0.01, * represents P < 0.05 and, ns represents no significance.

To identify cell types linked to OP, we analyzed the GSE147287 dataset and found monocytes significantly enriched in both genes, marking them as the key cell population (Figure 7A). We further divided the monocytes into seven subpopulations (Figure 7B), with Subpopulation 3 identified as the progenitor cell initiating the developmental trajectory, eventually differentiating into various monocyte types (Figure 7D). We then examined the expression of key genes during this differentiation, finding two distinct expression clusters across different stages (Figure 7C).

Figure 7 Identification of OP-related differential cells and pseudo-time analysis. (A) Expression profiles of PADI4 and TUBB2A across samples. (B) Clustering analysis of monocytes. (C) Heatmap showing gene expression dynamics over pseudo-time. (D) Pseudo-time trajectory analysis of monocytes. *** represents P < 0.001, * represents P < 0.05 and, ns represents no significance.

Immune Cell Infiltration

Using the COPD dataset GSE76925 and the OP dataset GSE56815 from the GEO database, we performed immune cell infiltration analysis to characterize the immune landscape of both diseases. The proportion of activated NK cells was higher in the COPD cohort than resting NK cells (Figure 8A). Both PADI4 and TUBB2A showed strong correlations with activated NK cells (Figure 8C and D). In the OP cohort, although monocytes did not differ significantly between disease and control groups, they comprised a high proportion of immune cells, suggesting their important role in OP (Figure 8B). Furthermore, PADI4 was also strongly correlated with monocytes in the OP cohort (Figure 8E). However, TUBB2A did not show a strong correlation with monocytes in the OP cohort (Figure 8F). These results corroborate our previous findings.

Figure 8 Disease‑related immune characteristics and immune infiltrating cells associated with key genes. (A) Analysis of 22 key immune infiltrating cell types in COPD. (B) Analysis of 22 key immune infiltrating cell types in OP. (C) Lollipop plot showing correlation between PADI4 and immune infiltrating cells in COPD. (D) Lollipop plot showing correlation between TUBB2A and immune infiltrating cells in COPD. (E) Lollipop plot showing correlation between PADI4 and immune infiltrating cells in OP. (F) Lollipop plot showing correlation between TUBB2A and immune infiltrating cells in OP. **** represents P < 0.0001, *** represents P < 0.001, ** represents P < 0.01 and, * represents P < 0.05.

Molecular Docking Analysis

The molecular binding affinities between the associated chemicals and target proteins are shown in Table 1. A lower binding energy indicates a stronger interaction. Generally, a value below −5.0 kcal/mol suggests good binding activity, while below −7.0 kcal/mol indicates strong binding activity.^28,29 Notably, all binding energies were below −5.0 kcal/mol, demonstrating strong binding affinity. Among them, the binding between TUBB2A and estradiol was below −7.0 kcal/mol, indicating an extremely stable interaction. Figure 9 presents detailed molecular interactions, including hydrogen bonds, π‑π stacking, alkyl interactions, and π‑alkyl interactions.

Table 1 Docking Results Between Target Proteins and the Associated Chemicals

Figure 9 Molecular docking and 2D interaction map of the associated chemicals with target proteins. (A) PADI4 binds to BPA. (B) TUBB2A binds to BPA. (C) PADI4 binds to TCDD. (D) TUBB2A binds to TCDD. (E) PADI4 binds to Estradiol. (F) TUBB2A binds to Estradiol.

Molecular Dynamics Simulation

Based on these docking findings, we selected the TUBB2A–estradiol complex for further molecular dynamics simulations, as it exhibited the lowest binding energy among all docking results, suggesting a potentially strong and functionally relevant interaction. RMSD and Rg analyses indicated limited structural fluctuation and a compact complex conformation (Figure 10B and E). RMSF values suggested that the amino acid residues remained stable throughout the simulation (Figure 10C). SASA analysis revealed changes in solvent interactions, while hydrogen bond analysis confirmed the stability of hydrophilic binding (Figure 10A and D). MM/GBSA analysis further quantified the binding free energy, showing that van der Waals forces and electrostatic interactions were the primary drivers of complex stability (Figure S5). Free energy landscape analysis identified stable conformations of the complex within specific ranges of RMSD and Rg (Figure 10F). Collectively, these results provide insights into the binding mechanism of the target protein–chemical complex.

Figure 10 Molecular dynamics simulation of the TUBB2A‑estradiol complex. (A) Hydrogen bond number over time. (B) RMSD curve. (C) RMSF curve. (D) SASA curve. (E) Rg curve. (F) Free energy landscape.

Discussion

Our study identifies PADI4 and TUBB2A as co-pathogenic genes associated with COPD and OP, offering novel insights into the pathogenesis of these conditions. The PADI4 gene encodes Peptidyl Arginine Deiminase 4 (PAD4), a critical member of the Peptidyl Arginine Deiminase (PAD) family. PADI4 is instrumental in various physiological and pathological processes, particularly in the post-translational modification of proteins through citrullination. Citrullination, which involves the conversion of arginine residues, is a crucial modification with significant implications for gene regulation and the formation of NETs.^30,31 PADI4 is implicated in the pathogenesis of several diseases by modulating NETs formation, including heart failure, acute lung injury, rheumatoid arthritis, lupus nephritis, asthma, and abdominal aortic aneurysm.^32–37 In the field of oncology, PADI4 has been identified as significantly linked to tumorigenesis and cancer progression. Through the process of citrullination, PADI4 modulates gene expression, thereby affecting tumor cell proliferation, migration, and invasiveness.^38,39 Additionally, the function of PADI4 under hypoxic conditions merits further investigation. In cases of breast cancer and hepatocellular carcinoma, PADI4 expression is regulated by hypoxia-inducible factors (HIFs) and contributes to the transcriptional activation of hypoxia response elements (HREs).⁴⁰ While current literature does not directly link PADI4 to COPD or OP, its role in inflammation, cell differentiation, and cancer progression implies possible indirect connections to these conditions’ pathophysiology.

The β-tubulin encoded by the TUBB2A gene is an essential component of microtubules, which are integral to the formation and maintenance of the cytoskeleton. Microtubules serve as essential components of the cytoskeleton and play critical roles in multiple cellular activities such as mitotic division, cellular motility, and vesicular trafficking.⁴¹ The relevance of this protein to disease has attracted considerable attention in recent research, particularly in the context of neurodevelopmental disorders. Mutations in the TUBB2A gene have been identified as exceedingly rare causes of childhood neurodevelopmental disorders. These mutations often lead to severe brain developmental abnormalities, such as microcephaly, global developmental delay, language impairment, seizures, and profound behavioral disorders.⁴² Moreover, several studies have associated TUBB2A mutations with infantile-onset epilepsy and developmental delay, highlighting the gene’s essential role in neurodevelopment.^43,44 Although no direct studies link TUBB2A to OP, its role in microtubule and cytoskeletal functions hints at a possible impact on bone demineralization. Future research should integrate genetic, molecular, and clinical methods to explore this further. Additionally, the influence of TUBB2A variants on the risk of developing COPD and its exacerbations is still uncertain, necessitating more comprehensive genetic research to clarify its role in COPD.

Our observations indicate that weight loss, along with chemical and drug-induced liver injury and necrosis, as well as exposure to substances such as BPA, TCDD, and estradiol, are frequently linked to COPD and OP. Effective weight management in COPD patients is a critical determinant of their prognosis and quality of life. Research suggests that a low BMI and weight loss significantly elevate the mortality risk in individuals with COPD.⁴⁵ Moreover, weight loss is strongly associated with frailty in COPD patients, which increases the likelihood of acute exacerbations and hospitalizations and significantly raises all-cause mortality.⁴⁶ In the context of OP, particularly among the elderly population, weight loss has been associated with a decrease in bone mineral density, thereby elevating the risk of fractures. This phenomenon is likely attributable to bone loss resulting from diminished mechanical loading.⁴⁷ Psoraleae Fructus (PF), a traditional Chinese medicinal remedy, is commonly employed in the treatment of OP. However, research suggests that PF and its constituents may provoke specific hepatotoxic effects by enhancing the activation of the NLRP3 inflammasome.⁴⁸ Furthermore, recent evidence highlights the pivotal role of TNF-α in the pathophysiology of COPD. Targeting TNF-α and its associated signaling pathways has been shown to effectively attenuate the inflammatory response in COPD.^49,50

BPA, recognized as an endocrine disruptor, has been significantly associated with oxidative stress and inflammation in patients with COPD.^51,52 In the context of OP, BPA primarily disrupts hormonal balance during bone remodeling by mimicking estrogenic effects, thereby compromising the structure and function of bone tissue.⁵³ Furthermore, BPA induces oxidative stress and inflammatory responses, leading to osteocyte apoptosis and autophagy, which further destabilizes bone homeostasis.^54,55 Additionally, BPA has been shown to induce ferroptosis in osteoblasts via the p53/SLC7A11 axis, thereby inhibiting osteoblast proliferation and differentiation and exacerbating osteoporosis.⁵⁶ Similarly, research has demonstrated that exposure to TCDD adversely affects bone strength, structure, and density, with particularly significant impacts on bone mineralization and geometric architecture. These changes may contribute to the development or worsening of OP.⁵⁷ Mechanistically, TCDD significantly and dose‑dependently decreases the mRNA levels of osteoblast differentiation markers such as RUNX2, alkaline phosphatase, and osteocalcin, while also reducing the number of TRACP+ multinucleated osteoclasts, F‑actin rings, and the area of bone resorption, effects that are mediated by the aryl hydrocarbon receptor (AhR).⁵⁸ Additionally, studies have shown that estrogen supplementation can effectively mitigate COPD-like symptoms induced by endocrine disruptors. This suggests that estrogen may have anti-inflammatory properties by modulating the expression of estrogen receptor alpha (ERα) and nuclear factor kappa-B (NF-κB).⁵⁹ Furthermore, the role of estrogen in COPD extends beyond the regulation of inflammation. Recent research indicates that estrogen may also influence depressive symptoms associated with COPD by modulating oxidative stress and neuroinflammation.⁶⁰ In the context of osteoporosis, estradiol has been shown to exert protective effects on bone by promoting osteoprotegerin (OPG) as a decoy receptor for RANKL, thereby inhibiting osteoclast activation, and by activating the Wnt signaling pathway to enhance osteoblast survival and function.⁶¹ These findings, together with molecular docking and molecular dynamics simulations demonstrating stable binding between the associated chemicals and target proteins, suggest that BPA, TCDD, and estradiol may contribute to COPD–OP comorbidity through their association with PADI4 and TUBB2A; however, this hypothesis requires further investigation.

Through single-cell analysis and immune cell infiltration analyses, our study corroborates that NK cells constitute a potentially significant cell type in the pathology of COPD, although their precise role necessitates further functional validation. As an integral component of the innate immune system, the functions and mechanisms of NK cells in COPD are progressively being elucidated. Current research demonstrates that NK cells exhibit notable activation in both peripheral blood and lung tissue of COPD patients, with a strong correlation to disease severity. Notably, in advanced stages of the disease, there is a marked increase in activated Nkp44+ NK cells.^62,63 Moreover, the role of NK cells in COPD extends beyond the maintenance of inflammation; they may also be involved in the immunopathological processes of the disease. Single-cell RNA sequencing analysis has identified COPD-associated NK cell signature genes, which could serve as potential targets for further mechanistic studies and the development of novel therapeutics.⁶⁴ Moreover, NK cells display increased proportions in the bronchoalveolar lavage fluid of patients with COPD, a phenomenon that persists even after smoking cessation. This observation suggests a potential role for NK cells in the ongoing progression of the disease.⁶⁵ Additionally, antiviral responses mediated by NK cells in COPD patients exhibit abnormalities. Research indicates that tissue-resident CD49a+ NK cells show dysregulated responses to viral infections, which may lead to excessive inflammatory reactions during viral episodes, thereby exacerbating COPD symptoms.⁶⁶ Furthermore, chronic smoking results in decreased interleukin-16 (IL-16) concentrations within NK cells, a change associated with the involvement of oxygen radicals, further underscoring the complex role of NK cells in the pathophysiology of COPD.⁶⁷

Through single‑cell and immune cell infiltration analyses, our study indicates that monocytes are a potentially relevant cell type in the pathogenesis of OP, although their precise role requires further functional validation. As precursor cells to osteoclasts—the primary cells responsible for bone resorption—monocytes are integral to the onset and progression of OP. Initially, monocytes contribute to the pathological process of OP through the secretion of various cytokines. Research indicates that monocytes can secrete pro-inflammatory cytokines such as TNF-α, IL-1, and IL-6, which facilitate osteoclast formation and activation, thereby accelerating bone resorption.^68,69 Furthermore, monocytes enhance osteoclast differentiation and bone resorption by expressing receptor activator of nuclear factor kappa-B ligand (RANKL), which interacts with osteoclast precursor cells.⁷⁰ Additionally, the involvement of monocytes in OP is associated with their activity within the circulatory system. Studies have demonstrated that circulating monocytes can impact bone density through specific gene expression profiles. For example, the upregulation of the ANXA2 gene in individuals with low bone density is significantly correlated with the onset of OP.⁷¹ Furthermore, the phosphorylation status of certain proteins in monocytes, such as HSP27, is implicated in the pathogenesis of OP, as these proteins facilitate monocyte migration and osteoclast formation.⁷² The role of monocytes in OP is also modulated by the immune system. Other immune cells, including T cells and B cells, affect monocyte function through cytokine secretion and direct cell-cell interactions.^73,74 Additionally, the interactions between monocytes and other immune cells may influence bone metabolism by altering the inflammatory environment within the bone.^75,76 In the field of OP treatment research, the regulation of monocyte differentiation and function has become a pivotal area of focus. Empirical studies have demonstrated that specific pharmacological agents can inhibit the differentiation of monocytes into osteoclasts through the modulation of gene expression. These findings have shown significant anti-osteoporotic effects in both in vitro and in vivo experimental settings.⁷⁷

COPD and OP are two clinically distinct diseases, yet their comorbidity is common and remains poorly understood at the molecular level. By identifying shared genetic determinants (PADI4, TUBB2A) and associated chemicals (BPA, TCDD, estradiol), our findings provide a basis for risk stratification and offer potential avenues for improving the management of COPD‑OP comorbidity through targeting these key genes and chemical factors. This study thus provides a framework for dissecting the molecular links between seemingly unrelated comorbid conditions. Despite these significant findings, several limitations warrant consideration. First, the sample size of the single-cell dataset used in this study is small, which restricts the reproducibility and generalizability of the results. Second, the single-cell data were obtained from individual nucleated cells in lung tissue and bone marrow from distinct patient cohorts. Although this cross-tissue, cross-individual design provides a promising approach for investigating comorbidity mechanisms, it does not fully account for potential confounding effects due to variations in tissue-specific biological contexts and inter-individual heterogeneity. Future research integrating matched multi-tissue samples from the same patient cohort will be essential for further validating and refining the core genes identified in this study. Additionally, the data utilized in the Mendelian randomization study were derived from European populations, which may restrict the generalizability of our findings to other ethnic groups. Furthermore, our MR analysis did not thoroughly assess the directionality between exposure and outcome; therefore, the possibility of reverse causality cannot be completely ruled out, and these limitations may affect the reliability of causal inference. Therefore, these findings should be regarded as exploratory.

Conclusion

This study combined single-cell data and MR analysis to uncover a shared genetic basis between COPD and OP, highlighting PADI4 and TUBB2A as key genes linked to both conditions. NK cells are primarily associated with COPD, while monocytes are linked to OP. The MAPK signaling pathway was the most significant according to KEGG enrichment analysis. In addition, our analysis indicates that BPA, TCDD, and estradiol may contribute to the pathogenesis and treatment of both diseases. Molecular docking and molecular dynamics simulations further complemented our findings. These findings provide new insights into the COPD‑OP comorbidity and suggest directions for future mechanistic research.