A Macrophage-Related Gene Signature for Identifying COPD Based on Bioinformatics and ex vivo Experiments

Introduction

Chronic obstructive pulmonary disease (COPD) is a respiratory disease characterized by progressive airflow limitation and chronic inflammation of the airway and lung parenchyma. It is the result of interactions between innate genetic factors, such as α1-antitrypsin deficiency, and acquired environmental factors, such as smoking.^1–4 Whether in developed or developing countries, the high incidence, disability and mortality of COPD have brought heavy economic burden to patients and society.^5,6 Diagnosis and timely treatment have critical clinical significance for slowing down the decline of pulmonary function and improving the long-term prognosis of patients.⁷ Nevertheless, relying solely on pulmonary function monitoring and imaging examination cannot fully identify the specific initial stage of COPD at present.^8–10 Undoubtedly, the accurate clinical diagnosis of COPD has yet to be resolved which is the utmost important.

In recent years, increasing studies have proved the essential role of immune cell infiltration in the occurrence and development of various diseases.^11,12 In the scope of COPD, Yang et al found that monocytes and M0 macrophages were up-regulated and CD8 T cells, activated NK cells, M2 macrophages, resting dendritic cells and resting mast cells were down-regulated in COPD lung tissues, while these differentially infiltrated immune cells were significantly associated with ferroptosis-related genes.¹³ Remarkably, even immune cells, such as monocytes, that have never been directly exposed to external stimuli, such as cigarette smoke or biomass smoke, have been observed to manifest defects in their function and regulation.^14,15 The expressions of SLC27A3 and STAU1 which were identified as the diagnostic markers of COPD by Zhang et al were enhanced in COPD models, accompanied by the activation of immune infiltration.¹⁶ However, although many biomarkers have been developed to identify COPD, a novel diagnostic model that focused on the immune microenvironment associated with COPD has not been constructed. With the development of technology, the fusion of bioinformatics analysis and machine learning (ML) has been applied to identify the diagnostic or prognostic genes related to COPD.^16–18 For example, Dhong et al identified TRPC6 as a diagnostic biomarker of Particulate matter- (PM-) induced COPD using machine learning models,¹⁹ but such studies still lacked validation of sufficient external datasets or relevant experiments to ensure stability.

In our investigation, two COPD microarray datasets were first extracted from NCBI Gene Expression Omnibus (GEO) database. Secondly, vital immune cells that were differentially infiltrated in COPD and healthy samples were identified and quantified by the cell-type identification by estimating relative subsets of RNA transcripts (“CIBERSORT”) algorithm, and a weighted gene co-expression network analysis (WGCNA) was used to screen out the key modules associated with these cells. Based on the candidate genes acquired from the intersection of key modules, we constructed a novel diagnostic prediction model using ML approaches and multivariate logistic regression analysis. Finally, we validated the diagnostic performance of our model on one internal test set, three external test sets, one single-cell transcriptome dataset, cigarette smoke extract- (CSE-) induced COPD models in vitro and lung tissue from COPD patients. This study aims to provide a new method for the diagnosis of COPD according to the changes in immune microenvironment, which is different from lung function and imaging examination.

Materials and Methods

Raw Data Acquisition

