An integrated investigation of mitochondrial genes in COPD reveals the causal effect of NDUFS2 by regulating pulmonary macrophages

Abstract

Background

Despite the increasing body of evidence that mitochondrial activities implicate in chronic obstructive pulmonary disease (COPD), we are still far from a causal-logical and mechanistic understanding of the mitochondrial malfunctions in COPD pathogenesis.

Results

Differential expression genes (DEGs) from six publicly available bulk human lung tissue transcriptomic datasets of COPD patients were intersected with the known mitochondria-related genes from MitoCarta3.0 to obtain mitochondria-related DEGs associated with COPD (MitoDEGs). The 32 hub MitoDEGs identified from protein-protein interaction (PPI) networks demonstrated superior overall diagnostic efficacy to non-hub MitoDEGs. Random forest (RF) analysis, least absolute shrinkage and selection operator (LASSO) regression, and Mendelian Randomization (MR) analysis of hub MitoDEGs further nominated NDUFS2, CAT, and MRPL2 as causal MitoDEGs for COPD, whose predominate expressions in pulmonary macrophages were revealed by an independent single-cell transcriptomic dataset of COPD human lungs. Finally, NDUFS2 was evaluated as the top-ranked contributor to COPD in the nomogram model and its downregulation in pulmonary macrophages could result in pro-inflammatory secretion, enhanced intercellular communications, whereas depressed phagocytosis of macrophages as revealed by gene set variation analysis (GSVA) and cell-cell interaction (CCI) analysis of single-cell transcriptomic dataset of COPD human lungs, which was later confirmed in COPD mouse model and macrophage cell lines.

Conclusions

Our study established the causal linkage between mitochondrial malfunctions and COPD, providing a potential therapeutic avenue to alleviate pulmonary inflammation accounting for COPD by targeting mitochondria-related genes. NDUFS2, a canonical component of mitochondrial electron respiratory chain, was highlighted instrumental for the susceptibility of risk-exposed individuals to COPD.

Background

As a heterogeneous lung condition characterized by chronic respiratory symptoms due to abnormalities of the airway and/or alveoli that cause persistent and progressive airflow obstruction [1], chronic obstructive pulmonary disease (COPD) is a leading cause of morbidity and mortality worldwide, imposing an enormous burden on individuals, families, and society [2]. Over the decades, our understanding of the pathogenesis of COPD has changed dramatically from a self-inflicted disease induced by cigarette smoking mainly occurring in the elderly to an outcome of many dynamic and cumulative gene-environment interactions that occur through the lifetime of the individual [3]. To date, the best documented genetic risk factor for COPD, albeit rare, are SERPINA1 mutations accounting for α-1 antitrypsin (AATD) deficiency [4]. Numerous other genetic variants associated with COPD have been identified [5]. However, their individual effect size is small with their causality of COPD largely uncertain [6].

While forming a dynamic and interconnected network intimately integrated with other cellular compartments to regulate multiple biological processes, including cellular metabolism, signal transduction, proliferation, migration, autophagy, and apoptosis, mitochondria also operate beyond the boundaries of the cell to regulate communication between cells and tissues, making them a critical signaling hub in COPD pathogenesis that bridges inflammation, bioenergetics, metabolism, ion homeostasis, and oxidative stress, as well as a potential therapeutic target for COPD [7]. Despite the revolutionary progress in mitochondrial biology, with the emergence of high-confidence database of mitochondrial proteins, such as MitoCarta3.0 [8] and MitoCoP [9], our understanding of mitochondria in COPD is still left behind [10]. On one hand, previous studies of mitochondrial dysfunction were mainly focused on epithelial cells due to their direct exposure to diverse COPD risk factors. Human lung tissue consists of over 40 cell types with different functions as well as energy demands and therefore different number as well as biology of mitochondria [11]. Mitochondrial dynamics of other cell types in the context of COPD, together with their potential influence on each other and their contribution to COPD development warrant further exploration. On the other hand, considering the remarkably sensitivity, responsiveness, and plasticity of mitochondria to various stressors [7, 12], it is currently impossible to distinguish whether the alterations in mitochondria observed in COPD lung tissues are merely adaptive responses or major factors that cause, promote, or aggravate COPD.

Nowadays, transcriptomic studies have uncovered enormous amounts of differentially expressed genes (DEGs) between healthy individuals and COPD patients [13, 14]. Functional annotations of these DEGs identified the biological processes significantly perturbed under COPD conditions, yet were insufficient to distinguish DEGs directly responsible for COPD from those merely markers of other causal genes. In recent years, Mendelian randomization (MR) stands out among analytic approaches used in association studies [15]. Following Mendel’s laws of inheritance, MR uses genetic variations strongly associated with exposure factors as instrumental variables (IVs) to assess the causality of an observed association between a modifiable exposure or risk factor and a clinically relevant outcome [16]. When combined with expression quantitative trait loci (eQTL) data, MR also does well inferring causal relationships between genes expression and disease outcomes [17].

In this study, by leveraging MitoCarta3.0 database, we identified mitochondria-related DEGs associated with COPD (MitoDEGs) from six bulk human lung tissue transcriptomic datasets of COPD patients. Through integrated analysis of these MitoDEGs, including machine learning and MR analysis, we found that NDUFS2, MRPL2, and CAT, predominately expressed in pulmonary macrophages, were causally involved in COPD. Functional analysis together with experimental investigation revealed the modulatory role of NDUFS2 in pulmonary macrophages, providing mechanistic linkage of NDUFS2 dysregulation with COPD pathogenesis.

Methods

Ethics statement

All animal experiments were approved by the Institutional Animal Care and Use Committee at Nanjing Medical University (approval number SYXK-2023-0029) and conformed to the ARRIVE guidelines. All datasets and summarized statistics utilized in the MR analyses were generated by previous studies, for which ethical approval and individual consent were obtained for all original studies.

Study design

We collected bulk human lung tissue transcriptomic datasets of COPD patients and control individuals from public repository NCBI GEO and carried out differential gene expression analysis to identify COPD-associated DEGs of every dataset, which were further intersected with mitochondria-related gene lists retrieved from the MitoCarta3.0 database to pick out MitoDEGs. Protein-protein interaction (PPI) networks of MitoDEGs were constructed for exploration of hub MitoDEGs, which were subsequently subjected to machine learning and MR analysis to screen causal MitoDEGs for COPD. After the cellular expression distribution of every causal MitoDEGs within the lung was demonstrated in single-cell RNA sequencing (scRNA-seq) dataset of COPD human lung tissues, a nomogram model was constructed to evaluate the contribution of their dysregulation to risk for COPD. Functional analysis was performed to infer the influence of the causal MitoDEG with the greatest contribution to COPD on biological processes, signaling pathways, and intercellular communications, which was preliminarily verified in COPD mouse model and in vitro experiments. The overall flowchart is shown in Fig. 1.

Fig. 1

Flowchart of the study. Six publicly availabe bulk human lung tissue transcriptomic datasets including 321 COPD patients and 164 controls individuals were collected and subjected to integrated analysis to screen for causal mitochondria-related genes for COPD and infer their functions underlying COPD, followed by experimental verification

Data resource and quality control

