SPG7-Mediated Regulation of mPTP and Mitochondrial Flickering in COPD: A Bioinformatics-Based Prediction of Mechanistic Framework

Introduction

Asthma and chronic obstructive pulmonary disease (COPD) constitute the two most critical components of the chronic respiratory disease (CRD) burden. In the context of global population aging, COPD is the primary driver of the CRD burden.¹ Despite deepening insights into the molecular mechanisms underlying COPD, disease progression remains closely associated with elevated mortality rates during acute exacerbations.² Deaths attributed to COPD rank third among chronic disease-related mortality, preceded only by cardiovascular diseases and stroke. Concurrently, COPD is the fifth leading disabling disease globally.³ Furthermore, our group’s previous epidemiological research¹ identified an accelerated upward trend in the prevalence and incidence of COPD among middle-aged and elderly populations in low socio-demographic index (SDI) and high SDI regions, respectively. This indicates that post-pandemic pulmonary damage and global population aging will become global public health priorities in the coming decades, necessitating targeted interventions to mitigate the sustained damage inflicted by COPD on patients worldwide during the process of pulmonary aging.⁴

COPD is a progressive aging-associated disease characterized by recurrent alternations between acute exacerbations and stable states as the disease progresses. Within the stable phase, the GOLD guidelines⁵ further classify stable COPD (STCOPD) into three stages: A, B, and E. Each stage signifies an aggravation of the patient’s clinical symptoms and reflects the progressive destruction of airway structures and the aging of pulmonary tissue in patients with COPD. Mounting evidence indicates that the stage-wise progression of COPD is closely associated with the process of pulmonary aging.^6–8 This process predisposes patients with COPD to at least two extrapulmonary comorbidities,⁹ including osteoporosis, skeletal muscle dysfunction, and cognitive impairment. This systemic aging process is highly correlated with the oxidative stress experienced by patients with COPD.¹⁰ It is primarily manifested by the exacerbating effect of low-grade inflammatory infiltration—driven by the senescence-associated secretory phenotype (SASP)—on COPD,¹¹ as well as a vicious cycle wherein oxidative stress further amplifies inflammation-induced cell death.¹² Therefore, mitigating oxidative stress, enhancing energy metabolism, delaying cellular senescence, and preventing the destruction of pulmonary architecture caused by acute exacerbations hold the potential to open a new window for the comprehensive management of COPD.¹³

Oxidative stress is widely recognized as a crucial factor driving COPD progression, and both exogenous and endogenous reactive oxygen species may be implicated in the pathogenesis of COPD.¹³ Although inflammatory responses induced by cigarette smoke are considered the primary source of oxidative stress in COPD, studies have shown¹⁴ that oxidative stress persists even after smoking cessation. This implies that ROS leakage originating from mitochondria exerts a profound impact on the pathological mechanisms of COPD. Mitochondrial DNA is 30 times more susceptible to oxidants than nuclear DNA.¹³ Consequently, mounting evidence indicates that ROS leakage from the respiratory chain also leads to extensive mitochondrial dysfunction in COPD, thereby triggering chronic pulmonary inflammation.¹⁰ Mitochondria are characterized by their dynamic stability and systemic nature. Therefore, mitochondrial dysfunction is also a systemic phenomenon¹³ that affects multiple cell types, resulting in skeletal muscle atrophy and respiratory muscle dysfunction. For patients with COPD, pulmonary cells are the first to bear the burden of mitochondrial oxidative stress. Furthermore, mitochondrial dynamics regulate the onset and progression of COPD. Numerous studies have demonstrated that distinct states of mitochondrial compensation or decompensation manifest across the early stages, late stages, and acute exacerbations of COPD (AECOPD).³ These state alterations are primarily reflected in aberrant changes in mitochondrial quality control, encompassing mitophagy, fusion, and fission.³ As the largest intracellular energy factories,¹⁵ proper mitochondrial quality control serves as a critical pathway for maintaining ionic homeostasis, averting oxidative stress, and preventing intrinsic cellular apoptosis.^15,16 However, mitochondrial quality control is an adenosine triphosphate (ATP)-dependent process; ultimately, cell fate is determined by the balance between the demand for substrate clearance and the cell’s tolerance capacity, and is constrained by the energy supply.^17,18 These lines of evidence also elucidate why the majority of studies failing to stratify by COPD status present a typical dual effect: the simultaneous presence of insufficient mitophagy alongside the impacts of excessive mitophagy^19,20 (Figure 1).

Figure 1 Key questions proposed based on mitochondria-related mechanisms and phenomena in COPD: COPD is a stage-specific age-related lung disease, characterized by stage-specific exacerbation and systemic aging. This alteration is highly associated with mitochondrial status. The role of mitochondrial membrane potential in maintaining mitochondrial quality control directly determines intrinsic cellular apoptosis. However, the direction of mitophagy remains controversial across all stages of COPD during this process. Based on this, we propose three key questions. ↑: Directional transitions/ directional progression. Created with BioRender.com.

Based on the aforementioned phenomena, our research group has raised three questions as illustrated in Figure 1, and formulated the hypothesis presented in Figure 2: “The fundamental purpose of mitophagy is to maintain cell survival at the expense of self-sacrifice. If mitochondrial function can be restored and mitochondrial quality control bypassed while maintaining cell survival, thereby preserving a high level of cellular metabolic activity, it would contribute to delaying COPD-induced pulmonary aging”. This conjecture is well-founded, as the majority of proper mitochondrial functions rely on the stability of the mitochondrial inner membrane potential.¹⁶ Once the membrane potential is lost, mitochondria do not immediately initiate mitophagy; rather, there is a delay of approximately 30 minutes in Parkin activation. During this period, if the membrane potential is restored, PINK1 activation will be rapidly reversed,²¹ and oxidative phosphorylation will resume. The mitochondrial membrane potential is strictly governed by the ion transport processes across the mitochondrial membrane,¹⁶ and the mitochondrial permeability transition pore (mPTP) serves as a critical regulatory factor in maintaining this potential.²² Such mPTP-mediated “mitochondrial flickering” has been proven to play a decisive role in cell survival and apoptosis. The reversible opening of the mPTP directly regulates ROS and Ca2+ homeostasis, and mitigates the formation of mtROS through transient depolarization, respiratory chain activation, and accelerated electron flow.²³ Therefore, as the initial step in a series of investigations, this study employs bioinformatics approaches (Figure 3) to examine the differential expression of mPTP-related genes across distinct disease states (STCOPD and AECOPD). While identifying the key genes that regulate “ mitochondrial flickering”, we also seek to evaluate the plausibility of the aforementioned hypothesis. Through this research, we anticipate discovering novel pathways that regulate the mitochondria-cell survival framework, thereby providing new insights and directions for future related studies.

Figure 2 Mitochondria-cell survival architecture based on the hierarchical theory of COPD: Based on previous controversies and existing evidence, we propose a novel theory. We hypothesize that mitochondrial flickering plays a key role in regulating mitophagy. Enhancing mitochondrial flickering maintains the stability of mitochondrial membrane potential and regulates mitochondria-cell ion homeostasis. This process exerts at least two functions. First, during the stage when cells are under excessive stress and on the verge of death, this function can sustain cell survival by promoting mitophagy. Second, during the chronic inflammatory phase, it protects mitochondria from mitophagic damage, thereby maintaining cells in a high-energy metabolism state. Based on this, we further propose targeted intervention strategies for different stages of COPD. ↑: the upregulation or downregulation of ROS. ΔΨm: directional transitions/ directional progression. Created with BioRender.com.

Figure 3 Analysis flow chart. ↑: Directional transitions/ directional progression.

Methods

Ethical Exemption

Ethical approval was waived for this study because all analyses were based on publicly available summary statistics databases, did not involve personal data, and were entirely anonymized. This study also obtained an ethical exemption in accordance with local policies and complies with Items 1 and 2, Article 32 of the Measures for Ethical Review of Life Science and Medical Research Involving Humans: (1) conducting research utilizing legally obtained public data or data generated through observation without interfering with public behavior; and (2) conducting research using anonymized information and data.

Data Sources