Datasets for three COPD airway tissues (GSE20257,²⁰ GSE10006²¹ and GSE8545)²², two COPD lung tissues (GSE57148^23–25 and GSE76925)^26,27 were downloaded from GEO (https://www.ncbi.nlm.nih.gov/geo/). All datasets above are gene expression arrays and were generated using the GPL570 (HG-U133_Plus_2) Affymetrix Human Genome U133 Plus 2.0 Array. GSE20257, GSE10006, GSE8545 contains 53 healthy and 23 COPD samples, 22 healthy and 27 COPD samples, 18 healthy and 18 COPD samples, respectively. GSE20257 was used as the training set for airway tissue, while GSE10006 and GSE8545 were used as the validation set for airway tissue. GSE57148 and GSE76925 were used to model the effect of external validation of lung tissue. All GEO raw datasets were pre-processed by “affy” in R, including background calibration, normalization and log2 transformation.²⁸ When multiple probes corresponded to a common gene, the average value was taken as its expression value.

Quantification of Immune Cells

The “CIBERSORT” inverse convolution algorithm (http://cibersort.stanford.edu/) was used to quantify 22 immune cells.^29,30 The “CIBERSORT” algorithm combined vector regression and LM22 gene set to estimate cell type proportions, incorporating samples with p < 0.05 into this work. We compared the proportion of immune cells between COPD and healthy samples. Immune cells with differential distribution and high predictive abilities were identified by the area under the curve (AUC) calculated via the “pROC” software package.³¹

Differentially Expressed Analysis

The “limma” package was applied to the GEO cohorts to screen for differentially expressed genes (DEGs). P-value < 0.05 and |log2FC| > 0.5 was set as the cut-off point.^32,33

Weighted Correlation Network Analysis

Weighted correlation network analysis (WGCNA) is a systems biology method for identifying gene relationship patterns across samples.³⁴ We chose the optimal soft threshold power β through the function “pickSoftThreshold”, where scale-free topological fitting index R2 = 0.9. aij = |Sij|β (aij: adjacency matrix between gene i and gene j; Sij: similarity matrix, which is obtained by Pearson correlation of all gene pairs) was used to calculate the matrix composed of weighted correlation values between genes and genes, ie, the adjacency matrix. Next up, the adjacency matrix was transformed into the topological overlap measure (TOM) matrix to reflect the similarity of co-expressed genes. We obtained the highly similar co-expressed genes and merged them into modules via the “dynamic tree cutting” algorithm and the “merged dynamic” algorithm. The fractions of vital immune infiltrating cells in COPD patients were used as sample traits. We identified the key modules according to a high correlation between these co-expression gene modules and sample traits. Additionally, the genes in the important modules were the key genes, and the intersected key genes between GSE20257 and GSE10006 were included as candidate genes in our work.

Functional Enrichment Analysis

Functional enrichment analysis, including Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) analysis, was performed by the R software packages “clusterProfiler”, “circlize”, etc.^35,36 The screening criteria for significance terms was P-value <0.05 after adjustment.

Construction of Protein–Protein Interaction (PPI) Network

We constructed the PPI network via STRING database (https://string-db.org/), and then the hub genes in the network were processed and visualized in “Cytoscape” software (version 3.9.1).³⁷

Construction and Validation of a Diagnostic Model

In GSE20257, the Boruta algorithm with 500 maxRuns was used to further screen for significant genes among the 41 candidate genes using the “Boruta” package. The genes were classified into three categories, including confirmed, tentative and rejected. Meanwhile, we put the 41 candidate genes into the least absolute shrinkage and selection operator (LASSO).³⁸ LASSO, a dimension reduction approach, shows superiority in evaluating high-dimensional data in comparison to regression analysis. LASSO analysis was implemented with a turning/penalty parameter utilizing a 10-fold cross-verification via glmnet package. The minimum lambda value was obtained by 10 times cross-validation. Nine genes identified as significant by Boruta as well as having non-zero regression coefficients were used for stepwise multivariate logistic regression analysis.³⁹ We constructed a score based on the 3 genes with P < 0.05 in the stepwise multivariate logistic regression. Score = Constant + Gene1 × Coef1 + Gene2 × Coef2 + Gene3 × Coef3 + - - - + Genen × Coefn (Constant, constant; Gene, gene expression; Coef, regression coefficient). The airway sequencing datasets GSE10006 and GSE8545 were used to verify the stability of this signature.

Processing of Single Cell Sequencing Data

We downloaded 10×scRNA-seq data containing three samples with COPD lung tissue from GSE167295.⁴⁰ Cells with fewer than 200 or more than 10,000 genes, fewer than 1000 molecules, and more than 20% mitochondrial and ribosomal genes were removed to retain high-quality cellular data. The remaining genes were expressed in at least three single cells. Above all, we used the “harmony” R package to remove batch effect between different samples⁴¹ and the “Seurat” R package to normalize the RNA-seq data via the function “NormalizeData”. After that, the standardized scRNA-seq data were transformed into Seurat objects and the top 2000 highly variable genes were identified using the function “FindVariableFeatures”. We applied the function “RunPCA” of the “Seurat” R package to perform principal component analysis (PCA) of the top 2000 genes based on dimensionality reduction of the RNA-Seq data. In the end, we employed “JackStraw” analysis to identify significant PCs and selected the top 15 PCs with the smallest p-value based on the proportion of variance for cell clustering analysis. The consolidated data were clustered using the functions “FindNeighbors” and “FindClusters”, and cell subsets were displayed via UAMP plot. The “SingleR” R package was used to annotate various cell subsets.⁴²

Cell Culture and Treatment

HBE cells, an SV40-transformed, normal HBE cell line, were obtained from Chi Scientific (Jiangsu, China). THP-1 cells were purchased from the National Collection of Authenticated Cell Cultures. THP-1 cells were treated with 10 µM PMA for 48 h to form THP-1 macrophages (THP-M cells), and THP-M cells were seeded in a six-well plate at a density of 1×106 cells per well. HBE cells were treated with 5% CSE for 24 h before they were cocultured with THP-M cells on 24 mm diameter inserts with 0.4 μm pores (#3412, Corning, USA).⁴³

Preparation of CSE

CSE was prepared as previously reported with some modifications.⁴³ Briefly, the smoke of a Da Qian Men (10 mg tar and 0.8 mg nicotine/cigarette, Shanghai, China) was bubbled through 10 mL of serum-free DMEM. The resulting suspension was adjusted to pH 7.4 and then filtered through a 0.22 μm pore filter (Merck Millipore, USA) to remove bacteria and large particles. The CSE was standardized by measuring the absorbance at 320 nm and defined as 100% CSE. The CSE was diluted to the desired concentration with medium and used in experiments within 30 min.

Determination of Cell Viability

Cell viability was assessed according to the instructions for Cell Counting Kit-8 kits (Beyotime, China). In short, 5000 cells per well were seeded in 96-well plates. 0, 2, 5, 10% CSE were added, and the preparations were incubated at 37°C and 5% CO2 for 0, 6, 12, 24, 48, 72 h. CCK-8 solution (10 μL) was added to each well, and, after incubation for 2 h, the absorbance was measured at 450 nm. For following experiments, we exposed HBE cells to 5% CSE for 24 h based on the outcome (Supplementary Figure 1).

RNA Extraction and Real-Time PCR

Total RNA was isolated by use of RNAiso Plus (9108/9109, Takara, Japan). Total RNA (1 μg) was transcribed into cDNA by use of PrimeScriptTM RT reagent kit with gDNA eraser (RR047A, Takara, Japan). The RT-PCR assay was performed with TB Green™ Premix Ex Taq™ II (RR820A, Takara, Japan) and ABI 9600 real-time PCR detection system (Applied Biosystems). GAPDH served as internal control for mRNA. Fold changes in expression of each gene were calculated by a comparative threshold cycle (Ct) method using the formula 2-ΔΔCT. The primer sequences are as follows (Table 1).

Table 1 Primer Sequences Used

Western Blot Assay

Total lysates were prepared according to the manufacturer’s recommendations (Beyotime Institute of Biotechnology, Shanghai, China). Protein concentrations were measured with the BCA protein assay (Beyotime Institute of Biotechnology, Shanghai, China). Equal amounts (20 μg) of protein were separated by 10% sodium dodecyl sulphate–polyacrylamide gel electrophoresis and were transferred to polyvinylidene fluoride (PVDF) membranes (Millipore, Billerica, MA). Membranes were incubated overnight at 4°C with mouse anti-GAPDH antibody (ab8245, Abcam), goat anti-CLEC5A antibody (AF2384, R&D systems), rabbit anti-FTL antibody (10,727-1-AP, Proteintech) and rabbit anti-SLC2A3 antibody (ab41525, Abcam). After several washing steps, the membrane was incubated with secondary horseradish peroxidase (HRP)-conjugated antibody at room temperature for one hour. Detection was performed with the Immobilon ECL system (Millipore, S.p.A., Italy). The densitometric analyses of the bands were performed with ImageJ software.

Immunohistochemical Staining

The lung tissues were fixed with paraformaldehyde, then were embedded in paraffin and sectioned into sections. The lung sections were deparaffinized and rinsed at room temperature. The primary antibody used was the same as in Western Blot. After incubated with the appropriate HRP-conjugated secondary antibodies, the samples were incubated with a 3.39-diaminobenzidine-4HCl (DAB; Sigma, Missouri, USA)-H2O2 solution to visualize immunolabeling. The lung sections were mounted with a coverslip and observed under a microscope.

Statistical Analysis

All statistical analyses were performed in R (version 4.2.1). The predictive power of immune cells and the efficacy of the diagnostic model were assessed using receiver operating characteristic (ROC) curves. The Wilcoxon analysis is used to analyze the differences between the two groups, including the comparison of 22 immune cells and candidate signatures between COPD patients and healthy controls. Statistical significance was identified based on P < 0.05.

Results

Identification of Important Immune Cells in COPD

In order to explore immune microenvironment landscapes in healthy controls and COPD patients, we used the “CIBERSORT” algorithm to quantify the relative proportion of infiltrating immune cells in GSE20257 and GSE10006 (Figure 1A and B). The correlation between immune cells may provide clues to better understand the composition of the immune microenvironment in COPD (Supplementary Figure 2A, B). In GSE20257, the fraction of resting dendritic cells and M0 macrophages were significantly higher in the COPD patients than those in healthy controls, but lower proportion of M2 macrophages was found in COPD patients (Figure 1C, P < 0.05). In GSE10006, we found the healthy controls had a higher immune infiltrating level of plasma cells, follicular helper T cells, and gamma delta T cells compared to those in COPD patients, while lower proportions of M0 macrophages and resting dendritic cells were found in healthy controls (Figure 1D, P < 0.05). Obviously, the expression of M0 macrophages and resting dendritic cells was consistent throughout the two cohorts. In addition, receiver operating characteristic (ROC) curve showed that M0 macrophages had high predictive accuracy in GSE20257 (Figure 2A, AUC > 0.7). In GSE10006, plasma cells, follicular helper T cells, M0 macrophages, and resting dendritic cells all had good predictive values (Figure 2B, AUC > 0.7). However, as shown in Figure 1, M0 macrophages accounted for a significantly larger proportion of infiltrating immune cells than dendritic cells, follicular helper T cells, etc. Therefore, M0 macrophages were considered to be the most important immune cells of COPD in this work.

Figure 1 Immune microenvironment of airway tissues in healthy controls and COPD patients. (A) Relative proportions of infiltrating immune cells in GSE20257. (B) Relative proportions of infiltrating immune cells in GSE10006. (C) Differences of immune cells composition between healthy controls and COPD patients in GSE20257. (D) Differences of immune cells composition between healthy controls and COPD patients in GSE10006. Blue represents healthy controls; red represents COPD patients. ns, no significance; *P < 0.05, **P < 0.01.

Figure 2 Identification of key immune cells and DEGs. (A) AUCs of three differentially distributed immune cells in GSE20257. (B) AUCs of five differentially distributed immune cells in GSE10006. (C) Volcano plot of DEGs in GSE20257. (D) Volcano plot of DEGs in GSE10006.

Determination of Key Modules Associated with Important Immune Cells

The R package “limma” was utilized to analyze differential expression genes in the metadata cohorts (GSE20257 and GSE10006). As shown in Figure 2C, 773 up-regulated and 726 down-regulated DEGs were identified in GSE20257. In GSE10006, 348 up-regulated genes and 449 down-regulated genes were visualized by volcano plot (Figure 2D). Next, we took advantage of the above DEGs to determine the key modules associated with important immune cells in two data sets, respectively. When the soft threshold power β was 8, scale-free topology fitting index R2 reached 0.9 in GSE20257, then we acquired 6 modules using the “merged dynamic” algorithm (Figure 3A and C). We found that the green module had a good correlation with M0 macrophages (Figure 3E, cor = 0.41, P < 0.001). Similarly, in GSE10006, we chose β = 7 as the optimal β and acquired 5 modules (Figure 3B and D). The brown module was highly correlated with M0 macrophages among the modules (Figure 3F, cor = 0.67, P < 0.001). Hence, the green module in GSE20257 and the brown module in GSE10006 were chosen as key macrophage M0-related modules.

Figure 3 Determination of key modules associated with important immune cells. (A) Scale Independence and average connectivity in GSE20257. (B) Scale Independence and average connectivity in GSE10006. (C) Cluster dendrogram in GSE20257. (D) Cluster dendrogram in GSE10006. (E) Heatmap of correlation between modules and important immune cells in GSE20257. (F) Heatmap of correlation between modules and important immune cells in GSE10006.

Enrichment Analysis of Key Genes and Construction of PPI Network

We took the intersection of the genes between the above two modules and ultimately obtained 41 candidate genes (Figure 4A). The results of GO enrichment analysis, consist of 240 biological processes (BPs), 9 cellular components (CCs) and 14 molecular functions (MFs), indicated that these candidate genes were mainly involved in immune-related activities, including regulation of cytokine production, regulation of immune effector process, regulation of innate immune response, specific granule, secretory granule membrane and immune receptor activity (Figure 4B, Supplementary Table 1). The KEGG enrichment analysis exhibited that these genes mostly participated in phagosome, lysosome, neutrophil extracellular trap formation, B cell receptor signaling pathway, necroptosis and even ferroptosis (Supplementary Figure 3A, Supplementary Table 2). Finally, we explored the interactive relationship among the proteins encoded by candidate genes in STRING database, and the hub genes in the network were screened and visualized using “Cytoscape” software. The results existed 29 nodes and 134 edges in total (Supplementary Figure 3B, Supplementary Table 3).

Figure 4 LASSO regression analysis for selection of trait genes. (A) 41 candidate genes were acquired via Venn diagram. (B) Donut plot of GO enrichment analyses. (C) Ten cross-validations of adjusted parameter selection in the LASSO model. Each curve corresponds to one gene. (D) LASSO coefficient analysis. Vertical dashed lines are plotted at the best penalty parameter.

Construction and Validation of a Macrophage M0-Related Diagnostic Model

Firstly, 41 candidate genes were downscaled by the “Boruta” algorithm to classify them into two categories, including 15 confirmed genes and 26 rejected genes (Table 2). Next, to obtain higher accuracy, the 41 candidate genes were selected after ten cross-validations by the “LASSO” algorithm with an optimal penalty parameter of 0.0038. We selected the minimum criteria for building the LASSO classifier and identified 14 feature genes (Figure 4C and D, Table 2). Two machine learning yielded 9 genes that were used in our next step of stepwise multivariate logistic regression analysis (Table 2). Finally, based on the results of the stepwise multivariate logistic regression analysis, the P-values of three genes were less than 0.05, and we constructed a new model based on the expression values of these three genes and the coefficients of the stepwise multivariate logistic regression analysis. Score = −43.5741 + CLEC5A * 1.5158 + FTL * 4.7296 - SLC2A3 * 1.4324.

Table 2 Detailed Results of the Boruta and LASSO Algorithms

Compared with healthy controls, COPD patients tended to have higher scores in training group (GSE20257) and two validation groups (GSE10006 and GSE8545). In other words, a higher score implied the higher expressions of CLEC5A, FTL and SLC2A3 and a greater likelihood of suffering from COPD (Figure 5A, D and G). Most of all, the AUC was 0.916, 0.873 and 0.932 in GSE20257 (Figure 5B), GSE10006 (Figure 5E) and GSE8545 (Figure 5H), respectively, indicating our model has an excellent value in diagnosing COPD. In these three GEO cohorts, the expressions of CLEC5A, FTL and SLC2A3 in COPD patients were all higher than the healthy controls (Figure 5C, F and I). To broaden the applicable value of this model, we used two lung tissue sequencing data (GSE57148 and GSE76925) for external verification. The AUCs were 0.675 in GSE57148 and 0.688 in GSE76925 (Supplementary Figure 4A–D). Obviously, our model still had certain predictive ability.

Figure 5 Identification of a 3-gene signature. Different results and expressions of CLEC5A, FTL and SLC2A3 between high- and low-score groups in GSE20257 (A), GSE10006 (D) and GSE8545 (G). AUC of predicted outcome in GSE20257 (B), GSE10006 (E) and GSE8545 (H). Expressions of CLEC5A, FTL and SLC2A3 in COPD patients and healthy controls in GSE20257 (C), GSE10006 (F) and GSE8545 (I). ns, no significance; **P < 0.01, ***P < 0.001.

Expression Levels of CLEC5A, FTL and SLC2A3 in Single Cell Sequencing Data

Based on the scRNA-seq data of GSE167295, we obtained gene expression profile data of peripheral lung parenchymal tissue from three COPD patient samples for further analysis. After stringent quality control filtering to remove low-quality cells, the data were removed from batch effects and normalized. We used the “harmony” R package for batch effects (Supplementary Figure 5A and B). The PCA method was used to downscale, and the top 15 PCs were selected for further analysis based on P-values (Figure 6A). After obtaining each cell cluster, we combined the “singleR” algorithm to annotate cell subpopulations, and UMAP was used to display specific cell subpopulations (Figure 6B). The heatmap shows the relative expression of marker genes in each cell subpopulation (Figure 6C). The expression of CLEC5A, FTL and SLC2A3 in each cell subpopulation was displayed using featureplot function. The results showed that these three genes were indeed expressed in the macrophages from lung tissue of COPD patients: the expression level of CLEC5A was very low and was almost only expressed in macrophages; although FTL was expressed in all cell subsets, the expression intensity was highest in macrophages; while the expression of SLC2A3 was the highest in monocytes (Figure 6D and E).

Figure 6 Analyses of single cell RNA sequencing data from lung tissue of COPD patients. (A) PCA was used for linear dimension reduction and the top 20 PCs were selected according to p-value. (B) UMAP plot of 9 cell subgroups shown by different colors. (C) Heat map showed the relative expression of marker genes for the 9 cell subgroups. Red represents high expression; blue represents low expression. (D and E) Expression of CLEC5A, FTL and SLC2A3 in the 9 cell subgroups.

Expression of CLEC5A, FTL and SLC2A3 in vivo and in vitro

To simulate CSE-induced airway immune microenvironment, THP-M cells were differentiated from THP-1 monocytes and then co-cultured with 5% CSE-treated HBE cells for 24 h. The consequences of qRT-PCR and Western blot showed that the mRNA and protein levels of CLEC5A and SLC2A3 were up-regulated in a time-dependent manner in THP-M cells co-cultured with CSE-treated HBE cells, while the expression of FTL in protein level was up-regulated, but there was no significant change in mRNA level (Figure 7A–E).

Figure 7 Expression of CLEC5A, FTL and SLC2A3 in THP-M cells co-cultured with CSE-treated HBE cells for 24 h. (A) Expression of CLEC5A mRNA. (B) Expression of FTL mRNA. (C) Expression of SLC2A3 mRNA. (D) Western blots. (E) Relative protein levels of CLEC5A, FTL and SLC2A3 were determined. The data are the mean ± SD (n=3). *P < 0.05.

Not only that, but we evaluated CLEC5A, FTL, and SLC2A3 expression levels in the lung tissue of non-smokers, smokers, and COPD patients. The expression levels of these three proteins in the lungs of COPD patients, especially CLEC5A, were prominently higher than those of non-smokers and smokers (Figure 8A and B). Immunohistochemical staining results showed that CLEC5A, FTL and SLC2A3 were expressed in pulmonary macrophages from COPD patients with elevated expression levels compared with non-smokers or smokers (Figure 8C).

Figure 8 Expression of CLEC5A, FTL and SLC2A3 in COPD patients. (A) Western blots. (B) Relative protein levels of CLEC5A, FTL and SLC2A3 in lung were determined. Con-NS, nonsmokers without COPD; Con-S, smokers without COPD; COPD, COPD patients. Data are mean ± SD (n = 3). *P < 0.05, ***P < 0.001. (C) Immunohistochemical staining for CLEC5A, FTL and SLC2A3 in non-smokers, smokers, and COPD patients. Positive protein staining appeared brown and nuclear staining appeared blue. Red arrows pointed macrophages. Scale bar = 50μm, original magnification: ×400.

Discussion

Global initiative for chronic obstructive lung disease (2023 REPORT) first wrote the concept of “preserved-ratio impaired spirometry” (PRISm), that is, the one-second rate is normal (FEV1/FVC ≥ 0.7 after inhalation of bronchodilators) but lung ventilation function is impaired (FEV1% and/or FVC%<80% after inhalation of bronchodilators).⁴⁴ Due to the mighty compensatory capacity of human lung and the absence of valuable biomarkers,85% of COPD patients in China have had moderate or above airflow limitation at the time of diagnosis. COPD will turn into the third leading cause of death in the world by 2030.⁴⁵ It follows that there is an urgent demand to develop early diagnosis models for COPD. Variations in the pulmonary immune microenvironment have been considered to play a crucial role in the immunomodulatory response and the occurrence and development of COPD at present.^12,46–48 Therefore, taking immune infiltration as the starting point, it has broad prospects to search for novel and effective signatures for COPD, which potentially consequence means of intervening early in diagnosis and treatment, thereby prolonging the survival of patients.

It is well known that pulmonary inflammation is the core pathogenesis of COPD, and smoking is considered to be the main incentive for COPD. Considering that the airway and even the alveolar are the sites directly exposed to cigarette smoke, we selected two airway tissue sequencing datasets (GSE20257 and GSE10006) for the study, where GSE20257 only contains sequencing data from small airways, while GSE10006 contains data from atmospheric ducts and small airways. We determined the infiltration of immune cells in samples from GSE20257 and GSE10006 via the “CIBERSORT” algorithm and the LM22 gene set. Through the differential analysis and evaluation of the ROC curve, we found that the fraction of M0 macrophages increased in COPD patients compared with healthy controls in both cohorts and M0 macrophages had satisfactory predictive value. Therefore, we considered M0 macrophages as the critical immune cells in COPD. Many studies have demonstrated that macrophages are the major inflammatory cell type in COPD, and the number of macrophages had a significant increase in sputum, bronchial biopsies and bronchoalveolar lavage fluid (BALF) of COPD patients.⁴⁹ Macrophages are the main effector cells for uptake and processing of respiratory pathogens and inhalable particles. While COPD macrophages had defects in activity, adhesion, and phagocytosis, which was likely to promote further pulmonary recruitment of macrophages.^50,51 Human rhinovirus (HRV) can impair the phagocytosis of COPD macrophages and enhance their responses to cytokines via the TLR3/IFN pathway, thereby preventing the resolution of inflammation and leading to prolonged exacerbations in COPD.⁵² Besides, there was also a great difference in the composition of macrophages between COPD patients, smokers, and non-smokers; the proportion of M1 macrophages in smokers and COPD patients was lower than that of M0 and M2 macrophages.⁵³ The findings above further illustrated the importance of M0 macrophages to COPD. It should be noted that the differences in polarized macrophages between COPD patients and healthy controls were not statistically significant in this work. We considered that this result may be related to the failure to analyze the COPD patients according to age, smoking history, GOLD grade, disease phenotype, presence or absence of acute exacerbation, on account of lack of basic information and clinical data.

Forty-one candidate genes were obtained from the intersection of two key modules that are highly correlated with M0 macrophages. Functional enrichment analysis of these candidate genes revealed that almost all the biological function items obtained, such as regulation of cytokine production, regulation of innate immune response, specific granule, phagosome, lysosome, ferroptosis, etc., were related to the occurrence and development of COPD. For instance, Kearley et al found that after cigarette smoke exposure or viral infection, IL-33 can amplify the function of macrophages, causing them to produce more cytokines to aggravate local inflammation.⁵⁴ CSE-treated neutrophils were effectively recognized and non-inflammatory phagocytosed by monocyte-derived macrophages, increased degranulation of which can affect the immune response of COPD patients.⁵⁵ An increased number of alveolar macrophages with abnormal phagosome and lysosome functions was defined as a characteristic change in COPD. HRV16 can perturb phagocytosis in macrophages by downregulating the expression of Arpin, a key factor in targeted phagocytosis of macrophages.⁵⁶ Furthermore, when macrophages ingested some toxic particles, such as cigarette smoke, it led lysosomes to release excessive proteases and other toxic enzymes into the cytoplasm, resulting in dying or dead macrophages.⁵⁷ Previous study by our research group also found that NCOA4-mediated ferroptosis was involved in the polarization of COPD macrophages.⁴³ Obviously, the functions of 41 candidate genes were associated with most macrophage M0-associated activities. To sum up, we speculated that these candidate genes may be involved in the occurrence and development of COPD through immune-related activities, and they may become the potential targets for diagnosing or treating COPD.

Finally, we constructed a predictive gene signature consisting of three genes (CLEC5A, FTL and SLC2A3) based on the gene expression profile in GSE20257, using ML approaches and multivariate logistic regression analysis. All three genes were highly expressed in COPD patients. C-Type lectin receptor 5A (CLEC5A), a spleen tyrosine kinase- (Syk-) coupled pattern recognition receptor to bacterial and viral infections, is involved in the innate immune response. Studies have revealed that extracellular vesicles released from activated platelets can amplify the inflammatory signal induced by dengue virus via CLEC5A/TLR2 on macrophages and promote macrophages to secrete proinflammatory mediators such as TNF-α and IL-6.^58,59 Importantly, even under sterile conditions, once CLEC5A on non-inflammatory macrophages was targeted for activation, it would elicit a response similar to that against dengue virus and lead to a shift of myeloid cells to a proinflammatory phenotype.⁶⁰ Besides, smoking can increase the expression of CLEC5A in immune cells.⁶¹ Wortham et al demonstrated that CLEC5A was highly expressed on alveolar macrophages in mouse models of COPD and human long-term smokers, while it was almost not expressed in macrophages of healthy controls,⁶² which was consistent with our findings. What’s interesting was that Sung et al found that thrombotic inflammation was attenuated dramatically in clec5a−/− mice,⁶³ based on which we dared to speculate that perhaps the high expression of CLEC5A is a major risk factor for embolism in COPD patients. Ferritin light chain (FTL) has a powerful immunomodulatory function, which can prevent organ injury caused by excessive inflammatory response by inhibiting NF-κB activation.⁶⁴ Meanwhile, FTL played a protective role in PM-induced lung injury by maintaining the stability of lysosomal function in macrophages.⁶⁵ Tumor cells would produce FTL-rich vesicles to affect the recruitment of macrophages in hypoxic environment, so as to better escape host immunity.⁶⁶ The encoding product of SLC2A3 gene (GLUT3) is a member of glucose transporter (GLUT) family, which is widely distributed in various human cells and closely related with glycometabolism. Numerous studies have confirmed that GLUT3 was highly expressed in pulmonary inflammatory diseases.^67–69 Hochrein et al found that GLUT3 controlled T helper 17 cell-mediated inflammatory responses via glycolytic-epigenetic reprogramming.⁷⁰ These results further emphasize the vital function of FTL and SLC2A3 in inflammation and immunity. Our macrophage M0-related diagnostic model achieved high AUCs in both the training set and the validation sets derived from airway tissue sequencing but was not well represented in the validation set from lung tissue. We considered that this may be related to the excessive complexity of cell groups in lung tissue and the small proportion of immune cells. To make the model more convincing, we further successfully demonstrated that CLEC5A, FTL and SLC2A3 are indeed expressed in the macrophages from COPD patients and at higher levels compared with smoking or non-smoking healthy controls by single-cell data analysis and experiments in vitro and in vivo. Overall, the three genes that make up our diagnostic model are able to serve as biomarkers to identify COPD patients.

Although the concept of M0/M1/M2 macrophages seems outdated in the omics era, the LM22 signature we used as a basis for the CIBERSORT algorithm to identify cell types does not include certain cell subsets that are known to play a role in COPD according to literature (like MDSCs) because it is derived from sorted human blood immune populations. We expect that more sophisticated concepts and analytical methods will emerge to construct relevant genetic models in the future. In addition, notwithstanding our model provides a method for the diagnosis of COPD that is different from lung function and imaging examination, the combination of genetic diagnosis model and analysis of important clinical data can effectively improve the sensitivity and specificity of disease identification.⁷¹ Limited by its lack of relevant information, our model should be further validated in large-scale clinical data to ensure the stability of the results and to determine whether it can be applied to the early diagnosis of COPD. As previous research has indicated that cigarette smoke can have a significant impact on the transcriptomics of tracheal epithelium. Notably, HBE cell lines are non-differentiated, and it would be valuable if we could continue to explore the impact of CSE on the expression of CLEC5A, FTL, and SLC2A3 in differentiated epithelium and compare it to undifferentiated epithelium (basal cells) in subsequent studies.⁷² Additionally, differences in gene expression between males and females in response to chronic injury in the airway epithelium could have implications on disease susceptibility, both in COPD and other respiratory diseases.⁷³ Considering that there are certain differences in the therapy of COPD patients with different phenotypes.^74,75 Our subsequent research will also focus on applying ML to the early differentiation of COPD phenotypes and assessing the sensitivity to glucocorticoid therapy in COPD patients with different phenotypes. At present, deep learning has been applied to drug development and vaccine design,⁷⁶ and we should also keep an eye on this in order to find potential drugs that can benefit COPD patients. Finally, although we have discovered that CLEC5A, FTL, and SLC2A3 help to identify COPD, the specific mechanisms by which these three genes impact COPD remain to be explored.

Conclusion

To sum up, our study clearly revealed the relationship between COPD and immune cells, in which M0 macrophages were essential. More than anything, we constructed and validated a novel macrophage M0-associated diagnostic model with satisfactory predictive value. Immunoinfiltration has always been a research hotspot of novel therapies for COPD. CLEC5A, FTL and SLC2A3 are expected to be promising targets of immunotherapy in COPD.