Six COPD lung tissue datasets (accession number: GSE103174, GSE106986, GSE151052, GSE38974, GSE47460, GSE57148) were obtained from public repository NCBI GEO (https://www.ncbi.nlm.nih.gov/) [18].Datasets were accessed by using GEO query package (version 2.68.0). Principal component analysis (PCA) was used to identify outlier samples of each dataset (Supplementary Figure S1a). Specific sample information for each dataset was provided in Supplementary Table S1.

Identification of DEGs and functional enrichment analysis

DEGs between COPD and normal groups were identified using the limma (version 3.56.2) package, with identification criteria of p < 0.05 and | logFC (Fold-change) |≥ 0.585), corresponding to a FC of 1.5 that is commonly used in transcriptomic studies to identify DEGs [19].Volcano and heatmap for visualization of DEGs via ggVolcano (version 0.0.2) and ComplexHeatmap (version 2.16.0) packages, which could better demonstrate the level of difference in DEGs between the COPD and normal groups. Gene set enrichment analysis (GSEA), Kyoto encyclopedia of genes and genomes (KEGG) pathway and gene ontology (GO) enrichment were performed to explore the potential biological properties of DEGs using the clusterProfiler (version 4.8.1) package, and the results were visualized by the GseaVis (version 0.0.5) and ggplot2 (version 3.4.1) packages.

PPI network construction and identification of hub MitoDEGs