The objective of this study is to reveal the alterations in how differentially expressed genes regulate mPTP function in STCOPD and AECOPD, thereby inferring the critical role of mPTP-mediated ΔΨm changes in cellular life-and-death decisions during COPD state transitions. Therefore, the data sources must explicitly indicate whether the samples were derived from patients with COPD in the stable phase or the acute exacerbation state. According to the GOLD guidelines, the deterioration of AECOPD symptoms occurs within 14 days;⁵ thus, we defined patients with COPD experiencing acute exacerbations within 14 days as being in the AECOPD state. We retrieved the expression data of COPD-related genes under different states from the Gene Expression Omnibus (GEO) database (https://www.ncbi.nlm.nih.gov/geo/), including GSE148871, GSE148877, and GSE112165. Among them, GSE148871 and GSE148877 originated from the same study by the same research team, containing 304 and 209 samples, respectively. GSE148871 is an intervention study on AECOPD, whereas GSE148877 is an intervention study on STCOPD. We selected the baseline sputum samples prior to intervention from these two datasets as the training set, ultimately obtaining 28 STCOPD samples and 44 AECOPD samples. Subsequently, we retrieved the baseline and day-84 post-exacerbation sputum and blood samples from GSE148871 as the validation set, which included 16 STCOPD sputum samples, 21 STCOPD blood samples, 44 AECOPD sputum samples, and 49 AECOPD blood samples. Furthermore, GSE112165 was also utilized for the validation set, providing 17 STCOPD sputum samples and 13 AECOPD sputum samples, respectively. All three aforementioned datasets utilized the same GPL570 platform. To investigate the differential expression of target genes at the cellular level, we conducted an in-depth analysis using the GSE268542 dataset, which is based on 10x single-cell RNAseq technology and the GPL21697 platform, comprising 1 STCOPD lung tissue cell sample and 1 Influenza A Virus (IAV)-infected COPD lung tissue cell sample. Finally, we retrieved 2549 mPTP-related genes through the GeneCards database (https://www.genecards.org/), and after selecting those with a relevance score higher than 20, we ultimately obtained 196 mPTP-related genes. Further details regarding sample sources and the inclusion information for group integration can be found in Figure 3 and Table S1.

Differential Expression and Enrichment Analyses

To identify the differentially expressed genes (DEGs) between STCOPD and AECOPD, the raw expression matrix of the training set was normalized prior to differential analysis to ensure the comparability of sample baselines. Given that the two datasets originated from parallel cohorts of the same research team and the disease groupings were completely mutually exclusive, conventional batch effect correction algorithms would eliminate genuine biological variations. Therefore, following the methodological design of the datasets’ source literature,²⁴ compulsory batch effect removal was not performed in this study. Subsequently, the limma package (v 3.64.3) was utilized to analyze the normalized matrix. DEGs were obtained based on the thresholds of |log2FC| > 1 and a false discovery rate (FDR, ie, adjusted P-value) < 0.05 following multiple testing correction via the Benjamini-Hochberg (BH) method. Thereafter, enrichment analysis of the DEGs was conducted utilizing the Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) databases, along with R packages such as clusterProfiler (v 4.16.0), to delineate the primary biological alterations during the transition from STCOPD to AECOPD. GO analysis revealed the significantly enriched biological processes (BP), cellular components (CC), and molecular functions (MF) of the DEGs, whereas KEGG analysis elucidated the signaling pathways involving the key target genes. The significance threshold for both enrichment analyses was set at an FDR (adjusted P-value) < 0.05 based on the BH correction method.

Screening of Key mPTP-Related Genes Among DEGs

