The scRNA dataset (GSE173896) and the RNA-sequencing (RNA-seq) dataset (GSE57148) were obtained from GEO (http://www.ncbi.nlm.nih.gov/geo). The samples in both datasets were lung tissues from COPD and control subjects.

Four COPD (GSM5282538, GSM5282539, GSM5282540, and GSM5282541) and three control nonsmokers (GSM5282546, GSM5282547, and GSM5282548) scRNA-seq lung tissue samples were obtained from GSE173896. The total number of cells in the COPD and control groups was approximately the same. Afterward, the Seurat package initially filtered scRNA-seq data to obtain higher-quality cells (11). For further analysis, cells and genes must fulfill the three quality control criteria: 1. The detected gene in at least three cells. 2. Detected at least 200 genes per cell. 3. No more than 15% of the genes are mitochondrial genes. With quality control, the overall count of cells in the COPD and control groups was 11,024 and 11,225, respectively. The canonical correlation analysis (CCA) algorithm integrated the seven samples and removed batch effects. Then, 3,000 highly variable genes were chosen for analysis, and the number of principal components (PCs) was programmed to 50 to get clusters of cells, which were visualized as a “UMAP” plot (12). All cells were subsequently labeled by the marker genes for each cell type. FindMarker was used to identify differential genes.

The Monocle2 was applied to perform single-cell trajectory analysis (13). All neutrophil objects were extracted, and the DDRTree method was used to reduce the dimension of cells. Then, the reduce Dimension function was implemented to determine the type of cell differentiation state. Finally, the “plot cell trajectory” function was employed to visualize the differentiation trajectory of cells.

Bioinformatics (https://www.bioinformatics.com.cn/) website was an open online platform for data processing and visualization. Gene Set Enrichment Analysis (GSEA) was used to identify enriched GO terms and KEGG associated with cell differentiation state. Sorting all of the differential genes according to logFC and GSEA was conducted on the website.

The level of gene set enrichment in each cell was assessed using the irGSEA package. This package was used to score individual cells and to generate numerous gene set enrichment score matrices. The half vlnplot and density scatterplot visualized some specific enriched pathways.

We calculated the extent of immune cell infiltration in COPD and control lung tissue samples in GSE57148 by using the CIBERSORT algorithm. This algorithm set the PERM to 1000 and set the cutoff to p < 0.05 (14). The Pheatmap package created a heat map of different immune cells, and the Vioplot package displayed the abundance. A correlation heat map was created to visualize the correlation among infiltrating immune cells. Furthermore, the correlation of critical genes with the differential immunity cells was analyzed.

Weighted gene co-expression network analysis (WGCNA) was performed to construct modules with similar expression pattern genes on GSE57148 (15). Firstly, outlier samples and low-expressed genes were eliminated. Secondly, the adjacency matrix was built based on an optimal soft threshold and converted into a topological overlap matrix. Then, a hierarchical clustering tree was constructed to classify the highly co-expressed genes into the same module. Finally, COPD and neutrophils were identified as characteristics to compute the correlation between gene modules and traits.

COPD samples (n = 96) and normal samples (n = 90) in the GSE57148 dataset were analyzed using the limma package to screen for differentially expressed genes (DEGs), with an adjusted p-value of < 0.05 and | log2(fold change, FC) | > 0.4 as thresholds.

GO analysis was a popular approach for massive functional enrichment, encompassing biological process (BP), cellular component (CC), and molecular function (MF). KEGG was a globally accessible database storing essential data on genomes, biological pathways, diseases, and drugs. GO terms and KEGG pathway analysis were annotated for gene function by the clusterProfiler package, and the top-ranked results were visualized (16).

The MR analysis was to investigate the causality between DEGs and the risk of COPD, where the instrumental variables (IVs) were single nucleotide polymorphisms (SNPs). The cis-eQTL data of DEGs were obtained from eQTLGen, and GWAS statistics for COPD outcome were obtained from integrated epidemiology unit (IEU) database (ukb-b-20464) (17, 18). The cis-eQTL data were used as exposure factors and SNPs with P-value < 5 × 10^–8 were used for the MR analysis. If no suitable SNPs were available at the stringent threshold, the threshold was loosened to 5 × 10^–6. There were still no relevant SNPs, so this cis-eQTL data was not used for MR analysis. All SNPs were identified 100 kb upstream and downstream of the corresponding gene location, with a threshold for Linkage disequilibrium (LD) of R2 < 0.2 and a physical distance threshold of 250 kb.

The TwoSampleMR package performed MR analysis. MR estimates were calculated for each SNP using the Wald ratio method. If more than one SNP was available, MR estimates were weighted by the inverse variance of the ratio estimates (inverse variance weighted, IVW). Cochran’s Q statistic was used to test for heterogeneity (19). When the number of SNPs was no less than three, MR-Egger regression evaluated potential horizontal pleiotropy (20).

The protein-protein interaction (PPI) network was built on the STRING online site (21). An interaction score greater than 0.4 (medium confidence score) was considered meaningful. Then, the PPI network was further analyzed in Cytoscape. CytoHubba, a Cytoscape plugin that discovered hub terms and sub-networks from interaction groups, was used to analyze PPI network hub genes. The top 10 genes were confirmed using the maximal clique centrality (MCC) computing method with the Cytohubba plugin.

Peripheral blood samples in this study were obtained from Tongji Hospital, Wuhan, China. All subjects signed an informed consent form. The Ethics Committee of Tongji Hospital, Huazhong University of Science and Technology had approved this research.

Peripheral blood was centrifuged to obtain blood cells, which were treated with erythrocyte lysis buffer (TBD, China) to obtain leukocytes for the assay. The total RNA of leukocytes was extracted with RNAiso plus kit (Takara, Japan), followed by reverse transcription to cDNA using cDNA RT-PCR kit (Takara, Japan). qRT-PCR was performed using a CFX Connect Real-Time System (Bio-Rad, USA) with SYBR Premix Ex Taq (Takara, Japan) and the specific primers. The primers (Sangon Biotech, China) were shown in Table 1 :

The NetworkAnalyst online tool was used to construct a gene regulation network of the hub gene–transcription factor interaction network. The transcription factor data were obtained from the JASPAR database to build the network (22, 23). Subsequently, the results were visualized using Cytoscape.

All data was analyzed by R (version 4.1.3) and GraphPad Prism 8 Software. For normally and non-normally distributed data, we used an unpaired t-test and a non-parametric test to compare differences between groups, respectively. Pearson correlation analysis was used to estimate the correlation. Unless explicitly stated, the p-value less than 0.05 was considered statistically significant.

The workflow of this study was presented in Figure 1 . A quality control procedure was performed on the scRNA-seq dataset. As shown in Figure 2A , specific low-quality cells were excluded based on their gene expression profile to avoid misinterpretation of the data. The mitochondrial proportion was used as a threshold to filter out low-quality cells. The CCA algorithm integrated the datasets from four COPD patients and three controls ( Figure 2B ). According to the expression level of marker genes ( Figure 2C ), all cells were classified as 16 clusters and annotated into 14 cell types based on the marker genes: T cells (marked with CD3D, CD2, and CD69), neutrophils (FCN1 and S100A9), dendritic cells (TGFBI and HLA-DMB), mast cells (KIT and TPSAB1), NK cells (CCL5, NKG7, GNLY, and CD247), B cells (CD79A, CD79B, and IGHG3), macrophages (MARCO, MSR1, and MRC1), alveolar type 1 cells (AGER, KRT7, and TSPAN13), alveolar type 2 cells (SFTPD, SFTA2, SLC34A2, and ABCA3), club epithelial cells (BIRC5, KRT19, and CXCL17), ciliated epithelial cells (CAPS, RSPH1, PIFO, and TSPAN1), endothelial cells (CALCRL and RAMP2), fibroblasts (COL1A2 and COL6A2), and smooth muscle cells (CALD1, TAGLN, and NOTCH3)( Figure 2D ).

According to the hallmark gene sets collected from the MsigDB database, we used the irGSEA package to calculate each cell’s gene set enrichment score. Overall, there were two hallmark gene sets upregulated in COPD samples: INFLAMMATORY-RESPONSE and TNFA-SIGNALING-VIA-NFKB. At the same time, the two gene sets were significantly downregulated for control samples ( Figure 3A ). Next, the same operation was on different types of cells. The two gene sets were both upregulated in neutrophils and DCs ( Figure 3B ). INFLAMMATORY_RESPONSE gene set was significantly enriched in neutrophils ( Figure 3C and Figure 3D ).

Neutrophilic inflammation predominates in the pathogenesis of COPD. The FindMarkers method was employed to identify genes differentially expressed in neutrophils from COPD patients and controls, revealing many inflammatory and immune responses-associated biological processes in COPD ( Supplementary Figure 1 ). Thus, we performed a simulation analysis on the cell trajectory differentiation of all neutrophils. It was shown that the darker the blue color, the earlier the cells differentiated, indicating that the neutrophils differentiate from the left to the right with the differentiation of time, and the latest differentiated cell was represented by the lightest blue ( Figure 4A ). Neutrophils labeled with different colors had three differentiated states, with the red one on the left being the earliest differentiated type ( Figure 4B ). The differentiation process between COPD and control neutrophils was then investigated, and the results showed that COPD neutrophils differentiated later, whereas control neutrophils differentiated earlier ( Figure 4C ). All the cells analyzed were neutrophils ( Figure 4D ). The GSEA was conducted in neutrophils of different differentiation states to investigate the relationships between pseudotime and gene functions and showed that the latest differentiated neutrophils were mainly associated with the upregulation of inflammatory and immune responses ( Figures 4E–G , Supplementary Figure 2 ).

The RNA sequencing dataset GSE57148 included 98 COPD and 91 normal samples. The cell type proportions in each sample were calculated by the CIBERSORT algorithm, and different colors represented different immune cells ( Figure 5A ). There were differences in the proportions of four types of immune cells, including neutrophils, M2-type macrophages, resting NK cells, and T follicular helper cells, in COPD samples and normal samples. Furthermore, the proportion of neutrophils was significantly increased in COPD samples ( Figure 5B ). In Figure 5C , the correlation heatmap indicated that neutrophils correlated negatively with Dendritic cells resting (p<0.01, r= -0.217), CD8 T cells (p<0.01, r= -0.216), T follicular helper cells (p<0.01, r= -0.305), NK cells activated (p<0.01, r= -0.191), and M2 type macrophages (p<0.01, r= -0.246). Conversely, a positive correlation between neutrophils and three cells was demonstrated, including Dendritic cells active (p<0.01, r= 0.236), monocytes (p< 0.01, r= 0.251), and Mast cells activated (p<0.01, r= 0.289).

As shown in Figure 6A , the clustering tree for samples screened out three abnormal samples, and 186 samples (96 COPD and 90 normal samples) remained for further analysis. Sequentially, based on the scale-free network, the best soft-threshold power of 14 was chosen. ( Figure 6B ). These genes were clustered into 12 modules based on their similar co-expression characteristics ( Figure 6C ). The correlation analysis of modules and traits (neutrophils and COPD) illustrated that the orangered4 module and lightssteelblue1 module were significantly associated with neutrophils (p<0.001, r= 0.5; p<0.001, r= -0.48) and COPD (p<0.001, r= 0.31; p<0.001, r= -0.46) ( Figure 6D ). Therefore, the genes in the orangered4 module and lightssteelblue1 module were selected as essential genes relevant to neutrophils for further analysis.

827 DEGs were identified between COPD and normal samples in GSE57148, including 390 upregulated and 437 downregulated genes ( Supplementary Figure 3 ). Venn diagrams were performed to obtain the DEG profiles intersection of scRNA-seq DEGs, RNA-seq, and WGCNA model genes, and 33 upregulated and 18 downregulated feature genes in total were remarkably expressed differentially in all 3 groups ( Figures 7A, D , Supplementary Table 1 ). KEGG pathway analysis revealed that the upregulated DEGs were mainly enriched in the IL-17 signaling pathway, NF-kappa B signaling pathway, and TNF signaling pathway. Meanwhile, the downregulated DEGs were primarily enriched in oxidative and collecting duct acid secretion ( Figures 7B, E ). The GO analysis showed that the upregulated DEGs were enriched in smooth muscle proliferation, secretory, and cytokine activity ( Figure 7C ). The downregulated DEGs were enriched in response to the metal ion, late endosome, and proton transmembrane transporter activity. ( Figure 7F ).

The MR analysis was used to explore the relationship between feature DEGs and COPD to verify whether the relationship was causal or correlative. Among these 51 feature DEGs, 9 genes were eliminated without suitable SNPs, and 2 genes had SNPs with P-value < 5 × 10^–6. Table 2 and Supplementary Table 2 displayed the causal effect of each feature DEG on COPD. IVW analysis reveals that NAMPT (OR=1.0102, p=0.0064), PTGS2(OR=1.0139, p=0.0161), SLC2A3(OR=1.0052, p=0.0014) and TRIB1(OR=1.0172, p=0.0405) levels were associated with an increased risk of COPD, while CDKN1A (OR=0.9962, p=0.0169), and CSRNP1(OR=0.9933, p=0.0384) were related to a decreased risk of COPD. Information on SNPs and other MR estimates for 6 genes was shown in Supplementary Table 3 . The Cochran’s Q test revealed no heterogeneity in these gene results. MR-Egger regression analysis indicated no horizontal pleiotropy in the results of these genes ( Table 2 ).

The protein–protein interaction (PPI) network of 51 feature DEGs was constructed based on the STRING online database ( Figure 8A ). The interaction networks were visualized in Cytoscape and contained 25 nodes and 47 edges. Furthermore, we implemented the cytoHubba plugin to search for the most significant hub genes by the maximal clique centrality algorithm ( Figure 8B ). The top 10 hub genes were identified, including IL1B, CXCL2, PTGS2, MCL1, SOD2, EGR1, CDKN1A, CXCL3, NAMPT, TNFAIP3, and ZFP36, and both NAMPT and MCL1 ranked 10th. The 11 hub genes were all upregulated genes. Correlation analysis revealed that these hub genes had a correlation with various immune cells, and all of them were significantly positively correlated with neutrophils (all P<0.05) ( Figures 8C–M ).

27 COPD patients and 27 healthy individuals provided peripheral blood samples. Since neutrophils make up about 70% of all leukocytes in peripheral blood, peripheral blood leukocytes were isolated from blood samples. The qRT-PCR in leukocytes was performed to verify hub genes. In Table 3 , the clinical characteristics of two group subjects were compared. Compared with normal subjects, the mRNA levels of IL1B, CXCL2, PTGS2, MCL1, SOD2, EGR1, NAMPT, CXCL3, and ZFP36 were markedly higher in COPD patients ( Figure 9 ). Moreover, correlation analysis revealed that IL1B, CXCL2, SOD2, NAMPT, and CXCL3 were inversely correlated with lung function in COPD patients (r = -0.4374, p<0.05; r = -0.5292, p<0.01; r = -0.4537, p<0.05; r = -0.5617, p<0.01; r = -0.5709, p<0.01; respectively) ( Figure 10 ).

The intersection of 6 genes identified by Mendelian randomization and 9 genes elevated in peripheral blood leukocytes of COPD patients identified 2 more crucial genes, NAMPT and PTGS2 ( Figure 11A ). Furthermore, NAMPT was negatively correlated with lung function in COPD patients. Using the JASPAR database, a network of transcription factors (TFs) for two genes was constructed. TF network showed that 10 TFs were obtained, of which 2 were GATA3 and FOXL1 with degree ≥ 2. ( Figure 11B ).