MitoDEGs were obtained by intersecting DEGs from each dataset with 1,136 mitochondria-associated genes from MitoCarta3.0 (http://www.broadinstitute.org/mitocarta) [8] and visualized by the VennDiagram (version 1.7.3) package generated venn diagram. After all the MitoDEGs were imported into STRING (https://string-db.org/) to obtain the PPI networks and the results were uploaded into Cytoscape software, the modules were grouped by MCODE (criteria of degree cut-off = 2, node score cut-off = 0.2, k-core = 2, and max depth = 100) and those with high scores in each group were taken as the identified hub MitoDEGs [20].

Receiver operating characteristic (ROC) analysis

The ability of hub MitoDEGs and non-hub MitoDEGs to predict COPD was examined and compared using the pROC (version 1.18.5) package. Specifically, logistic regression of the two models was first established separately by the glm () algorithm, the curves were established with the false positive rate represented by 1-Specificity as the x-axis and the sensitivity representing the true positive rate as the y-axis, and an area under the curve (AUC) greater than 0.7 was considered to have good predictive ability. The roc.test () function was then used to compare the strength of the two models in their ability to predict COPD.

Machine learning

Two machine learning algorithms, random forest (RF) and least absolute shrinkage and selection operator (LASSO) regression analyses, were used to further screen for signature genes from the 32 hub MitoDEGs. First, the six datasets were merged and normalized using the sva (version 3.50.0) package. Then, 472 samples were divided into a training set (330 samples) and a test set (142 samples) in a ratio of 7:3.

The RF model establishes the optimal number of variables by calculating the average error rate of candidate pivotal genes. The importance of the genes was ranked using the randomForest (version 4.7.1) package and the genes with MeanDecreaseAccuracy greater than 5 and MeanDecreaseGini greater than 5 were selected. To validate the RF model, we constructed ROC curves for internal validation using the GSE57148 dataset, which contains the largest number of samples, and for external validation using the GSE76925 dataset, which includes 110 COPD samples and 40 control samples, obtained from GEO [21]. LASSO represents a regression analysis algorithm that applies regularization for variable selection. Next, we performed LASSO regression analyses of the RF screened genes using the glmnet (version 4.1.8) package to further identify the characterized genes.

MR analysis

To assess whether the 32 hub MitoDEGs were causally associated with COPD, we performed MR analyses using eQTL data of the 32 genes (exposure) with COPD (outcome). First, the eQTL data for each gene were obtained from eQTLGen (https://eqtlgen.org/), which collects genetic data on gene expression of 31,684 individuals from 37 datasets derived from blood samples from mostly healthy European individuals [22]. We downloaded the updated statistically significant blood circ-eQTL data for 2019 on 10 August 2023. To confirm the genetic variant reliably associated with risk factor, the cis-eQTL with p<5×10 ^− 8, minor allele frequency (MAF) > 0.01, and F value [calculate as (beta/SE)²] > 10 [23] were selected as IVs for analysis and performed clumping by setting a pairwise linkage disequilibrium (LD) cutoff value of r² < 0.01 and clump_kb=10000 [24]. Phenoscanner V2 was utilized to eliminate single nucleotide polymorphisms (SNPs) with potential associations with confounding variables or outcomes [25]. The outcome GWAS data of COPD (“ukb-d-COPD_EARLYANDLATER”) (https://gwas.mrcieu.ac.uk/), comprising 1897 cases and 359,297 controls, were sourced from the IEU OpenGWAS database on 20 August 2023. After identifying the causal MitoDEGs, we performed an additional MR analysis to validate the relationship between the three causal MitoDEGs and lung function (FEV1/FVC) using the GWAS dataset from ebi-a-GCST007431, which includes 321,047 samples and was sourced from the IEU OpenGWAS database. We harmonized the exposure-outcome datasets (i.e., consistent direction of effect between exposure and outcome associations) and performed two-sample MR analyses using the TwoSampleMR (version 0.5.7) package [26] and sensitivity analyses using a variety of analytical methods. In MR results, Wald Ratio method was selected for assessment if there was only one SNP, and inverse variance weighted (IVW) method was selected as the primary assessment approach if there were multiple SNPs. For MR results with statistical significance (p<0.05), further sensitivity analyses were conducted, and results that did not meet the sensitivity analysis criteria were excluded. Specifically, heterogeneity between individual estimates of genetic variance was assessed using the MR egger and IVW methods of Cochran’s Q test [27]. Horizontal pleiotropy of IVs was tested using the MR-Egger intercept method, where the intercept indicates the average pleiotropic effect between genetic variants (the average direct effect of the variants on the outcome), and evidence of directional pleiotropy exists if the intercept is different from zero (MR-Egger test) [28]. The leave-one-out sensitivity test was used to check whether the results were caused by any single SNP, and to generate forest plots.

scRNA-seq data processing

The COPD scRNA-seq dataset as we previously reported (accession number GSE171541) [29] was analyzed to characterize cellular localization by visualizing the uniform manifold approximation and projection (UMAP) information of genes using the “FeaturePlot” function in the Seurat (version 5.0.3) package. To better exclude the effect of smoking, we reorganized the smoking history and disease diagnosis. Briefly, all information was classified into four groups based on whether they smoked and whether they had COPD, and the expression of genes under different groupings was visualized using the ggplot2 package as mentioned above.

Gene expression correlation analysis

To estimate gene expression correlations in macrophages across individuals in scRNA-seq dataset GSE171541, we calculated the average expression levels of each gene in macrophages in each individual using the AverageExpression () function and then performed Pearson correlation analysis between average expression values of one gene of 9 individuals in macrophages and average expression values of another gene of 9 individuals in macrophages. Highly correlated genes are defined as those with Pearson correlation coefficient ≥ 0.8 or ≤ − 0.8. The top 50 genes with the highest correlation coefficient were visualized using the ggplot () function.

Nomogram construction

Alveolar macrophage (AM) bulk RNA sequence (bulk RNA-seq) dataset was also obtained from NCBI GEO (accession number GSE130928), containing 22 COPD samples and 66 Control samples [30]. First, we showed the gene expression using the ggplot2 (version 3.4.1) package as mentioned above. Then, a nomogram predicting the risk of COPD based on the expression of CAT, MRPL2, and NDUFS2 in alveolar macrophages was constructed using logistic regression with the rms (version 6.7.1) package. A bootstrap with 1000 repetitions using the calibrate function was used for bias correction to assess the validity of the nomogram.

Gene set variation analysis (GSVA) and GSEA based on scRNA-seq data

For further exploring the function of NDUFS2 in macrophages, we analyzed it based on scRNA-seq data. Specifically, we used the "subset" function to propose macrophages for further analysis. Then, macrophages were divided into NDUFS2 low expressed (macrophages ^NDUFS2low) and NDUFS2 high expressed (macrophages ^NDUFS2high) groups based on an NDUFS2 mean value of 0.5940285. Next, we compared the differences in pathway activity in the macrophages ^NDUFS2low group relative to the macrophages ^NDUFS2high group by GSVA (version 1.48.3), and the differences in biological function by GSEA.

CellChat analysis

To infer the potential impact of NDUFS2 expression levels on intercellular communication in lung macrophages, we performed CellChat analyses between the macrophages ^NDUFS2high and macrophages ^NDUFS2low. Specifically, the overall sample grouping information was redefined based on the NDUFS2 expression in macrophages, dividing the samples into macrophages ^NDUFS2high and macrophages ^NDUFS2low groups. The “netVisual_diffInteraction” function was used to visualize the signaling differences incoming to and outgoing from macrophages between the two groups, and the “netVisual_heatmap” function was used to visualize the overall communication differences between cells. The visualization of all signaling pathways was achieved using the “rankNet” function. Considering the downregulation of NDUFS2 expression in COPD macrophages, comparing communication analyses based on disease context would allow us to better explore the impact of NDUFS2 on cell communication. Therefore, we compared the communication differences between macrophages with high NDUFS2 expression in the control group (macrophages ^NDUFS2high in Ctrl) and macrophages with low NDUFS2 expression in the COPD group (macrophages ^NDUFS2high in COPD). All visualizations were achieved using the above mentioned functions.

Mice and cigarette smoke (CS) exposure

Twelve 8-week-old male C57BL/6J mice were purchased from the Model Animal Research Center of Nanjing University, National Resource Center for Mutant Mice, and randomly divided into two groups (six mice per group) for exposure to either CS or room air (Air) for 6 months. Huangshan brand cigarettes (China Tobacco Anhui Industrial Co., Ltd., Bengbu, Anhui, China) were used in a CSM-100C inhalation exposure system (TOW-Int Tech, Shanghai, China) for CS exposure, with total particulate matter (TPM) monitored in real-time using a Microdust Pro real-time dust monitor (Casella, Germany). Based on previous studies, the CS exposure environment was set with an average TPM of 500 mg/m³ and a carbon monoxide (CO) concentration of 300 parts per million (ppm), with exposure occurring twice daily for one hour each session, over a period of 6 months. At the end of the exposure period, lung tissues from the mice were collected for pathological analysis to confirm the successful establishment of the COPD model in mice.

Isolation of primary mouse AMs

On the final day of the 6-month exposure regimen, after completing the exposure, we isolated primary AMs from bronchoalveolar lavage fluids (BALFs) as previously described [31]. Specifically, mice were anesthetized with avertin (Sigma Aldrich, St. Louis, MO, USA) and dissected to expose the trachea, which was then cannulated with a 26 G syringe needle (BD, Franklin Lakes, NJ, USA). The lungs were lavaged with a total of 15 ml of ice-cold PBS, and gently massaged to ensure thorough recovery, in 20 aliquots of approximately 0.75 ml each. The BALFs were then centrifuged at 300 g for 15min at 4°C, and the cell pellets were resuspended in RPMI 1640 medium (Gibco, Burlington, ON, USA) containing 100 μg/ml streptomycin, 100 U/ml penicillin, 10 mmol/l HEPES, and 50 μmol/l β-mercaptoethanol (Sigma Aldrich). The cells were incubated at 37°C in a 5% CO₂ atmosphere for one hour. Subsequently, non-adherent cells were carefully washed away with sterile PBS, leaving the adherent, purified primary mouse AMs.

Histopathological analysis

The pathological characteristics of COPD were determined through hematoxylin and eosin (H&E) and Masson staining. After lung lobes were extracted from 4% paraformaldehyde fixation, they were dehydrated in 30% sucrose, followed by embedding in optimal cutting temperature (OCT) compound. Subsequently, according to the manufacturer's instructions, H&E and Masson staining were performed to assess the pathological features of COPD. Subsequently, we quantified the pathological changes using mean linear intercept (MLI) and destructive index (DI) as previously described [32]. Briefly, a crosshair grid was drawn at the center of each field of view, and the total length of each line divided by the number of alveolar intercepts yielded the MLI, which represents the average distance between alveolar surfaces. Supplementarily, a grid with 42 points was superimposed on the lung fields to record the number of normal (N) and destroyed (D) alveoli. The DI was calculated using the formula: DI = D / (D + N) × 100.

Immunofluorescence staining

Immunofluorescence staining on frozen sections was performed to determine both the differential gene expression and cellular localization, allowing for the identification of gene expression variances and their cellular distribution. For immunofluorescence staining, primary antibodies targeting NDUFS2 (Boster Biological Technology, Cat#A05618-3), F4/80 (Abcam Cat# ab204266), and HSP60 (Proteintech Cat# CL647-15282) were used. Secondary antibodies conjugated with fluorescent markers corresponding to the respective attributes were then employed for incubation. Finally, fluorescence images were collected using a laser scanning microscope and analyzed using ImageJ (National Institutes of Health, Bethesda, MD, USA).

Cell culture and transfection of small interfering RNAs (siRNAs)

The immortalized human-derived monocyte cell line THP-1 and the mouse-derived macrophage cell line RAW264.7 were cultured in 1640 medium containing 10% fetal bovine serum supplemented with 100 μg/ml streptomycin, 100 U/ml penicillin, 0.25μg/ml amphotericin B, and 10 mmol/l HEPES and incubated at 37°C in a humidified atmosphere with 5% CO₂. THP-1 cells were differentiated into macrophages by incubation with phorbol 12-myristate 13-acetate (PMA) (MCE, Cat#HY-18739) as previously described [33]. To silence NDUFS2 in macrophages, siRNAs against NDUFS2 were designed and synthesized by Ribobio (Guangzhou, Guangdong, China). For transient transfection of siRNAs, Lipofectamine 3000 reagent (Invitrogen, Grand Island, NY, USA) was mixed with siRNAs at a ratio of 12 μl: 10 μl 20 μM siRNAs storage solutions and Lipofectamine 3000 reagent according to the manufacturer’s protocol. Target sequences of all siRNAs against NDUFS2 used in this study were provided in Supplementary Table S2.

RNA extraction and quantitative real-time PCR (qPCR)

Total RNA was extracted using Trizol reagent (Sigma, Cat#T9424). Then, cDNA was obtained by reverse transcription of RNA (Vazyme, Cat#R211-01). Finally, qPCR was conducted using SYBR Green dye (Vazyme, Cat#Q341-03). The relative expression levels of the target genes were calculated using the Ct method (2^-ΔΔCt), normalized using the reference gene RPLP0. The primer sequences are provided in Supplementary Table S3.

Western blot analysis

Protein samples were collected, and the concentration was determined using BCA (ThermoFisher, Cat#23225). Proteins were separated using 10% SDS-PAGE and transferred to PVDF membrane. The membrane was incubated with the primary antibody overnight at 4°C. The following primary antibodies were used: anti-NDUFS2 (Boster Biological Technology, Cat#A05618-3), anti-GAPDH (Proteintech Cat# CL650-60004). The PVDF membranes were incubated with secondary antibodies for 1h at room temperature and then subjected to assay using the ECL system.

Phagocytic function of macrophages

The phagocytic function of macrophages was measured using the Cell Meter Fluorescence Phagocytosis Assay Kit (AAT Bioquest, Cat#21225) according to the manufacturer's instructions. Texas Red signal indicated the endocytosis of Protonex 600 Red-Latex beads within macrophages, while CytoTrace Green was used to assess cell numbers as an internal control [34].

Statistical analysis

Unless otherwise stated, each experiment was carried out in triplicate. All analyses were run in R (version 4.2.2) and R Studio (version 1.0.143). The data were expressed as mean ± standard deviation. Student’s t-test was used to compare the two groups with GraphPad Prism software, considering p < 0.05 as the threshold for significance. The legend contained specific statistical information.

Results

Functional enrichment analysis of COPD-associated DEGs emphasizes the mitochondrial alterrations

We first obtained 6 bulk RNA-seq datasets of COPD lung tissues (GSE103174, GSE106986, GSE151052, GSE38974, GSE47460, GSE57148) from GEO database, including 321 COPD samples and 164 normal samples altogether (Supplementary Table S1). After eliminating outlier samples by PCA (Supplementary Figure S1a), we conducted differential gene expression analysis of every dataset (Fig. 2a and Supplementary Figure S1b). To independently assess molecular alterations associated with COPD in these datasets, we carried out functional enrichment analysis of these COPD-associated DEGs separately. As expected, biological ontologies determined by GSEA (Fig. 2b) and KEGG pathway analysis (Supplementary Figure S1c) appeared to be heterogeneous among these datasets, which could be ascribed to the heterogeneity of COPD, yet they still highlighted inflammation, tissue fibrosis, and metabolic remodeling as the general status of COPD lung tissues. Notably, biological ontologies associated with mitochondria (such as mitochondrial inner membrane, mitochondrial large ribosomal, mitochondrial respiratory complex I/IV, regulation of mitochondrion organization, mitochondrial ATP synthesis coupled electron transport, release of cytochrome C from mitochondria, mitochondrial electron transport, and NADH to ubiquinone) were among the top dysregulated GO terms associated with COPD in all the datasets (Fig. 2c), suggesting the involvement of mitochondria in COPD pathogenesis.

Fig. 2

Functional enrichment analysis of COPD-associated DEGs emphasizes the mitochondrial alternations. a Volcano plots for COPD-associated DEGs in six bulk RNA-seq datasets of human lung tissues from COPD patients and control individuals with the number of up- and down-regulated DEGs of each dataset displayed in the plot. One dot represents a gene. Blue dots represent down-regulated DEGs. Red dots represent up-regulated DEGs. Gray dots represent genes showed no significant differential expression. b GSEA analysis of COPD-associated DEGs in six bulk RNA-seq datasets. The activated pathways were arranged in the left column and the suppressed ones were in the right column. c Functional enrichment analysis showing molecular function (green), cell component (red), and biological process (blue) in GO categories of COPD-associated DEGs in six bulk RNA-seq datasets

PPI network analysis of MitoDEGs identifies 32 hub MitoDEGs associated with COPD

To further explore the potential role of mitochondria in COPD pathogenesis, we focused on transcriptomic changes of mitochondria in the context of COPD. COPD-associated DEGs of each dataset were intersected with mitochondria-related gene lists retrieved from the MitoCarta3.0 database to pick out MitoDEGs from the dataset, which added up to 353 MitoDEGs among all the datasets (Fig. 3a). To better characterize interaction topologies, we performed PPI network analysis of all these MitoDEGs, among which genes composing mitoribosomes, peroxisomes, and electron transport chain (ETC), as well as those participating in mitochondrial fusion/fission automatically clustered (Fig. 3b). Subsequently, using MCODE following the filter criteria as degree cut-off = 2, node score cut-off = 0.2, k-core = 2, and max depth = 100, a total of 32 hub MitoDEGs associated with COPD were identified from 4 significant modules, namely MRPL2, MRPL13, MRPL27, MRPL17, MRPL36, MRPL32, MRPL33, MRPL38, MRPL40, MRPL54, NDUFS2, NDUFS3, NDUFB7, NDUFAB1, NDUFA9, NDUFA7, NDUFA12, NDUFSA6, COX5A, UQCR10, EPHX2, ECH1, CAT, SCP2, MLYCD, HSD17B4, HMGCL, ABCD3, ABCD2, ATAD1, GDAP1, and FIS1 (Fig. 3c). Intriguingly, using ROC analysis of these hub MitoDEGs and other non-hub MitoDEGs from the 4 significant modules identified by MCODE, we noted that the 32 hub MitoDEGs were significantly better at distinguishing between COPD patients and control individuals compared to non-hub MitoDEGs in the cohort incorporating the 6 RNA-seq datasets of COPD lung tissues (GSE103174, GSE106986, GSE151052, GSE38974, GSE47460, GSE57148) (Fig. 3d), which strongly implied their potential in COPD diagnosis.

Fig. 3

PPI network analysis of MitoDEGs identifies 32 hub MitoDEGs associated with COPD. a Venn plot showing the intersection of mitochondria-related genes retrieved from the MitoCarta3.0 database (the area defined by the red line) with COPD-associated DEGs from six bulk RNA-seq datasets (6 areas defined by different background colors). The numbers of COPD-associated DEGs of each dataset intersected with mitochondria-related genes were marked in red. All the genes marked in red were defined as MitoDEGs. b The PPI network of all MitoDEGs. Genes (rectangular boxes) enriched in the same biological processes were arranged together in the same color. c The key modules identified by MCODE with 32 hub MitoDEGs (red). d ROC analysis of hub MitoDEGs (red) and other non-hub MitoDEGs (green) from the key modules identified by MCODE in (c) in the cohort incorporating the 6 RNA-seq datasets of COPD lung tissues (GSE103174, GSE106986, GSE151052, GSE38974, GSE47460, GSE57148)

Machine learning and MR analysis of hub MitoDEGs nominate three causal MitoDEGs for COPD

Next, to better clarify the potential association between mitochondria and COPD with higher credibility, we conducted an integrated analysis of the hub MitoDEGs combined with machine learning and MR. Specifically, the 18 genes first picked out as feature genes from the 32 hub MitoDEGs by RF analysis (Fig. 4a) displayed a strong discriminative power as reflected by the AUCs in ROC curve analyses based upon training set (AUC = 0.916), test set (AUG = 0.902), internal dataset GSE57148 (AUG = 0.881), and external dataset GSE76925 (AUG = 0.828) (Figure 4b). Subsequently, the 18 feature genes were further reduced to 13 by LASSO regression following the filter criteria as Log(λ) = 5.726 (Fig. 4c). Meanwhile, the eQTL data of the 32 hub MitoDEGs representing the genetic variants associated with the expression of these genes were downloaded from eQTLGen [22] and used as IVs to infer causal relationships between hub MitoDEGs (exposure) and COPD (outcome) by a two-sample MR analysis (Supplementary Tables S4 and S5). Among the 6 significant MR genes, including MRPL2 (95%CI: 1.000-1.002), MRPL27 (95%CI: 1.001-1.004), CAT (95%CI:0.999-0.999), NDUFS2 (95%CI:0.998-0.999), HSD17B4 (95%CI:1.000-1.001), and SCP2 (95%CI:1.000-1.00), only SCP2 exhibited significant horizontal pleiotropy (p = 0.029), and thus was excluded (Fig. 4d). The causal effects of each genetic variation for MRPL2, MRPL27, CAT, NDUFS2, and HSD17B4 on COPD were shown in the scatter plots (Fig. 4e). Notably, using the IVW method, only NDUFS2 and CAT were associated with COPD with OR < 1 as reflected by the slopes of the lines fitted, indicating the protective role of NDUFS2 and CAT against COPD. Whereas the other 3 causal genes (MRPL2, MRPL27, and HSD17B4) with OR > 1 may be detrimental in COPD pathogenesis. After removing each genetic variation for every causal MitoDEG, we systematically performed the MR analysis on the remaining genetic variations again and the results remained consistent with the overall association for each causal MitoDEG, further supporting the causality significant of every causal MitoDEG (Fig. 4f). However, we also noticed that a few error bars in the leave-one-out plots of CAT, NDUFS2, and HSD17B4 extended beyond 0, suggesting that there may be other potential influential genetic variations driving the causal links. Therefore, to improve the reliability of our analysis, we intersected the causal MitoDEGs for COPD nominated by MR analysis with the aforementioned 13 feature MitoDEGs successively identified by RF analysis and LASSO regression. Finally, only NDUFS2, MRPL2, and CAT appeared simultaneously in the results of machine learnings and MR analysis (Fig. 4g). Their contributions to COPD were further confirmed by another MR analysis using lung function (FEV1/FVC) as outcome measure (Supplementary Figure S2), aligning with standard COPD diagnostic criteria [35].

Fig. 4

Machine learning and MR analysis of hub MitoDEGs nominate three causal MitoDEGs for COPD. a The mean decrease in accuracy (MeanDecreaseAccuracy) (left) and mean decrease in Gini (MeanDecreaseGini) (right) were used to rank the relative importance of the 18 selected feature MitoDEGs for COPD which were selected by RF analysis following the filter criteria as both MeanDecreaseAccuracy and MeanDecreaseGini > 5. b Assessing the performance of the RF model using the receiver operating characteristic curve (ROC), the area under the curve (AUC) values were 0.916 in the training set, 0.902 in the test set, 0.881 in the internal validation dataset (GSE57148), and 0.828 in the external validation dataset (GSE76925), reflecting the overall stability and reliability of the RF model in (a). c LASSO regression showed log(λ) = 5.726 when the error of the model is minimized (top), and 13 MitoDEGs from 18 feature MitoDEGs identified by RF analysis in (a) were selected (bottom). d Forest plot showing the MR estimates of 6 hub MitoDEGs (MRPL27, MRPL2, NDUFS2, CAT, SCP2, HSD17B4) and COPD risk using the inverse variance weighting (IVW) method together with the results of sensitivity analyses. Pval represents p-value for MR analysis, P.pleio represents p-value for horizontal pleiotropy analysis, P.heter.ivw represents p-value for heterozygosity with IVW method, P.heter.mr represents p-value for heterozygosity with MR-Egger method. e Scatter plots showing the causal effect of MRPL27, MRPL2, NDUFS2, CAT, and HSD17B4 on COPD. f Leave-one-out plot demonstrating the causal effect of MRPL27, MRPL2, NDUFS2, CAT, and HSD17B4 on COPD when leaving one SNP out. g Venn diagram showing the intersection of 18 feature genes identified by RF analysis (green), 13 signature genes selected by LASSO regression analysis (blue), and 5 causal genes nominated by MR analysis (red). Only NDUFS2, MRPL2, and CAT appeared simultaneously in the results of three analyses

Reanalysis of causal MitoDEGs reveals their expression distribution and risk prediction value

To figure out the expression distribution of NDUFS2, MRPL2, and CAT, we enquired their expression levels in our previously published scRNA-seq dataset of human lung tissues from COPD patients and control individuals with different smoking histories (GSE171541) (Fig. 5a) [29]. Analysis of our scRNA-seq data revealed enrichment of NDUFS2, MRPL2, and CAT in mononuclear phygocyte system, especially macrophages (Fig. 5b). Nevertheless, the overall expression level of MRPL2 was not as high as those of NDUFS2 and CAT and a remarkable abundance of CAT was also observed in alveolar type 2 cells. Given the pivotal role of macrophages in orchestrating pulmonary inflammation during COPD development [36] and the proposed causal role in COPD pathogenesis of NDUFS2, MRPL2, and CAT identified by MR analysis, we sought to determine whether the alterations in the expressions of NDUFS2, MRPL2, and CAT in pulmonary macrophages were consistent with their predicted causal effect on COPD. According to the results of gene expression correlation analysis of NDUFS2, CAT, and MRPL2, none of them was correlated with the other two in macrophages (Supplementary Figure S3), therefore preliminarily ruling out the potential epistatic effects among MRPL2, CAT, and NDUFS2. The downregulation of NDUFS2 and CAT in AMs from COPD patients as compared with those from control individuals consisted with their protective roles predicted by MR analysis, whereas in an apparent paradox, the potential detrimental causal gene, MRPL2, was also decreased in AMs from COPD patients (Fig. 5c). Therefore, we constructed a nomogram model by logistic regression with the expression levels of NDUFS2, MRPL2, and CAT in pulmonary macrophages for the risk of COPD (Fig. 5d), the accuracy of which evaluated by the calibration curve indicated a good predictive capacity of this nomogram model (Fig. 5e). According to our nomogram model, the risk of COPD increased with the downregulation of NDUFS2, MRPL2, and CAT in pulmonary macrophages, again revealing the inconsistency between results of MR analysis and MRPL2 expression as regards risk for COPD. Note that NDUFS2 was characterized as a protective factor for COPD, whose downregulation not only was predicted to have a causal link with COPD by MR analysis but also contributed the most to risk for COPD among others in nomogram model, it was the focus of our following analyses.

Fig. 5

Reanalysis of causal MitoDEGs reveals their expression distribution and risk prediction value. a UMAP plot of the 14 identified cell types is used to visualize the result of clustering of cells from the scRNA-seq dataset GSE171541. One dot represents a cell. One color represents a cell type. b Feature expression of NDUFS2, MRPL2, and CAT on UMAP plots of all the cells. The color ranging from light red to dark red represents the expression level of the gene from low to high. c Violin plots showing the expression of NDUFS2, MRPL2, and CAT in human AMs isolated from COPD patients and control individuals of bulk RNA-seq dataset GSE130928. d The nomogram model for predicting risk for COPD with expression levels of NDUFS2, MRPL2, and CAT in pulmonary macrophages. The expression level of each gene is located on each variable axis, and a line is drawn on top to determine the number of points received for each variable value. The total points projected at the bottom scale indicate the risk for COPD. e The calibration curve of nomogram in (d). B = 1000 repetitions boot. Mean absolute error = 0.04

Functional analysis suggests the modulatory role of NDUFS2 in macrophage function and communication

To gain insight into the expression properties of NDUFS2 in pulmonary macrophages, a deeper examination of NDUFS2 expression in pulmonary macrophages taking smoking history into account in the context of COPD was conducted based upon our previous scRNA-seq dataset of COPD human lung tissues. Surprisingly, NDUFS2 increased significantly in pulmonary macrophages from non-COPD smokers compared to non-COPD non-smokers. However, under COPD condition, this upregulation of NDUFS2 in pulmonary macrophages in response to smoking was substituted by remarkable downregulation (Fig. 6a). GSVA highlighted the activation of pathways related to inflammation (such as TLR signaling pathway, TGF-β signaling pathway, NLR signaling pathway, chemokine signaling pathway, and JAK-STAT signaling pathway) as well as suppressed physiological immune function (such as antigen processing and presentation) in pulmonary macrophages with low level of NDUFS2 (macrophages ^NDUFS2low) compared to those with high NDUFS2 expression (macrophages ^NDUFS2high) (Fig. 6b). These differences were supported by GSEA, as inflammation response was enriched in macrophages ^NDUFS2low, whereas biological processes essential for physiological immune function of macrophages, including phagocytosis and endocytosis, were enriched in macrophages ^NDUFS2high (Fig. 6c). Considering the central role of macrophages in COPD pathogenesis as well as the multicellular microenvironment within lung tissues, we next investigated the intercellular communications of macrophages that were associated with NDUFS2 expression level by CellChat [37, 38]. Taking macrophages as sender cells (signal source) that drive outgoing intercellular communications, we first explored whether and how other cell types were influenced by macrophages ^NDUFS2low. Notably, outgoing cell-cell interactions (CCIs) of macrophages ^NDUFS2low with all the cell types including macrophages were increased compared with those of macrophages ^NDUFS2high (Fig. 6d), especially CCIs of macrophages ^NDUFS2low with innate immune cells, indicating the recruiting effect of macrophages ^NDUFS2low on other immune cells when orchestrating pulmonary inflammation. On the contrary, we also assessed the incoming CCIs of macrophages ^NDUFS2low under the condition that macrophages ^NDUFS2low were target cells (signal receivers) of intercellular signaling sent by other cells. Overall, incoming CCIs of macrophages ^NDUFS2low with all the cell types except for mast cells were increased compared with macrophages ^NDUFS2high, which was particularly obvious for those with structural cells and those with macrophages (Fig. 6d), suggesting the enhanced chemotactic activity of macrophages ^NDUFS2low in response to structural cells. Alterations of CCIs of macrophages ^NDUFS2low as signal source and receivers with other cell types respectively were more visualized after quantification (Fig. 6e). Meanwhile, the comparison of information flows in the CCIs networks for macrophages ^NDUFS2low and macrophages ^NDUFS2high showed considerable differences (Fig. 6f). Typical pro-inflammatory signals such as IL-6, IL-1, and TGF-β as well as chemokines such as CCL and CXCL families were dominant in CCIs networks for macrophages ^NDUFS2low, probably underling the increased outgoing CCIs of macrophages ^NDUFS2low with immune cells as well as incoming ones with structural cells. Intriguingly, signaling pathways of MHC-I and MHC-II were depressed in macrophages ^NDUFS2low compared with macrophages ^NDUFS2high, consistent well with the suppressed antigen processing and presentation in macrophages ^NDUFS2low highlighted by GSVA (Fig. 6b). Furthermore, altered outgoing and incoming CCIs of macrophages ^NDUFS2low from COPD group as compared with macrophages ^NDUFS2high from control group were similar to those between overall macrophages ^NDUFS2low and macrophages ^NDUFS2high (Fig. 6g, h). Changes of information flows accounting for these changed CCIs of macrophages ^NDUFS2low from COPD group as compared with macrophages ^NDUFS2high from control group also exhibited well consistency with those between overall macrophages ^NDUFS2low and macrophages ^NDUFS2high (Figure 6i). These findings implied that instead of being the outcome resulting from COPD development, altered CCIs of macrophages according to NDUFS2 expression level were very likely to mediate the vital role of macrophages in orchestrating the inflammation core to COPD pathogenesis, which corresponded well to the contributions of NDUFS2 downregulation to COPD inferred by MR analysis and nomogram model.

Fig. 6

Functional analysis suggests the modulatory role of NDUFS2 in macrophage function and communication. a Split violin plot showing the expression distributions of NDUFS2 in pulmonary macrophages among control and COPD patients with or without smoking histories. b Bar plot showing GSVA scores of the top altered KEGG pathways in macrophages ^NDUFS2low (blue) and macrophages ^NDUFS2high (green). c GSEA analysis highlighting activated inflammatory response and suppressed phagocytosis as well as endocytosis in macrophages ^NDUFS2low as compared with macrophages ^NDUFS2high. Circle plots (d) showing the altered outgoing (left) and incoming (right) intercellular communications of macrophages ^NDUFS2low as compared with macrophages ^NDUFS2high (the thicker the line, the greater the change of interaction), which was overall quantified and visualized as a heatmap (e), coupled with the significant signaling pathways being ranked based on their differences of overall information flow within the outgoing (left) and incoming (right) communication networks between macrophages ^NDUFS2low and macrophages ^NDUFS2high (f). Circle plots (g) showing the altered outgoing (left) and incoming (right) intercellular communications of macrophages ^NDUFS2low from COPD patients as compared with macrophages ^NDUFS2high from control individuals, which was overall quantified and visualized as a heatmap (h), coupled with the significant signaling pathways being ranked based on their differences of overall information flow within the outgoing (left) and incoming (right) communication networks between macrophages ^NDUFS2low from COPD patients and macrophages ^NDUFS2high from control individuals (i)

Expression of NDUFS2 was decreased in AMs from COPD mice

As unveiled by our analysis based on scRNA-seq dataset of human lung tissues that the expression level of NDUFS2 only significantly decreased in pulmonary macrophages from COPD smokers (Fig. 6a), we therefore verified the downregulation of Ndufs2 in cCS-induced COPD mouse model. After histopathological changes of COPD in lung tissues from CS-exposed mice being confirmed (Supplementary Figure S4a-4d), we measured the protein levels of Ndufs2 in lungs from control and COPD mice in situ. Immunofluorescence staining showed that NDUFS2 decreased significantly in pulmonary macrophages marked by F4/80 in lung sections from COPD mice as compared with control ones (Fig. 7a). A similar decline in the mRNA level of Ndufs2 was confirmed in AMs isolated from COPD mice by qPCR (Fig. 7b). Considering that NDUFS2 encodes a core subunit of the mitochondrial membrane respiratory chain NADH dehydrogenase (complex I), HSP60, a typical chaperone protein functions in mitochondrion) [39, 40], was used as mitochondrial marker to demonstrate the subcellular localization of NDUFS2 in AMs by immunofluorescence staining, which clearly showed the decrease of NDUFS2 in the mitochondria of AMs from COPD mice (Fig. 7c).

Fig. 7

Expression of NDUFS2 was decreased in AMs from COPD mice. a Eight-week-old C57BL/6 mice were exposed to CS to construct COPD mouse model, while the littermates exposed to room air for the same period served as control (n = 6). The protein levels of NDUFS2 (red) and F4/80 (green) in these mice lungs were measured by immunofluorescence staining (Scale bars = 50 μm). The overall expression levels of NDUFS2 and F4/80 in mouse lung sections were quantified as the red and green fluorescence intensity, respectively. The expression level of NDUFS2 in pulmonary macrophages (marked by F4/80) was quantified as the proportion of NDUFS2⁺F4/80⁺ cells in total F4/80.⁺ cells. **P < 0.01 vs. Ctrl (Student’s t-test). Scale bars = 50 μm. b, c Eight-week-old C57BL/6 mice were exposed to CS to construct COPD mouse model, while the littermates exposed to room air for the same period served as control (n = 6). The mRNA (b) and protein (c) levels of Ndufs2 were measured by qPCR and immunofluorescence staining in the primary AMs isolated from each mouse. The overall expression levels of NDUFS2 and HSP60 in AMs were quantified as the red and green fluorescence intensity, respectively. The representative images of AMs at low magnification were displayed (Scale bars = 20 μm) in the white solid box at bottom left corner of every enlarged image at higher magnification (Scale bars = 10 μm). The protein levels of NDUFS2 (red) in AMs (marked by HSP60, green) was quantified as red fluorescence intensity. **P < 0.01 vs. Ctrl (Student’s t-test)

Silencing NDUFS2 leads to pro-inflammatory phenotype and disrupted phagocytosis in macrophages

In vitro, we also confirmed the expression and localization of NDUFS2 in human macrophage cell line differentiated from THP-1 cells (hereafter referred to as “TP cells”) as well as mouse macrophage cell line RAW264.7 cells by immunofluorescence staining (Fig. 8a, b). We next silenced NDUFS2 in TP cells and RAW264.7 cells by siRNAs (Fig. 8c, d) to validate the effect of NDUFS2 downregulation on macrophage pro-inflammatory phenotype and phagocytosis. As expected, the increased expression of pro-inflammatory cytokines (IL6, STAT3, IL1B, TGFB, NLRP3, and IL23) as well as chemokines (CCL2, CXCL1, CXCL10, CXCL11, CXCL12, and CD80) were detected in both NDUFS2-silenced TP cells (Figure 8e) and Ndufs2-silenced RAW264.7 cells (Figure 8f). In addition, along with the dysregulation of genes involved in phagocytosis (MACRO, CD209, and SIRPA) (Figure 8e and 8f), NDUFS2 downregulation also impaired phagocytosis of TP cells and RAW264.7 cells which was assessed by Protonex^TM 600 Red Latex Beads (Fig. 8g, h). Together, these results supported the role of NDUFS2 in macrophages predicted by our bioinformatic analyses.

Fig. 8

Downregulation of NDUFS2 leads to pro-inflammatory phenotype and disrupted phagocytosis in macrophages. The protein levels of NDUFS2 (red) and HSP60 (green) were measured by immunofluorescence staining in TP cells (a) and RAW264.7 cells (b). Scale bars = 20 μm for representative images at low magnification in the white solid box. Scale bars = 10 μm for enlarged image at higher magnification. The protein levels of NDUFS2 (red) in AMs (marked by HSP60, green) was quantified as red fluorescence intensity. TP cells (c) and RAW264.7 cells (d) were transfected with or without 3 different pieces of siRNAs against NDUFS2/Ndufs2 for 48 h, when the mRNA and protein levels of NDUFS2/Ndufs2 were measured by qPCR (left) and western blotting (right). siNDUFS2-2 and siNdufs2-3 were selected for further experiments in TP cells and RAW264.7 cells, respectively. *P < 0.01, **P < 0.01 vs. siNC (Student’s t-test). e, f TP cells (e) and RAW264.7 cells (f) were transfected with siNDUFS2-2 and siNdufs2-3 for 48 h respectively y, when the mRNA levels of IL1B, IL6, TLR4, TGFB1, NLRP3, IFNG, IL23, STAT3, TNF, CD80, CXCL8, CXCL10, CXCL11, CXCL12, CCL2, CD209, MACRO, SIRPA were measured by qPCR. *P < 0.01, **P < 0.01 vs. siNC (Student’s t-test). g, h TP cells (g) and RAW264.7 cells (h) were transfected with siNDUFS2-2 and siNdufs2-3 for 48 h respectively, when the phagocytosis was assessed by Cell Meter™ Fluorimetric Phagocytosis Assay Kit. The images were taken under fluorescence microscopy (Scale bars = 20 μm). The red fluorescence intensity of phagocytosed beads reflects the digestion extent of the ingested particles in macrophages, while the average number of engulfed beads within every cell was calculated to indicate the engulfing ability of macrophages. **P < 0.01, ***P < 0.001 or vs. siNC. (Student’s t-test)

Discussion

In the present study, by adopting a multi-dimensional and layered integrated analytic strategy, together with experimental verification, we, to the best of our knowledge, for the first time nominated risk genes whose genetically regulated expression associated with mitochondria may have a causal role in COPD, providing compelling evidence that dysregulation of mitochondrial network in pulmonary macrophages underlies pulmonary inflammation in COPD pathogenesis.

We observed a ubiquity of dysregulated mitochondrial biology in the heterogeneous context of COPD by means of functional enrichment analysis of COPD-associated DEGs obtained from six COPD human lung tissues transcriptomic datasets respectively. Further intersecting these COPD-associated DEGs with the inventory of known human mitochondria-related genes from the MitoCarta3.0 database allowed us to reconstruct and unravel the disrupted mitochondrial network via PPI analysis of MitoDEGs. Notably, MitoDEGs automatically clustered into 4 subclusters related to mitoribosomes, mitochondrial fusion/fission, peroxisomes, and ETC, indicating the dysregulation of mitochondrial biogenesis, mitochondrial dynamics, mitophagy, and oxidative phosphorylation (OXPHOS) in COPD lung tissues, which was further emphasized by MCODE as key modules among PPI network of MitoDEGs. While in good agreement with previous studies of the mitochondrial malfunctions in COPD lung tissues [41], our findings were also conducive to pick out hub genes fundamental to these mitochondrial malfunctions responsible for COPD.

Deficiency in mitophagy, accounting for intracellular accumulation of aberrant mitochondria, was reported to increase reactive oxygen species (ROS) level and trigger inflammatory responses underling COPD [42]. Among the significant feature genes for COPD identified by RF and LASSO simultaneously, we noticed EPHX2 from the key module of peroxisomes, which encodes soluble epoxide hydrolase (sEH) in both cytosol and peroxisomes [43]. Capable of converting beneficial epoxylipids (such as anti-inflammatory epoxyeicosatrienoic acids) into corresponding diol metabolites with reduced bioactivity [43], sEH was found to adversely associated with multiple age-related cardiovascular diseases [44], neurodegenerative diseases [45], and inflammatory diseases [46], including COPD [47]. The protective role of sEH deletion or inhibition in these diseases were ascribed to attenuated inflammation [45, 47], less autophagy [47], preserved mitochondrial function [48], and alleviated cellular senescence [49]. It has been reported that CS-exposed Ephx2^−/− mice exhibited less pronounced pulmonary inflammation and autophagy accompanied by mild distal airspace enlargement, restored lung function, and steady weight gain compared with CS-exposed wild type mice [47]. Nevertheless, our analyses prompt questions about the role of sEH-mediated epoxylipids hydrolysis in mitophagy and COPD-associated mitophagy deficiency, which warrants further studies.

Altered mitochondrial morphology and dynamics, which was heavily intertwined with mitophagy, has been documented in the alveolar epithelial cells of COPD patients [50, 51]. Imbalanced mitochondrial fusion/fission leads to hyper-fused mitochondria, further interrupting mitophagy [52]. Unexpectedly, MitoDEGs such as MFN1 [53], MFF [54], and OPA1 [52, 55], which are well-known for their vital roles in mitochondrial quality control and widely accepted as dysregulated mitochondrial dynamic regulators in COPD lung tissues, were not identified as feature hub MitoDEGs for COPD in our analyses. Conversely, an unappreciated gene, ATAD1, was brought into our sight. Localized on the mitochondrial outer membranes (MOMs) with conserved specialized structural elements [56], ATAD1 helps to maintain mitochondrial proteostasis in various biological contexts by clearing out mistargeted proteins from mitochondria [57, 58] and extracting overload mitochondrial precursor proteins stuck in the protein translocase channel [59]. While the role of ATAD1 in COPD pathogenesis is still missing, we notice that the mitochondrial quality control exerted by ATAD1 can be extended to antiapoptotic and antiviral functions through the removal of pro-apoptotic protein BIM [60] and viral protein NS5B [61] from the MOMs, respectively. Thus, an interesting question to be investigated by future studies is whether and how dysfunction of ATAD1-mediated mitochondrial quality control contributes to pulmonary inflammation and alveolar destruction in COPD pathogenesis.

Indispensable to mitochondrial translation machinery, mitoribosome is essential for synthesizing ETC of the OXPHOS machinery [62]. Here, we identified many mitoribosomal proteins in MitoDEGs associated with COPD, reminiscent of the changed expression levels (biogenesis) as well as subcellular localization (assembly) of mitoribosomal proteins in COPD lung tissues revealed by previous studies [63, 64]. Despite the lack of direct link of mitoribosome proteins with COPD pathogenesis, we provided strong evidence that dysregulation of mitoribosome proteins may play an critical role in COPD development, not only because they served as hub genes (including MRPL2, MRPL13, MRPL27, MRPL17, MRPL36, MRPL32, MRPL33, MRPL38, MRPL40, and MRPL54) in the PPI network of MitoDEGs, but also because their genetic variations were found to have causal effect on COPD by our MR analysis. Thus, more investigations are required to delineate the crosstalk between mitochondrial translation relying on mitoribosome components biogenesis and assembly, OXPHOS machinery, and COPD pathogenesis.

As a causal gene with OR less than one by using the IVW method in MR analysis, CAT was proposed to be a protective factor for COPD, whose downregulation in macrophages was regarded as a risk for COPD in our nomogram model. Encoded by CAT, catalase is a key antioxidant enzyme defense against oxidative stress, which is predominantly located in peroxisomes [65]. Analogous to bronchiolar epithelium of human lung tissue from smokers and COPD patients, bronchiolar epithelium of CS-exposed mice exhibited a similar diminished expression of catalase [66]. However, up till now, no strong evidence for an association of polymorphisms in CAT accounting for decreased catalase expression and/or activity with COPD risk has been provided [67, 68]. In addition to peroxisomes, catalase is also localized in mitochondria [65] and mitochondrial catalase (mCAT) overexpression enables a 20% extension of lifespan in mice together with less mitochondrial deletion and delayed onset of heart diseases as well as cataracts [69]. Similar protective role of mCAT has been reported in alveolar epithelial cells against oxidative stress-induced mitochondrial DNA (mtDNA) damage and apoptosis [70, 71]. To note that mCAT is not limiting with regard to alveolar epithelial cells and may have a broader role in modulating pulmonary inflammation outside of antioxidation. For instance, taking advantage of reducing equivalents in the form of NADH, macrophages from transgenic mice overexpressed mCAT exhibited a sustained increase in NF-κB activation and expression of NF-κB dependent inflammatory mediators in response to LPS [72]. Indeed, we are still far from an entire understanding of mCAT in the lung. Nor is mechanistically clear how, precisely, mCAT regulate COPD pathogenesis.

NDUFS2, a subunit of Complex I in ETC, forms ubiquinone binding sites rather than FeS clusters [73]. Several clinical mutations in NDUFS2 [74,75,76], such as R138Q, Y141C, and R228Q, have been reported in patients and modeled in the homologous E. coli enzyme, which demonstrated reduced oxidase activity and aberrant assembly [73]. Although loss of one allele of Ndufs2 in C57BL6/J mice had no beneficial nor detrimental influence to healthspan and lifespan [77], another study showed, using a HEK293 cell line clone missing NDUFS2, that disruption of NDUFS2 resulted in decreased cell growth, Complex I specific respiration, glycolytic capacity, ATP pool, and cell-membrane integrity, together with significantly increased Complex II respiration, ROS generation, apoptosis, and necrosis [78]. Moreover, in pulmonary artery smooth muscle cells redox-sensitive NDUFS2 becomes chemically reduced during acute hypoxia, which consequently leads to hypoxic pulmonary vasoconstriction [79], while NDUFS2-dependent NAD⁺ regeneration is essential for directing cell fate during postnatal alveolar development [80]. Recent studies have suggested the link of NDUFS2 and inflammation, with mechanisms less clear though. Dexamethasone elevates the levels of NDUFS2 in the hippocampus following traumatic brain injury while reducing the expression of IL1 and inflammatory infiltration [81]. Nanoparticles exposure downregulates expression of subunits of ETC including NDUFS2, disturbing cellular respiration and redox balance, which may contribute to the increased release of inflammatory molecules from macrophages due to oxidative stress, DNA damage, inflammation responses, and cytotoxicity [82]. Herein, for the first time, by machine learning and MR analysis, NDUFS2 was highlighted as a causal gene for COPD, whose reduction in pulmonary macrophages was predicted by nomogram model to have the greatest contribution to COPD. It is intriguing that the expression of NDUFS2 in response to CS exposure exhibited remarkable distinction between pulmonary macrophages from COPD patients and those that from non-COPD individuals. Given the pro-inflammatory phenotype and depressed phagocytosis observed in NDUFS2-silenced macrophages, we took the upregulation of NDUFS2 in response to deleterious stimuli such as CS in pulmonary macrophages as an essential mechanism by which macrophages subdues inflammation and oxidative stress while sustaining physiological functions (including phagocytosis, efferocytosis, and antigen presentation) via enhancing energy production and maintaining redox balance. Thus, we reasoned that the intrinsic unresponsiveness of NDUFS2 expression against CS in pulmonary macrophages from smokers might be an unanticipated culprit for macrophage-orchestrated pulmonary inflammation that contributes to smoking-triggered COPD.

There are certainly limitations of our current work. First, despite the use of six public datasets of COPD human lung tissues, it is far from enough for stratified analysis based on factors such as age, smoking history, gender, and disease severity. Second, reanalysis of hub MitoDEGs should be conducted in more scRNA-seq datasets to fully investigate their roles in different cell types during COPD development. Third, in vivo experiments are required to validate the effect of NDUFS2 in pulmonary macrophages on intercellular communications and pulmonary inflammation during COPD development. Last but not least, different from machine learning alone, our bioinformatic analysis strategy that integrates machine learnings, MR analysis, and nomogram construction has a prominent advantage in screening out causal genes for COPD with genetic variations, which however also means that it trends to filter out potential causal genes without genetic variations. Therefore, other COPD feature hub MitoDEGs excluded by MR analysis should also be considered in future studies.

Conclusions

In summary, the current body of work reveals the causal relationship of mitochondrial malfunctions with COPD and nominates genetically dysregulated NDUFS2, CAT, and MRPL2 from hub MitoDEGs as significant risk factors for COPD. Combining the single-cell transcriptomic analysis with experimental validation allows us to mechanistically dissect the contribution of NDUFS2 in pulmonary macrophages to smoking-triggered COPD.