To identify the key genes capable of regulating mPTP during disease progression, we first intersected the DEGs with the 196 genes obtained from the GeneCards database. Subsequently, the identified differentially expressed mPTP-related genes (DE-MPTPRGs) were integrated into the Search Tool for the Retrieval of Interacting Genes/Proteins (STRING) database (https://cn.string-db.org/) to construct a protein-protein interaction (PPI) network. A confidence score threshold was set to 0.9 to ensure the reliability of the included interactions.

To ensure the clinical significance of the key genes among the DE-MPTPRGs, we employed a robust Mendelian randomization (MR) framework for causal inference. In this analysis, DE-MPTPRGs were regarded as the exposure, whereas acute hospitalization for COPD served as the outcome. The cis-expression quantitative trait loci (cis-eQTL) data for the DE-MPTPRGs were derived from whole blood samples in the eQTLGen Consortium (https://www.eqtlgen.org/) and lung tissue samples from the Genotype-Tissue Expression (GTEx) database (Version V11; https://gtexportal.org/). The COPD acute hospitalization events were acquired from the IEU OpenGWAS database (https://gwas.mrcieu.ac.uk/), specifically utilizing the whole blood sample dataset labeled “finn-b-COPD HOSPITAL”. This dataset encompasses 16,380,466 single-nucleotide polymorphisms (SNPs) across 218,792 individuals, including 6500 COPD hospitalization events and 212,292 healthy European controls. Further information regarding data sources can be found in Table S2.

The MR analysis must satisfy three core assumptions: (1) a significant and robust association exists between the instrumental variables (IVs) and the exposure; (2) the IVs are independent of confounders; and (3) the IVs affect the outcome solely through the exposure pathway. Therefore, to screen for IVs, we initially utilized the TwoSampleMR package (v 0.6.22),²⁵ setting a significance threshold of . Then, the linkage disequilibrium (LD) effects of each SNP were removed based on the intervals of and . The LD reference panel used was “EUR”. Finally, the robustness of the IVs was verified by calculating the F-statistic, with a value exceeding 10 indicating reliable IVs. Subsequently, to satisfy the exclusivity assumption, we set an exclusion threshold of . Concurrently, non-palindromic SNPs were harmonized, and all SNPs containing palindromic sequences were excluded to ensure proper allele alignment between the exposure and outcome datasets. Notably, when processing the lung tissue samples derived from the GTEx database, we set the intervals to and to remove the LD effects of each SNP.

Various causal inference methods were employed for the MR analysis, including the MR-Egger, weighted median, inverse variance weighted (IVW),²⁶ simple mode, and weighted mode methods, although the conclusive results primarily relied on the IVW method. To control for the risk of false positives induced by multiple hypothesis testing, the BH method was applied for multiple testing correction on the causal inference P-values of all candidate genes. Under this method, an adjusted P-value (FDR) <0.05 was established as the criterion for determining a significant causal relationship, and the effect model was selected based on the results of the heterogeneity test. If the heterogeneity P-value was <0.05, a (multiplicative) random-effects model was employed; otherwise, a fixed-effects model was utilized. Furthermore, sensitivity analyses were required to ensure the robustness of the results, which encompassed heterogeneity assessment,²⁷ pleiotropy testing,²⁸ and the leave-one-out (LOO) method.²⁹ P-values <0.05 in the heterogeneity and pleiotropy tests indicated the presence of heterogeneity among the IVs and horizontal pleiotropy between the exposure and outcome variables, respectively. The LOO method sequentially removed each SNP and re-evaluated the results to ensure minimal variation in the outcome, thereby accounting for the impacts introduced by confounders.³⁰

The selected DE-MPTPRGs, while validated by sensitivity analysis to have a causal relationship with COPD hospitalization events, must simultaneously exhibit a consistent expression trend with that in the differential gene analysis. Only the genes ultimately satisfying these criteria will be designated as the core genes among the DE-MPTPRGs.

After identifying the core genes among the DE-MPTPRGs, we performed a colocalization analysis on the key biomarkers utilizing the coloc R package.³¹ To ensure the inclusion of all phenotype-associated eQTLs, we expanded the eQTL screening window to 500 kb upstream and downstream. For each locus, a Bayesian approach evaluated the degree of support for the following five mutually exclusive hypotheses:³² (1) neither trait has an association; (2) only trait 1 has an association; (3) only trait 2 has an association; (4) both traits have an association, but each corresponds to a different causal variant; and (5) both traits have an association and share the same causal variant. This analysis provided a posterior probability for each hypothesis tested (H0, H1, H2, H3, and H4). We set the prior probability of a SNP being associated only with trait 1 () to ; the prior probability of a SNP being associated only with trait 2 () to ; and the prior probability of a SNP being associated with both traits () to . If the posterior probability of a shared causal variant () is , the two signals are considered to demonstrate strong evidence of colocalization. Moderate evidence of colocalization was defined as Further information regarding the datasets utilized for the Mendelian randomization and gene colocalization analyses can be found in Table S3.

Validation of Predictive Models and Functional Annotation

After identifying the core genes among the DE-MPTPRGs, we validated the model across three datasets. The rms package (v 8.0–0) was utilized to construct a diagnostic nomogram incorporating the key genes, and calibration curves were established to evaluate the predictive performance of the nomogram. Subsequently, we plotted receiver operating characteristic (ROC) curves using the pROC package (v 1.19.0.1) to quantify the overall discriminative ability of the nomogram across all possible thresholds. Furthermore, decision curve analysis (DCA) was performed to determine the clinical utility of the nomogram. To further validate the expression trends of the key genes across different datasets, boxplots were employed for visual representation, and Pearson correlation coefficients were calculated using the psych package (v 2.5.6) to assess the synergistic interactions among the key biomarkers in the distinct validation sets.

To evaluate the biological significance of the alteration trends in the key genes, Gene Set Enrichment Analysis (GSEA) was conducted for each biomarker utilizing the clusterProfiler package, based on the MF dataset from the GO database. The enriched pathways were ranked according to their normalized enrichment scores (NES), with thresholds of and P<0.05 established as the criteria for statistical significance. This approach facilitated the identification of molecular functions significantly associated with the expression profiles of the key genes.

Immune Infiltration Analysis

Regardless of the COPD stage, the disease is invariably accompanied by persistent inflammatory responses and continuous immune stress. Therefore, to evaluate the immunological differences in the AECOPD state compared to the STCOPD state, we performed immune cell profiling utilizing the single-sample Gene Set Enrichment Analysis (ssGSEA) algorithm within the GSVA package (v 2.2.0).³³ This method computes the immune cell infiltration score for each sample based on gene sets comprising immune cell markers.³⁴ Subsequently, immune cell types exhibiting significant variations were identified by comparing the immune cell scores between the STCOPD and AECOPD groups. To control for multiple comparison errors, the BH method was applied for multiple testing correction of the P-values representing inter-group differences, followed by appropriate data visualization. Concurrently, Spearman correlation analysis was employed to clearly elucidate the correlations between the key genes and the significantly distinct immune cells.

Single-Cell RNA Sequencing Analysis

To investigate in depth the differential expression of key genes in lung tissue cells between STCOPD and AECOPD, we analyzed the single-cell dataset GSE268542 utilizing the Seurat (v 5.3.0) and SingleR (v 2.8.0) packages. Prior to the analysis, strict quality control was performed on the data to evaluate the impacts of doublet effects, the cell cycle, and batch effects. Subsequently, batch effect correction and gene normalization were executed utilizing the Harmony (v 1.2.3) and SCTransform algorithms. During the analysis, the resolution was set to 0.7, and the Uniform Manifold Approximation and Projection (UMAP) algorithm was employed to visualize the cell clusters. To identify the specific marker genes for each cell subpopulation, the default data matrix was set to the variance-stabilized SCT dataset, and the Wilcoxon rank-sum test was performed using the FindAllMarkers function within the Seurat package. We established stringent screening criteria to ensure the specificity of the marker genes: only significantly upregulated positive genes were extracted, requiring an adjusted P-value <0.05 following multiple testing correction via the BH method, a log-fold change threshold , and a difference in the expression proportion of candidate genes between the target cell population and the background cell population . Thereafter, the cell types were manually annotated and calibrated by referencing previous literature^35,36 and the CellMarker 2.0 database (http://117.50.127.228/CellMarker/CellMarkerBrowse.jsp). Then, heatmaps and stacked bar charts were utilized to present the cell annotation profiles and the proportional changes of cells across different COPD states, respectively. Finally, bubble plots and UMAP plots were employed to visualize the expression of key genes in distinct COPD states and their distribution within the cell clusters, respectively.

Additionally, the Monocle3 package (v 1.4.26) was further utilized to construct a pseudotime analysis model for immune cells, delineating their developmental trajectories and the expression levels of key genes during the developmental process. Concurrently, module clustering was performed on the developmental trajectories of the immune cells to highlight the biological functions of each module. Subsequently, we analyzed the level of intercellular communication utilizing the CellChat package (v 1.6.1) by calculating the total interaction strength of all ligand-receptor pairs within each signaling pathway, thereby inferring the mitochondrial functional demands among cells under different COPD states.

Molecular Docking

To infer the potential roles of the key genes in mitochondria based on the validated functions of existing compounds, molecular docking was performed. First, the potential associated compounds for the core genes were retrieved from the Drug Signatures Database (DSigDB) via the Enrichr website (https://maayanlab.cloud/Enrichr/) and ranked in ascending order according to their adjusted P-values. Subsequently, the 3D structures of the small molecule compounds were obtained from the PubChem database (https://pubchem.ncbi.nlm.nih.gov/). The protein structures of the core genes were acquired via the UniProt website (https://www.uniprot.org/); notably, only protein structure models determined by X-ray crystallography were selected for subsequent analyses. AutoDock Vina (v 1.5.7) was utilized to perform the molecular docking. During this process, the proteins. Plus, a web server (https://proteins.plus/) was employed to predict the active pockets of the macromolecular proteins, followed by the use of the Protein-Ligand Interaction Profiler (PLIP) tool (https://plip-tool.biotec.tu-dresden.de/plip-web/plip/index) to predict the interaction forces. Finally, PyMOL (v 2.6.0) and Discovery Studio 2021 were used to visualize the 3D and 2D results of the molecular docking.

Statistical Analysis

Except for the single-cell sequencing analysis, which was conducted using R software (version 4.4.3), all other analyses were performed using R software (version 4.5.0). Differences between groups were evaluated utilizing the Wilcoxon test. In analyses involving multiple hypothesis testing—including the screening of differentially expressed genes, enrichment analysis, comparison of immune infiltration differences, identification of single-cell marker genes, and Mendelian randomization causal inference—the Benjamini-Hochberg (BH) method was uniformly employed to calculate the false discovery rate (FDR) for multiple testing correction, unless otherwise specified. For conventional single-hypothesis tests, a P-value < 0.05 was considered statistically significant; for analyses subjected to multiple testing correction, an adjusted P-value (FDR) < 0.05 was defined as statistically significant.

Results

Initial Functional Enrichment of DEGs

In the differential analysis, we identified a total of 9137 DEGs (Figure 4A and B). Among them, 856 genes were upregulated, and 8281 genes were downregulated (Table S4). Notably, because the datasets were derived from parallel cohorts of the same research team and the disease groupings were completely mutually exclusive, conventional batch effect correction was not performed. Concurrently, constrained by public databases, the sample size of the training set could not be expanded with additional matched data, which introduced confounding factors and caused a certain degree of bias in the differential results. To identify the most prominent directions of biological changes from the aforementioned variations, we conducted GO and KEGG enrichment analyses. In the KEGG analysis, we found that the DEGs were highly associated with multiple neurodegenerative diseases, such as Alzheimer’s disease and Parkinson’s disease. They were also associated with muscle atrophy. Notably, these diseases are highly correlated with mitochondrial behavior. Furthermore, the KEGG analysis results indicated associations with viral and bacterial infections, elucidating the potential etiology of AECOPD (Figure 4C).

Figure 4 Preliminary screening and functional enrichment of differentially expressed mPTP-related genes in STCOPD and AECOPD: (A) Volcano plot. (B) Heatmap. (C) KEGG enrichment analysis of differentially expressed genes. (D) GO enrichment analysis of differentially expressed genes. (E) Venn diagram of the intersection between differentially expressed genes and mPTP-related genes. (F) PPI network diagram of differentially expressed mPTP-related genes.

In the GO enrichment analysis, a total of 1257 biological terms were obtained, comprising 849 BPs, 251 CCs, and 157 MFs. Subsequently, we discovered that these biological processes were closely related to mitochondrial behavior. The GO analysis suggested that, through the activation of ubiquitin-protein binding processes and gene translational and transcriptional activities, biological processes, including ubiquitin-proteasome system (UPS) metabolism, mitochondrial gene expression, macroautophagy, and vesicular transport, occurred within the organelle membrane, mitochondria, or nuclear regions (Figure 4D).

Identification of Core Genes

Intersecting the DEGs with the MPTPRGs yielded a total of 106 DE-MPTPRGs (Figure 4E). Subsequently, leveraging the PPI network provided by STRING, we identified 248 unique interaction pairs (Figure 4F). Within this PPI network, we observed that SPG7 exhibited high-intensity interactions with AFG3L2 and PPIF. This suggests that these three proteins might constitute essential components of the mPTP and regulate its transient opening function.

To further identify the DE-MPTPRGs causally associated with COPD hospitalization events, we conducted MR analysis. Following a preliminary MR analysis utilizing whole blood-derived eQTL data, we identified eight key genes—COQ2, DIABLO, IMMT, PARP1, SEC13, SLC25A3, SPG7, and TNF—that satisfied the core assumptions of the MR analysis. However, PARP1 and SLC25A3 were excluded as their regulatory trends were inconsistent with the findings of the differential analysis. Ultimately, we retained six key DE-MPTPRGs (details regarding their corresponding SNP counts and further information are presented in Table 1). However, an MR analysis utilizing the lung tissue-derived eQTL dataset did not yield additional or more meaningful results compared to those derived from whole blood. Details are provided in Table S5.

Table 1 Causal Relationships Between Differentially Expressed mPTP-Related Genes in STCOPD/AECOPD and COPD Hospitalization Events

During the MR analysis of the six key DE-MPTPRGs, after excluding SNPs significantly associated with COPD hospitalization events, 131 IVs were ultimately obtained. The F-statistic for each SNP exceeded 10, confirming their suitability for the MR analysis (Tables S6 and S7). Finally, employing the IVW method, we elucidated the significant causal associations between the key genes and COPD hospitalization events.

Forest plots and scatter plots demonstrated that DIABLO, IMMT, COQ2, SEC13, and TNF were negatively associated with COPD hospitalization events, whereas SPG7 exhibited a positive association (Figures S1 and S2). Funnel plots revealed no evident asymmetry, suggesting that this analysis adhered to the principle of random assortment dictated by Mendel’s Second Law (Figure S3). In the sensitivity analyses, both the heterogeneity test and the horizontal pleiotropy test yielded P-values greater than 0.05, indicating the absence of significant heterogeneity and confounding bias, thereby validating the robustness of the analytical results (Tables S8 and S9). Nevertheless, the potential for genetic pleiotropy may still exist. The LOO plots demonstrated that the exclusion of any single SNP did not significantly alter the overall effect of the remaining SNPs on any of the key genes, further supporting the reliability of the MR analysis results (Figure S4). Finally, we conducted a colocalization analysis to evaluate the association between the expression of the key protein SPG7 and COPD hospitalization events. Regrettably, we did not find sufficient evidence to support the presence of a significant shared causal variant between SPG7 expression and the occurrence of COPD hospitalization events . Further details are provided in Figure S5. We speculate that this outcome was influenced by the eQTL datasets employed, as neither the tissue- nor blood-derived eQTL data could comprehensively represent the pathological state of STCOPD. Secondly, the MR analysis results may imply that the exacerbation of COPD driven by SPG7 is not a local effect induced by a single potent mutation, but rather stems from the cumulative impact of multiple variants with minor effects. This is because, unlike colocalization analysis, MR has the capacity to aggregate multiple SNPs exhibiting consistent directions but weak individual effects, thereby capturing the accumulation of polygenic minor effects.

Validation and Functional Annotation of Key Genes

We utilized external validation sets to verify the diagnostic efficacy of the key genes, thereby further evaluating their critical roles in the transition from STCOPD to AECOPD. We integrated the key genes into a predictive nomogram and plotted calibration curves (Figure 5A–C). The calibration curves demonstrated a high degree of consistency between the predicted and actual observed outcomes, indicating acceptable predictive accuracy of the nomogram (Figure 5D–F). Furthermore, in the ROC curve evaluation, all areas under the curve (AUC) were greater than 0.7, further suggesting robust predictive performance (Figure 5G–I). To evaluate the clinical utility of the nomogram, decision curve analysis (DCA) was performed. The results revealed that the nomogram yielded a net benefit across a wide range of decision thresholds, indicating its superiority over the application of traditional single biomarkers in clinical scenarios (Figure S6).

Figure 5 Validation of key biomarkers: (A) Development of a nomogram utilizing the expression levels of each biomarker for the GSE148871 Sputum cohort. (B) Development of a nomogram utilizing the expression levels of each biomarker for the GSE148871 Blood cohort. (C) Development of a nomogram utilizing the expression levels of each biomarker for the GSE112165 Sputum cohort. (D) Calibration curve of the nomogram for the GSE148871 Sputum cohort. (E) Calibration curve of the nomogram for the GSE148871 Blood cohort. (F) Calibration curve of the nomogram for the GSE112165 Sputum cohort. (G) Receiver operating characteristic curve and area under the curve analysis to evaluate the predictive performance of the nomogram for the GSE148871 Sputum cohort. (H) Receiver operating characteristic curve and area under the curve analysis to evaluate the predictive performance of the nomogram for the GSE148871 Blood cohort. (I) Receiver operating characteristic curve and area under the curve analysis to evaluate the predictive performance of the nomogram for the GSE112165 Sputum cohort. ↑: Indicate specific numerical values on the scale.

Following the re-validation of the diagnostic efficacy, we evaluated the expression profiles of the key genes across the three validation sets. Boxplot analyses revealed that the expression levels of the six key genes did not all exhibit significant alteration trends between STCOPD and AECOPD, which may be attributed to the limited sample sizes of the validation sets. However, their overall median variation trends were fundamentally consistent with the findings from the training set. The primary discrepancy lay in the expression alterations of TNF, which were inconsistent across the three validation sets (Figure 6A–C). This phenomenon was even more pronounced in the correlation analysis (Figure 6E and F). We observed that across the three validation sets, the expressions of COQ2, DIABLO, IMMT, and SEC13 demonstrated positive synergistic correlations, whereas the expression of SPG7 was negatively correlated with all four aforementioned genes. As for TNF, its expression continued to exhibit significant disparities across the three validation sets.

Figure 6 Functional expression of Biomarkers: (A) Expression levels of each biomarker between STCOPD and AECOPD in the sputum samples of the GSE148871 validation set. (B) Expression levels of each biomarker between STCOPD and AECOPD in the blood samples of the GSE148871 validation set. (C) Expression levels of each biomarker between STCOPD and AECOPD in the sputum samples of the GSE112165 validation set. (D) Correlation analysis of each biomarker between STCOPD and AECOPD in the sputum samples of the GSE148871 validation set. (E) Correlation analysis of each biomarker between STCOPD and AECOPD in the blood samples of the GSE148871 validation set. (F) Correlation analysis of each biomarker between STCOPD and AECOPD in the sputum samples of the GSE112165 validation set. *P < 0.05, **P < 0.01, ***P < 0.001.

To comprehensively analyze the key molecular biological functions, GSEA was conducted. The GSEA results demonstrated that SPG7, SEC13, IMMT, DIABLO, COQ2, and TNF were significantly enriched in 108, 95, 86, 79, 77, and 76 molecular functions, respectively. Based on the expression trends of the key genes, we further discovered that during the transition from STCOPD to AECOPD, these molecular functions were predominantly concentrated in “voltage-gated monoatomic ion channel activity involved in regulation of presynaptic membrane potential”, “bitter taste receptor activity”, “serotonin receptor activity”, “neuropeptide receptor activity”, and “G protein-coupled amine receptor activity”. This implies a widespread activation of G protein-coupled receptors (GPCRs) (Figure 7A–F). The activation of these molecular functions further underscores the intimate relationship between the regulation of mitochondria-associated proteins and ion channel homeostasis, while also implying that the mPTP serves as a pivotal ion channel protein complex capable of executing mitochondrial quality control.

Figure 7 Gene set enrichment analysis of key biomarkers: (A) GSEA results for SPG7. (B) GSEA results for DIABLO. (C) GSEA results for SEC13. (D) GSEA results for IMMT. (E) GSEA results for COQ2. (F) GSEA results for TNF.

Immune Responses Under Changes in Key Genes

We evaluated the infiltration scores of 28 immune cell types in the STCOPD and AECOPD samples, visualizing the results utilizing stacked bar charts (Figure 8A). The analysis revealed significant differences in the infiltration abundances of 20 immune cell types between the two groups (P<0.05). Among these, only activated B cells and monocytes exhibited significantly increased infiltration in the AECOPD sputum samples (P<0.05). Conversely, the infiltration of the remaining 18 immune cell types—including T helper cells, natural killer (NK) cells, CD4⁺T cells, and CD8⁺T cells—was significantly decreased in the AECOPD sputum samples (P<0.05). This further reflects the state of immune activation present in the sputum of patients with AECOPD compared to those with STCOPD (Figure 8B). Furthermore, we performed a correlation analysis between the 20 differentially infiltrated immune cell types and the key genes. The results demonstrated (Figure 8C and D) that the infiltration of activated B cells and monocytes was significantly and positively correlated with SPG7. In stark contrast to the alterations of SPG7, IMMT and SEC13 exhibited a significant negative correlation with activated B cell infiltration, while showing a significant positive correlation with the infiltration of the majority of the downregulated immune cells. COQ2, DIABLO, and TNF demonstrated weak correlations with the infiltration of most immune cells, lacking representativeness. Such synchronous alterations between the immune infiltration patterns and the expression trends of key genes further suggest that the immune infiltration microenvironment might exert an influence on the mPTP through interactions with SPG7, IMMT, and SEC13.

Figure 8 Immune infiltration analysis results: (A) Stacked bar chart of single-sample GSEA scores for 28 immune cell types. (B) Differences in immune cell infiltration abundance between STCOPD and AECOPD. Red indicates upregulation in AECOPD, while blue indicates downregulation in AECOPD. (C) Correlation heatmap and network graph illustrating the internal correlations among the immune cells and their interaction network with the key genes. (D) Heatmap detailing the specific expression correlations between the key genes and the infiltrating immune cells. *P < 0.05, **P < 0.01, ***P < 0.001, ns: P > 0.05.

Single-Cell Analysis and Specific Expression of Key Genes

We further validated the expression alterations of the key DE-MPTPRGs between STCOPD and AECOPD at the cellular level through single-cell RNA sequencing (scRNA-seq) analysis. Following quality control, we obtained a total of 4730 cells and 19,194 genes (Figure S7A). Subsequently, we screened the top 3000 highly variable genes and utilized the top 41 principal components (PCs) for downstream analysis (Figure S7B and C). Through dimensionality reduction and clustering, 15 distinct cell clusters were generated (Figure S7D). By referencing the marker genes of pulmonary lineage cells from existing literature, we initially plotted the UMAP maps for 8 cell types and labeled their specific genes (Figure 9A and B). Similarly, for subsequent in-depth analysis, we further annotated 14 cell types (Figure 9C and D). However, the fifth cell cluster lacked specific cellular marker genes. To evaluate the type and state of this cell cluster, we discovered that within the fifth cell cluster (Table S10), MT-ND5 and MTRNR2L10 were activated, indicating enhanced mitochondrial respiratory activity in this cell population.³⁷ The elevated expression of FUS and HNRNPH1 implies that the cells within this cluster are undergoing transcriptome remodeling^38,39 in response to the altered inflammatory environment induced by IAV. Furthermore, the expression of STK4 was also upregulated, demonstrating that this cell cluster is highly associated with immune cells, including various T lymphocytes.⁴⁰ Therefore, we deduced that the fifth cell cluster primarily represents a mixed population of activated stress-state lymphocytes (ASL). In the stacked bar charts (Figure 9E and F), we observed that although the overall proportion of T lymphocytes remained virtually consistent, the proportion of ASL exhibited a multifold increase in AECOPD, which aligned with the variation trend of neutrophils. Additionally, we observed that the decline in the proportion of epithelial cells was predominantly driven by a multifold reduction in the proportion of alveolar type II epithelial cells.

Figure 9 Single-cell sequencing analysis results: (A) UMAP plot visualizing 8 distinct cell clusters identified through clustering analysis. (B) Heatmap displaying the cell annotation results and marker gene expression for the 8 distinct cell clusters. (C) UMAP plot visualizing 14 distinct cell clusters identified through clustering analysis. (D) Heatmap displaying the cell annotation results and marker gene expression for the 14 distinct cell clusters. (E) Stacked bar chart illustrating the proportional distribution of the 8 annotated cell types between the Stable COPD and AECOPD groups. (F) Stacked bar chart illustrating the proportional distribution of the 14 annotated cell types between the Stable COPD and AECOPD groups.

Subsequently, we investigated the expression profiles of the key genes across cell types between the groups. We observed that the expressions of COQ2, DIABLO, and TNF were not significant across the cells, whereas expression was predominantly observed for SPG7, IMMT, and SEC13 (Figure 10A), which could potentially be attributed to the limited sample size. The UMAP plots demonstrated the overall distribution and inter-group variation trends of the key gene expressions. The results indicated that although the inter-group variations in the expression of each key gene were not pronounced, the overall trends remained fundamentally consistent with the findings from the training set (Figure 10B–G). The expressions of COQ2, DIABLO, TNF, and SEC13 were relatively dispersed. However, SPG7 and IMMT were distinctly concentrated within the ASL cell cluster, implying that they might play crucial roles within this specific subpopulation.

Figure 10 Expression levels of biomarkers in each cell type: (A) Expression levels of each biomarker among various cells, AECOPD, and STCOPD. (B) UMAP plots showing the single-cell expression distribution of the biomarker SPG7 in the STCOPD and AECOPD groups. (C) UMAP plots showing the single-cell expression distribution of the biomarker DIABLO in the STCOPD and AECOPD groups. (D) UMAP plots showing the single-cell expression distribution of the biomarker SEC13 in the STCOPD and AECOPD groups. (E) UMAP plots showing the single-cell expression distribution of the biomarker IMMT in the STCOPD and AECOPD groups. (F) UMAP plots showing the single-cell expression distribution of the biomarker COQ2 in the STCOPD and AECOPD groups. (G) UMAP plots showing the single-cell expression distribution of the biomarker TNF in the STCOPD and AECOPD groups.

We observed that the key genes were primarily enriched in T lymphocytes. Therefore, we delineated the developmental trajectories of T lymphocytes utilizing pseudotime analysis. During T lymphocyte development, differentiation primarily bifurcates into two trajectories: one direction entails the differentiation of naive T cells into cytotoxic T cells, while the other involves differentiation into regulatory T cells, which subsequently progress to become the ASL subpopulation. Notably, compared with STCOPD, the abundance of the ASL cell cluster in AECOPD increased significantly, which potentially suggests a crucial role for ASL in immune activation and inflammatory responses (Figure 11A and B). Furthermore, Monocle3 clustered all T lymphocytes into a single partition (Figure 11C), further indicating that these cells belong to the same developmental lineage and are “closely related cells” with continuous differentiation relationships. Subsequently, we further analyzed the alterations in the expression levels of the key genes along the developmental trajectory. We found that, with the exception of TNF and SEC13, the expression levels of the remaining key genes were significantly elevated within the ASL cell cluster (Figure 11D), particularly SPG7. This further implies that SPG7 might be a pivotal gene regulating AECOPD via the mPTP. To determine the cellular functions within the ASL cluster, we constructed four co-expression gene modules (Figures 11E and S8). Given that ASL belongs to Module 4, we subsequently performed GO functional enrichment analysis on the genes within this module (Figure 11F). We discovered that the ASL cell cluster is highly associated with cellular respiration, oxidative phosphorylation, mitochondrial proton transmembrane transport, and mitochondrial membrane potential. This further suggests that the immuno-inflammatory state of AECOPD is influenced by the ASL cell cluster, which in turn is profoundly regulated by mitochondria. The underlying key mechanisms are likely highly correlated with the opening and closing of the mPTP.

Figure 11 Pseudotime analysis results of T cells: (A) UMAP plot displaying the clustering results of T cell subclusters. (B) Pseudotime trajectory analysis of T cells, comparing the developmental trajectories between the AECOPD and Stable COPD groups. (C) UMAP plot showing the partition of T cells used for the developmental trajectory construction. (D) Expression trends of the six key biomarkers mapped onto the T cell developmental trajectory. (E) UMAP plots illustrating the distinct gene expression modules along the T cell developmental trajectory. (F) Dot plot showing the gene enrichment analysis results associated with Module 4.

Finally, in the intercellular communication analysis, we observed that compared with STCOPD, the total number of interactions increased in AECOPD. This was predominantly characterized by an elevated number of interactions among immune cells, alongside a decline in interactions involving epithelial and endothelial cells (Figure 12A). Conversely, in stark contrast to the interaction count, the overall interaction strength among cells was significantly attenuated (Figure 12B). This indicates that during the AECOPD stage, the cellular architecture of the pulmonary tissue is disrupted, accompanied by the massive recruitment of immune cells; however, these recruited immune cells exist in a dysfunctional state. We postulate that this phenomenon might be attributed to the presence of dysfunctional mitochondria within these cells that cannot be effectively eliminated. This impairs ATP synthesis, thereby inducing oxidative stress, activating inflammatory pathways, and recruiting further immune cells, while simultaneously manifesting a state of diminished communication. This further corroborates the immuno-inflammatory infiltration microenvironment present in patients with AECOPD.

Figure 12 Analysis results of intercellular communication: (A) Number of interactions. (B) Interaction strength.

Regarding the overall findings of the single-cell analysis, constrained by an excessively small sample size, the expression of several key genes lacked statistical significance and exhibited unpronounced trends. Consequently, the single-cell analysis results serve merely as a supplementary reference to the preceding analyses.

Validation of Molecular Docking

The purpose of this study is not to screen for potential therapeutic agents for AECOPD, but rather to further validate whether the corresponding proteins of the key genes could serve as potential regulators of the mPTP, leveraging the confirmed functions of existing drugs. Regrettably, the database did not provide protein structure models for COQ2 and IMMT derived from X-ray crystallography; consequently, molecular docking for these two proteins was abandoned. Nevertheless, we obtained two potentially associated compounds for each core protein from the DSigDB database (Table 2). Finally, molecular docking was performed on the remaining four protein models, and the small molecule compounds bound to them all exhibited moderate to high binding affinities (Figure 13A–H).

Table 2 Association Between Core Biomarkers and Potential Compounds: Derived From DSigDB

Figure 13 Molecular docking results: (A) Molecular docking results showing the binding interaction between SPG7 and Oxidopamine. (B) Molecular docking results showing the binding interaction between SPG7 and Allococaine. (C) Molecular docking results showing the binding interaction between DIABLO and Cerivastatin. (D) Molecular docking results showing the binding interaction between DIABLO and Pterostilbene. (E) Molecular docking results showing the binding interaction between SEC13 and Aspirin. (F) Molecular docking results showing the binding interaction between SEC13 and Vitamin E. (G) Molecular docking results showing the binding interaction between TNF and Butein. (H) Molecular docking results showing the binding interaction between TNF and Ginsenoside Rh1.

From the presentation of all the aforementioned results, it is evident that the ultimate target gene of this study has been definitively identified. SPG7 demonstrated a crucial role across multiple key findings, which was further supported by the molecular docking results. Oxidopamine is a selective neurotoxin that induces dopaminergic neuronal damage via the generation of ROS and is widely utilized to establish animal models of Parkinson’s disease.⁴¹ The presence of potential binding sites between Oxidopamine and the SPG7 protein suggests that SPG7 might be a critical factor regulating mitochondrial function. Allococaine has been validated as a sodium channel inhibitor with the capacity to block voltage-gated sodium channels, thereby disrupting the stability of the membrane potential.⁴² Furthermore, the presence of even stronger binding sites between Allococaine and SPG7 indicates that SPG7 is highly correlated with membrane potential levels.

Discussion

COPD is an inflammatory disease characterized by progressive and irreversible structural and functional damage to lung tissue,⁴³ primarily manifesting as airflow obstruction, mucus hypersecretion, alveolar wall destruction, and airway smooth muscle cell proliferation. Among these structures, alveolar type II (AT2) epithelial cells play a critical role in maintaining the structural and functional integrity of lung tissue.⁴⁴ The extensive destruction and impairment of AT2 cells will lead to aberrant injury repair in COPD.⁴⁵ Therefore, maintaining a stable cellular functional state to protect lung tissue architecture and delay the decline in pulmonary function represents the core interventional objective in the treatment of COPD. However, for cells to exert their normal physiological functions, they predominantly rely on the energy supplied by mitochondria. Consequently, the cellular state is intimately associated with the functional regulation of mitochondria. Once mitochondrial regulatory mechanisms are disrupted, the most severe consequence is the mediated initiation of intrinsic cellular apoptosis.

Across different stages of COPD, this dysregulation of mitochondrial homeostasis serves as a key factor mediating lung tissue damage, among which mitophagy represents the most critical mitochondrial regulatory mechanism. The precise regulation of mitophagy has become a primary focus for numerous researchers, as well as one of the most controversial subjects in the field. Some studies suggest that cigarette smoke extract (CSE) exacerbates mitochondrial damage and depolarization by increasing mitophagy,⁴⁶ proposing that roflumilast exerts a protective effect on lung bronchial epithelial cells by inhibiting PINK1 expression.⁴⁷ Conversely, other studies^20,48,49 completely refute this perspective. These investigations argue that mitophagy attenuates the production of reactive oxygen species (ROS) in airway epithelial cells, thereby demonstrating that adequate mitophagy helps alleviate cigarette smoke-induced oxidative stress in lung tissue and mitigates COPD-associated airway inflammation. This contradiction is similarly evident in research concerning COPD-mediated cellular senescence. Several studies contend^50–52 that impaired mitophagy is a primary cause of cellular senescence, reflecting that insufficient mitophagy leads to the accumulation of damaged mitochondria. Alternatively, other research proposes that if autophagy is persistently activated and unregulated under stress conditions, it preferentially damages lung epithelial cells, driving pulmonary inflammation and injury in COPD by compromising epithelial cell function.⁵³ Furthermore, overactivation can lead to a decline in intracellular ATP levels, subsequently triggering programmed cell death or necrosis.^54–56 Regarding the contradictory phenomena across these studies, some scholars propose that the physiological demands or pathological alterations in autophagy levels might correspond to distinct COPD phenotypes.^17,57 We postulate that this phenomenon is associated not only with the COPD phenotype but also with the dynamic staging states of the disease.

Synthesizing the contradictory phenomena of the aforementioned previous studies and the preliminary results of the current study, we exploratively propose a mitochondria-driven hypothetical framework (Figure 2, left panel). The core of this hypothetical framework posits that “mitochondria strive to promote cell survival to maintain the structural and functional stability of pulmonary tissue, thereby protecting lung function and delaying pulmonary senescence”. We contend that for the structural cells of lung tissue, if cell death occurs, the physiological significance of all organelles is rendered obsolete, and pulmonary function will be further compromised. Therefore, this hypothesis initially proposes that upon exposure to massive external stimuli, cells preferentially initiate mitophagy; by sacrificing mitochondria, they attenuate intracellular energy metabolism to sustain cell survival. Mitophagy is contingent upon membrane potential depolarization, whereas the mitochondrial membrane potential is strictly governed by the ion transport processes of the mPTP.^16,22 Accordingly, the hypothesis further suggests that mitophagy is not the preferred choice for cellular autoregulation.²¹ Instead, cells are more inclined to maintain a high-energy state while sustaining survival, thereby executing complete cellular metabolic functions. Finally, the hypothesis proposes that “mitochondrial flickering” might act as the trigger for maintaining the stability of the mitochondrial membrane potential; specifically, enhancing the capacity for “flickering” can not only rescue cells on the verge of intrinsic apoptosis by augmenting autophagy, but also reverse non-essential autophagy and restore oxidative phosphorylation, thereby stabilizing cellular function.

Therefore, based on the aforementioned hypothetical framework, this study further deduces a novel COPD management strategy (Figure 2, right panel). We postulate that acute exacerbations of mild COPD are frequently accompanied by massive transient external stimuli. At this stage, reactive cellular autophagy is upregulated but remains insufficient to counteract the extensive damage, thereby predisposing cells to death. Consequently, by enhancing “mitochondrial flickering”, we can enable cells to maintain viability via mitophagy during the acute phase and autonomously restore cellular function following the resolution of the stimuli. Conversely, during the pulmonary senescence phase of STCOPD, we contend that chronic inflammatory stimuli drive mitophagy into a state of persistent overactivation, and the ensuing accumulation of senescent cells accelerates tissue senescence. At this juncture, we can similarly enhance “mitochondrial flickering” to avert tissue senescence induced by excessive mitophagy. Finally, in the end-stage of COPD, the organism exhibits a vicious cycle characterized by a low-energy state and hyperinflammation. At this point, solely enhancing “mitochondrial flickering” is ineffective. We should initially provide exogenous energy supplementation to increase mitophagy substrates. Subsequently, following anti-inflammatory interventions and the amelioration of the hyperinflammatory tissue microenvironment, enhancing “mitochondrial flickering” can prevent membrane potential collapse. This, in turn, activates mitophagy to sustain cell survival, ultimately allowing for the gradual restoration of cellular energy metabolism.

As the inaugural study to validate the aforementioned hypothetical framework, a pivotal task of this research is to verify the plausibility of this hypothesis. Through bioinformatic analyses, this study revealed that, compared with STCOPD, AECOPD is characterized by extensive AT2 cell death, accompanied by pronounced mitophagy and substantial energy metabolic demands, ultimately manifesting as widespread cellular dysfunction. Furthermore, we discovered that these alterations are highly correlated with mitochondrial membrane potential, ion channel activity, and transmembrane transport functions. This implies that the mPTP might play a crucial role in these processes. Subsequently, we identified six biomarkers—COQ2, DIABLO, IMMT, SEC13, SPG7, and TNF—which might establish causal relationships with AECOPD through the regulation of the mPTP. Finally, we postulate that SPG7 is a pivotal gene involved in the regulation of mitochondrial function in AECOPD by modulating the mPTP. Moreover, the SPG7 protein might interact with AFG3L2 or CypD to mediate protein modifications, thereby influencing the opening of the mPTP.

The Demand for Mitophagy is Prevalent in the AECOPD Stage

A minority of studies⁵⁸ propose that short-term exposure to PM2.5 or CSE increases susceptibility to AECOPD by elevating ROS production and subsequently activating mitophagy. These studies contend that excessive mitophagy, rather than its deficiency, contributes to the acute exacerbation of COPD. However, their arguments were established by directly targeting and inhibiting ROS. Consequently, we have reasonable grounds to question whether such alterations in mitophagy are not merely a survival requisite for cells to cope with oxidative stress. Encouragingly, a growing body of in vivo and in vitro research has validated a ubiquitous demand for mitophagy during the AECOPD stage. Initially, some scholars^3,48,59 postulated that mitophagy, initiated by CSE-induced dysregulation of mitochondria-associated protein homeostasis, exerts a protective role in AECOPD. Furthermore, an increasing number of studies^59–64 collectively assert that short-term exposure to CSE impairs mitophagy flux in pulmonary epithelial cells, which manifests as the accumulation of LC3-II, p62, and polyubiquitinated proteins. Conversely, enhancing mitophagy can alleviate this mitochondrial dysfunction.

The results of the present study further corroborate the plausibility of the hypothesis proposed in this section. We discovered that during the AECOPD stage, extensive AT2 cell death occurs, manifesting as a concurrent decline in both the number and strength of intercellular communications. This phenomenon is attributed to the fact that epithelial cells constitute the primary barrier between the organism and the external environment.⁶⁵ Infections triggered by external stimuli, such as CSE and IAV, primarily drive epithelial cells to experience earlier and more pronounced cellular energy depletion and mitochondrial dysfunction.⁶⁶ Consequently, to clear the infection and protect the cells, the organism manifests a robust energy demand. On the one hand, the tissues develop a hyperinflammatory microenvironment characterized by immune infiltration. This is evidenced by the profound activation of B cells in the sputum of patients with AECOPD. Initially, the activation of B cells requires energy; subsequently, activated B cells rely on oxidative phosphorylation to exert their immune functions.⁶⁷ On the other hand, the rapid activation of immune cells is observed within the pulmonary tissues of patients with AECOPD. This activation is predominantly concentrated within the ASL cell cluster, which concurrently exhibits a high demand for oxidative phosphorylation. Furthermore, this process is accompanied by massive neutrophil infiltration.

The energy requirements of immune cells and the maintenance of cell viability within infectious and hyperinflammatory microenvironments impose a tremendous energy burden on the organism, a phenomenon similarly reflected by mitochondrial alterations at the cellular level. The present study discovered that during the AECOPD stage, GPCRs are extensively activated, thereby providing an energy gain for the cells. However, the activation of these receptors is not entirely beneficial to mitochondria; rather, it likely maintains a state of equilibrium. Bitter taste receptors (TAS2Rs), a class of G protein-coupled receptors, have been demonstrated to possess anti-inflammatory and bronchodilatory activities.⁶⁸ Nevertheless, these effects are mediated through increased mitochondrial calcium uptake.⁶⁹ Consequently, TAS2Rs-mediated mitochondrial calcium overload leads to the loss of membrane potential, activates apoptotic pathways, and disrupts mitochondrial energy metabolism.⁷⁰ Although the activation of TAS2Rs elevates the risk of cell death, this risk is potentially mitigated and balanced by the concurrent activation of other GPCRs. For instance, the serotonin receptors identified in this study are functionally expressed within mitochondria, participating in the maintenance of mitochondrial ROS and calcium homeostasis, and enhancing ATP production efficiency to counteract stress and fluctuations in oxygen tension.⁷¹ This further corroborates the functional role of GPCRs localized to mitochondria.

The augmented energy production following stress is utilized for mitochondrial quality control to restore cellular function and promote cell survival. The present study reveals that in AECOPD samples, mitochondrial gene expression, macroautophagy, and vesicular transport are initiated, accompanied by the concurrent activation of the UPS. This implies the presence of abundant damaged proteins and mitochondria within the cells, alongside the widespread activation of mitophagy. Proteostasis is crucial for maintaining normal cellular processes and overall health.⁷² Fueled by adequate ATP,⁷³ the UPS can mediate the clearance of oxidatively damaged mitochondria under both aerobic and anaerobic conditions.⁷⁴ This process further accelerates the rate of mitophagy, thereby ensuring cell survival.

Nevertheless, the level of mitophagy remains insufficient. At the macroscopic level, it continues to exhibit similarities with mitochondrial dysfunction-related diseases, such as Alzheimer’s disease, Parkinson’s disease, and amyotrophic lateral sclerosis.⁷⁵ Concurrently, various intercellular communications also manifest as dysfunctional. Therefore, we contend that, at least for acute exacerbations in patients with non-end-stage COPD, enhancing mitophagy represents a priority strategy that aligns with the physiological demands of the organism.

SPG7 is a Potential Regulator of Mitochondrial Flickering

At the genetic level in AECOPD, the alterations in the six biomarkers further elucidate the cellular regulatory mechanisms following stress. The expression changes of these six genes can be categorized into three distinct profiles manifesting post-cellular stress. The first category is represented by SPG7 and DIABLO, which promote cell survival. SPG7 is postulated to regulate the transient opening of the mPTP to mitigate the collapse of oxidative phosphorylation, thereby preserving cell viability under continuous Ca²⁺ stress.^76,77 Although colocalization analysis does not support a significant shared causal variant between SPG7 and the occurrence of COPD hospitalization events, this does not negate the intrinsic impact of SPG7 on COPD-related hospitalizations. This is because COPD hospitalization events typically signify acute exacerbations or refractory severe COPD, wherein the confounding factors are not directly determined by a single gene or protein. If the hypothesis of the present study holds true, the protective upregulation of SPG7—acting in response to human energy demands and the decline in mitochondrial membrane potential—should manifest more prominently in patients with moderate-to-severe stable COPD. Furthermore, it should similarly exhibit systemic mitochondrial characteristics, thereby accumulating more extensively in whole blood. Only when the mitochondrial regulatory flux is thoroughly devastated, thereby inducing cell death and triggering an inflammatory burst, might this cumulative effect transition into AECOPD, culminating in a COPD hospitalization event. The downregulation of DIABLO, which encodes a mitochondria-released pro-apoptotic protein, indicates that cells are attempting to maximize survival through mitochondrial regulation.⁷⁸ Molecular docking further revealed that the DIABLO protein possesses binding sites for Cerivastatin and Pterostilbene. Cerivastatin can activate the mitochondrial apoptotic pathway,⁷⁹ whereas Pterostilbene can protect cells from damage.⁸⁰ This corroborates the biological significance of DIABLO downregulation. The second category involves detrimental changes manifesting as stress-induced injury, encompassing the downregulation of SEC13, IMMT, and COQ2. SEC13 is primarily involved in protein transport and inter-organelle communication; its downregulation under infection or oxidative stress attenuates the antiviral immune response and increases cellular susceptibility to oxidative damage.⁸¹ IMMT primarily maintains mitochondrial structural integrity and functional homeostasis. Its downregulation signifies that oxidative damage compromises mitochondrial function, subsequently impacting cell migration and survival.⁸² COQ2 primarily maintains mitochondrial energy metabolism by participating in coenzyme Q biosynthesis and exerts antioxidant effects. Its downregulation indicates that the mitochondrial oxidative phosphorylation pathway is impaired.⁸³ The final category is represented by TNF, reflecting the level of inflammation. The present study revealed inconsistent expression trends for TNF across different datasets. We postulate that this phenomenon might stem from the disparate inflammatory levels of COPD patients within the distinct datasets. Multiple studies have indicated that TNF-related protein levels are significantly elevated in patients with AECOPD.^84,85 Compared with non-COPD individuals, patients with COPD exhibit higher levels of TNF-α, IL-8, and hs-CRP. However, following therapeutic interventions, the levels of TNF-α, IL-8, and hs-CRP in the peripheral blood universally decrease.⁸⁶ Therefore, this phenomenon likely reflects the varying inflammatory states of COPD patients across the different datasets.

GSEA and single-cell sequencing analyses indicate that the alterations of these biomarkers correlate with the mitochondrial membrane potential via inner mitochondrial membrane ion channels, mitochondrial proton transmembrane transporters, and GPCRs, concurrently participating in the energy-yielding process of oxidative phosphorylation. Mitochondria maintain normal biological functions through numerous ion channels.¹⁶ Among these ion channels, the present study found that the regulation of Ca²⁺ channels plays a paramount role. Both TAS2Rs and serotonin receptors exert their biological functions by regulating the intramitochondrial Ca²⁺ concentration.^68,87 If Ca²⁺ uptake mediated by the mitochondrial calcium uniporter complex¹⁶ represents the maintenance of cellular ionic homeostasis, then the “mitochondrial flickering” mediated by the mPTP constitutes the autonomous homeostatic regulation of the mitochondria. This is because the opening of the mPTP does not merely precipitate membrane potential collapse;⁷⁶ an increasing body of evidence supports that mitochondrial flickering plays a vital role in mitochondrial Ca²⁺ homeostasis.⁷⁷ Consequently, we contend that mitochondrial quality control during the AECOPD stage is highly correlated with the mPTP.

Among the six key biomarkers, only the expression trend of SPG7 consistently exhibits a high correlation with the regions characterized by a demand for oxidative phosphorylation. Furthermore, in the PPI network, we discovered that SPG7 exhibits robust interactions with AFG3L2 and PPIF. CypD is a peptidyl-prolyl isomerase encoded by the PPIF gene, which is located in the mitochondrial matrix and plays a critical regulatory role in the mPTP.⁷⁶ To date, the exact molecular structure of the mPTP remains elusive, making CypD the sole identified critical regulatory factor of the mPTP.⁸⁸ This implies that SPG7 very likely regulates mPTP function via CypD.

We postulate that the mPTP is regulated via at least two distinct pathways. The first, conventional pathway involves mitochondrial Ca²⁺ overload or oxidative stress inducing CypD to trigger the prolonged opening of the mPTP, which ultimately culminates in cellular damage.⁷⁶ In the second pathway, we speculate that the SPG7 protein likely engages in a specific protein modification relationship with CypD. This interaction enables CypD to activate the transient opening function of the mPTP, thereby precluding its prolonged and detrimental activation.

Currently, an accumulating body of evidence substantiates our hypothesis. In 2015, the hypothesis postulating SPG7 as a core component of the mPTP was comprehensively proposed; this study⁸⁹ indicated that the depletion of SPG7 leads to higher mitochondrial Ca²⁺ retention and exhibits characteristics similar to CypD knockdown, namely preventing Ca²⁺- and ROS-induced ΔΨm depolarization and cell death. However, the overexpression of SPG7 did not increase mitochondrial sensitivity to Ca²⁺-induced mPTP opening in the control group. This suggests that SPG7-mediated mPTP opening might represent a regulatory mechanism distinct from its full and sustained opening. In January 2020, the hypothesis of SPG7 acting as a core mPTP component was refuted; a subsequent study⁸⁸ reported that no Ca²⁺-induced mPTP response was observed following SPG7 knockdown. Nevertheless, the study discovered that ATP levels increased significantly upon SPG7 overexpression. This further implies the stabilizing effect of SPG7 on the mitochondrial membrane potential. In October of the same year, another study⁷⁷ directly confirmed that the functional expression of SPG7 is a prerequisite for CypD-mediated⁷⁶ transient opening of the mPTP, and verified for the first time that SPG7 deficiency leads to neurodegeneration due to impaired mitochondrial flickering. In 2022, researchers⁹⁰ further proposed that when AFG3L2 and SPG7 form a heterohexamer, the mPTP opens abnormally and induces cellular apoptosis. This m-AAA complex executes protein quality control within the inner membrane by selectively degrading unassembled and damaged proteins, while concurrently participating in the assembly of respiratory chain complexes.⁷⁷ Therefore, the opening of the mPTP at this juncture is highly likely attributed to SPG7 losing its regulatory capacity over the mPTP following the formation of the heterohexamer. Accordingly, we can deduce that SPG7 is likely a pivotal factor regulating mitochondrial flickering, functioning independently of the traditional pathway of persistent mPTP activation.

Mitochondrial flickering serves as an efficacious pathway facilitating rapid Ca²⁺ release from mitochondria; furthermore, it can scavenge mtROS, activate the respiratory chain, and subsequently promote ATP synthesis.^77,91 This aligns perfectly with the observed upregulation of SPG7. Consequently, we postulate that the protective effects elicited by SPG7 knockdown or CypD knockdown are merely transient. Once intramitochondrial oxidative stress and calcium overload completely overwhelm the autophagy flux, cell death remains inevitable—it is merely a matter of time. Only by restoring or augmenting the intrinsic ion regulatory mechanisms of mitochondria can cell survival be durably sustained.

Limitations of the Study

Although this study proposes several enlightening scientific hypotheses, as an exploratory investigation, it inevitably possesses certain limitations. First, constrained by the available resources in public databases, the sample size of this study is relatively limited (particularly the small scale of the training set and the single-cell transcriptome dataset), which, to some extent, reduces statistical power and may increase the risk of sampling error. Second, racial and population heterogeneity exists among the different datasets included in the analysis. Concurrently, eQTL data derived from healthy populations may not fully map the true pathological states and phenotypic characteristics of patients with STCOPD. Furthermore, microenvironmental discrepancies among different biological specimen sources (eg, peripheral blood, lung tissues, induced sputum), coupled with the diversity of triggering factors for AECOPD, may all act as potential confounding factors influencing the robustness of the analytical results. Third, the colocalization analysis failed to confirm the existence of a significant shared causal variant between SPG7 and COPD hospitalization events. Therefore, the specific biological functions of SPG7 as a potential key regulatory factor still urgently require validation through future in vivo and in vitro experiments as well as prospective cohort studies. Finally, during the processes of data mining and statistical fitting, multiple hypothesis testing and potential unmeasured confounding factors might elevate the risk of false-positive results.

Conclusion

In summary, this study exploratively proposes a “mitochondria-mediated cell survival hypothesis” and preliminarily corroborates its plausibility through bioinformatic analyses. Specifically, our research reveals widespread demands for energy and mitophagy during the AECOPD stage. Cells may satisfy this physiological demand and maintain viability by enhancing “mitochondrial flickering”, thereby preserving the structural integrity of the tissues. Additionally, we identified six biomarkers causally associated with acute exacerbations of COPD, among which SPG7 exhibited a profound correlation with mitochondrial flickering. However, the non-significant colocalization analysis attenuates the robustness of this hypothesis. Although the integrated results of this study suggest that SPG7 might be a potential key protein regulating the mPTP, the direct associations among these speculations, underlying mechanisms, and clinical evidence remain insufficiently explored. Given that persistent pulmonary senescence, frequent acute exacerbations, and multiple pulmonary infections alongside other related complications constitute the primary survival burden for patients with COPD, validating the mitochondria-related hypotheses of this study and integrating mitochondrial dysfunction with clinical outcome data represent a vital direction for future research endeavors.