Single cell transcriptomics in blood of patients with chronic obstructive pulmonary disease

Abstract

Background

Chronic obstructive pulmonary disease (COPD) is a leading cause of morbidity and mortality worldwide. Single-cell RNA sequencing (scRNA-seq) provides gene expression profiles at the single-cell level. Hence, we evaluated gene expression in the peripheral blood of patients with COPD.

Methods

Peripheral blood samples from seven healthy controls and eight patients with COPD were obtained in this study. The 10X Genomics Chromium Instrument and cDNA synthesis kit were utilized to generate a barcoded cDNA library for single cell RNA-sequencing. We compared the scRNA-seq data between the COPD and control groups using computational analysis. Functional analyses were performed using Gene Ontology and Kyoto Encyclopedia of Genes and Genomes pathway analyses.

Results

scRNA-seq was used to analyze the transcriptome of peripheral blood mononuclear cells from seven normal controls and eight patients with COPD. We found an increased number of monocyte/macrophages in the COPD group compared to the normal control group. Among the differentially expressed genes (DEGs) in monocyte/macrophages, we identified 15 upregulated genes (EGR1, NR4A1, CCL3, CXCL8, PTGS2, CD83, BCL2A1, SGK1, IL1B, BTG2, NFKBIZ, DUSP2, MAFB, PLAUR and CCL3L1) and 7 downregulated genes (FOLR3, RPS4Y1, HLA-DRB5, NAMPT, CD52, TMEM176A and TMEM176B) in the COPD group compared to the normal control group.

Conclusions

Using scRNA-seq, we found differences in cell type distribution, especially in monocyte/ macrophages. Several upregulated and downregulated genes were found in the monocyte/macrophages of the COPD group.

Introduction

Chronic obstructive pulmonary disease (COPD) is a leading cause of morbidity and mortality worldwide. It is characterized by irreversible expiratory airflow limitation due to obstruction in the small airways and emphysematous destruction of distal airspaces in the terminal bronchioles [1]. Globally, COPD leads to a reduced quality of life and a high disease burden [2]. According to the World Health Organization, COPD affects more than 65 million people worldwide; in 2005, 5% of the global population (more than 3 million people) died from COPD-related illness, and COPD is expected to be the fourth most common cause of death by 2030 [3]. In South Korea, COPD is the seventh leading cause of death and is expected to result in higher medical costs, as well as social and economic burdens, due to continuous exposure to risk factors and an aging population [4]. Although many studies have been conducted to identify biomarkers for early detection and diagnosis of COPD, the biomarkers responsible for COPD remain poorly understood. To better understand how biomarkers can be used to diagnose COPD,, it is vital to identify gene expression-based molecular predictors using single-cell RNA sequencing (scRNA-seq) [5, 6].

Many genetic studies have been performed on lung tissue or bronchoalveolar lavage (BAL) fluid. However, because sample collection for these approaches is invasive, it is important to identify more easily accessible biomarkers. Therefore, we used peripheral blood samples to study transcriptomic profiles because of their ease of use in biomarker screening. As scRNA-seq can resolve cell expression variability within cellular populations, it could enhance our understanding of COPD biomarkers.

The development of scRNA-seq has greatly facilitated transcriptomic studies, leading to insightful findings regarding cellular expression variability and dynamics [7]. scRNA-seq provides unprecedented opportunities for exploring gene expression profiles at the single-cell level [6, 7]. These studies have revealed significant cell-to-cell gene expression variability. To discover blood biomarkers of COPD, we conducted scRNA-seq at the single-cell level on the peripheral blood of patients with COPD and normal controls. We analyzed differences in cell type distribution and differentially expressed genes (DEGs). Finally, we performed functional analysis using Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway analyses for COPD group compared to normal control group.

Method

Study design and population

Data from the Chronic Obstructive Pulmonary Disease in Dusty Areas (CODA) cohort were analyzed. The CODA cohort enrolled participants living in six cities (Gangneung, Donghae, Samcheok, Yeongwol, Danyang, and Jecheon) near cement plants in the Kangwon and Chungbuk provinces of South Korea between 2012 and 2017. Details of the design and materials of the CODA study have been described previously [8].

We obtained data from seven healthy participants and eight participants with COPD from the CODA cohort. Peripheral blood samples were obtained from patients living in two cities (Samcheok and Yeongwol), who were divided into two groups: COPD and normal control. Peripheral blood mononuclear cells (PBMCs) were isolated by density centrifugation and processed either directly or after frozen storage. The cells were then labeled with cell hashing antibodies and loaded onto droplet-based (10x) scRNA-seq platforms (Fig. 1).

Fig. 1

Overview of study design. Overview of study design for normal and COPD blood samples using scRNA seq

Written informed consent was obtained from all the participants. This study was approved by the Institutional Review Board (IRB) of Kangwon National University Hospital (IRB No. 2012-06-007-061).

scRNA-seq

We utilized the 10X Genomics Chromium Instrument and copy DNA (cDNA) synthesis kit (10x Genomics: Chromium Next GEM Automated Single Cell 5′ Library, Gel Bead Kit v1.1 and Human TCR/BCR amplification kit) to generate a barcoded cDNA library for scRNA-seq from approximately 3,000 cells. The quality of the cDNA library was assessed using an Agilent Bioanalyzer. Using this library, we performed two paired-end 100 bp Flow Cells using an Illumina Novaseq 6000 whole-genome sequencing system.

scRNA-seq analysis

Sequencing reads were processed with Cell Ranger v6.0.2 software, using the human reference genome, GRCh38. Feature-barcode matrices were generated, and downstream analyses were performed using the Seurat v4.2.0 package in R v4.0.3 software [9]. Quality control, cell filtering, normalization, and clustering were performed using Seurat. We annotated annotated cell type-specific markers using PanglaoDB [10] and confirmed them with a dot plot in Supplementary Fig. 1. Candidate marker genes with high expression levels in each cluster were selected and cell types were defined according to the criteria in Supplementary Table 1. Data integration was achieved through FindIntegrationAnchors and IntegrateData functions with default options. The integrated data was further processed using the ScaleData and RunPCA functions. We performed normalization and principal component analysis on the preprocessed matrix, using the first 16 principal components for clustering. The shared nearest neighbor algorithm was used to perform clustering with a resolution of 0.4, and visualization was performed using uniform manifold approximation and projection. To identify DEGs between the normal control and COPD groups, we used the default test in the “FindAllMarkers” function of Seurat (i.e., the Wilcoxon rank-sum test), which is the most widely used approach. scProportionTest v0.0.0.9000 was used to confirm differences between COPD and the normal control by cell type [11]. scProportionTest analyzes the differences between cell proportions corresponding to cell types between different groups. For DEGs, GO and KEGG pathway analyses were performed using the “ClusterProfiler v3.18.1” package in R, which supports statistical analysis and visualization of functional profiles for genes and gene clusters.

Results

Study characterization

The mean age was 78 years in the normal control group and 80.6 years in the COPD group (Table 1), and 73.3% of the participants were male. In the control group, there were four former smokers and three current smokers, while in the COPD group, there were five former smokers and one current smokers. The emphysema index was 7.8 in the control group and 6.44 in the COPD group. The FEV1/FVC value was normal at 0.75 in the normal control group but decreased to 0.58 in the COPD group. In the COPD group, there was one participant in GOLD stage 1, six participants in GOLD stage 2, and one participant in GOLD stage 3. There were no significant differences in age and sex between the normal control and COPD groups. The COPD group showed lower lung function compared to the control group. Using available clinical data, we compared the levels of peripheral blood cells (Fig. 2; Table 2).There were no significant differences in the blood cell differential count between peripheral blood samples from the control and COPD groups.

Fig. 2

Results of white blood counts in peripheral blood from COPD and normal groups. Definition of abbreviations: WBC white blood cell, LYMP lymphocyte, MONO monocyte, SEG segmented neutrophil, EOSINO eosinophil, BASO basophil

Table 1 Baseline characteristics of study participants

Table 2 Comparison of blood cell differential counts between the COPD and normal control groups

Cell subpopulation identification

We used Uniform manifold approximation and Projection (UMAP) for scRNA seq profiling. UMAP was used for reference-based cell type annotation in the COPD group. The UMAP visualization of scRNA-seq profiles (10x) included 56,097 cells from 15 samples (26,462 cells from seven normal control samples and 29,635 cells from eight COPD samples) (Fig. 3).

Fig. 3

UMAP visualization of 56,097 cells from 15 samples. 26,462 cells from seven normal control samples and 29,635 cells from eight COPD samples. Nine cell types were identified based on cluster-specific marker genes

We found an increased number of monocyte/macrophages in the COPD group compared to the normal control group (Fig. 4; Table 3). Among the cell types, only monocyte/macrophages showed a significant difference in the COPD group compared with the normal control group (FDR < 0.05) (Fig. 5).

Fig. 4

Comparison of cell ratio contained in each cell type in COPD and control groups

Fig. 5

Comparison of cell type significance between COPD and normal control groups. The average cell expression values of different COPD groups for the same cell type are compared, and the statistical significance of these comparisons is analyzed using advanced scProportionTest. Only monocyte/macrophages showed a significant difference in the COPD group compared with the normal control group (FDR < 0.05)

Table 3 Percentage of cells contained in each cell type in COPD and control groups. Cell subpopulation identification using scRNA sequencing of COPD and normal control groups

Analysis of DEGs

Bioinformatics analysis revealed the presence of DEGs in the COPD group. Volcano plots were used to visualize the identified DEGs in each cell type within the COPD group (Fig. 6). Based on the proportion of cell types, we focused on DEGs in monocyte/macrophages. For DEGs in other cell types, refer to Table 4.

Fig. 6

Volcano plot of monocyte/macrophage. |Fold change|>=1.2 and p-value < 0.05 as criteria for selecting DEGs. Red indicates upregulated genes in the COPD group compared to the normal control group. Blue indicates downregulated genes

Table 4 DEGs of monocyte/macrophage with |Fold change|>=1.2 and p-value < 0.05

Among the DEGs in monocyte/macrophages, we found 15 upregulated genes (EGR1, NR4A1, CCL3, CXCL8, PTGS2, CD83, BCL2A1, SGK1, IL1B, BTG2, NFKBIZ, DUSP2, MAFB, PLAUR and CCL3L1)and 7 downregulated genes (FOLR3, RPS4Y1, HLA-DRB5, NAMPT, CD52, TMEM176A and TMEM176B) in the COPD group compared to the normal control group (Fig. 6; Table 4).

Enrichment analysis of DEGs in the COPD group

Analysis of signaling pathways (GO, KEGG)

We found an increased ratio of monocyte/macrophages in the COPD group compared to the normal control group. Therefore, we performed further functional analysis of monocyte/macrophages.

GO analysis revealed that the upregulated DEGs in monocyte/macrophages of the COPD group were enriched in cytoplasmic translation, DNA replication, histone modification, ribonucleoprotein complex biogenesis, ribosome biogenesis, peptidyl-lysine modification, non-coding RNA processing, non-coding RNA metabolic processing, DNA-dependent DNA replication, and double-strand break repair. The downregulated DEGs in monocyte/macrophages of the COPD group were enriched in aerobic respiration, cellular respiration, T cell activation, adenosine triphosphate metabolic process, positive regulation of lymphocyte activation, generation of precursor metabolites and energy, adaptive immune response based on somatic recombination of immune receptors built from immunoglobulin superfamily domains, oxidative phosphorylation, immune response-regulating signaling pathways, and energy derivation by oxidation of organic compounds (Fig. 7).

Fig. 7

Enrichment analysis of monocyte/macrophage for Gene Ontology in the COPD and normal control groups. (Red for the COPD group, and blue for the normal control group)

The results of the KEGG pathway analysis revealed that the upregulated DEGs in monocyte/macrophages of the COPD group were enriched in the ribosome, cell cycle, nucleocytoplasmic transport, spliceosome, sphingolipid signaling pathway, nucleotide excision repair, and autophagy. The downregulated DEGs in the monocyte/macrophages of the COPD group were enriched in oxidative phosphorylation and proteasomes (Fig. 8). Both GO and KEGG pathway analyses indicated that the upregulated DEGs in monocyte/macrophages of the COPD group were enriched in ribosome biogenesis, while the downregulated DEGs in monocyte/macrophages of the COPD group were enriched in oxidative phosphorylation.

Differentially expressed genes were categorized and integrated into the transduction signaling map using the KEGG pathway database. KEGG pathway analysis was performed on monocyte/macrophages (Fig. 8).

Fig. 8

Enrichment analysis of monocyte/macrophage for KEGG in the COPD and normal control groups. (Red for the COPD group, and blue for the normal control group)

Discussion

In this study, we identified genes in each cell type of the COPD group compared to the control group through scRNA analysis, uncovering novel insights. Our study had two major findings. First, we found an increased ratio of monocytes/macrophages in the COPD group compared to those in the control group. Although this has been previously reported in studies on COPD [12], out findings provide new context and significance. Second, among the DEGs in monocytes/macrophages, we found 15 upregulated genes (EGR1, NR4A1, CCL3, CXCL8, PTGS2, CD83, BCL2A1, SGK1, IL1B, BTG2, NFKBIZ, DUSP2, MAFB, PLAUR and CCL3L1) and 7 downregulated genes (FOLR3, RPS4Y1, HLA-DRB5, NAMPT, CD52, TMEM176A and TMEM176B) in the COPD group compared to the control group.

Neutrophils, the most abundant cells in the sputum or BALF of patients with stable or worsening COPD, are associated with the severity of airflow obstruction in COPD [13]. Studies have shown that inhibiting these cells or their functions can improve symptoms and reduce lung damage in chronic respiratory diseases [14, 15]. Although not statistically significant, our study found increased numbers of neutrophils, eosinophils, and white blood cells, which is consistent with previous research. Neutrophils have low RNA content, which limits the number of transcripts that can be analyzed within these cells compared to other cell types. This necessitates a separate protocol or additional collection of neutrophil cells for in-depth analysis during the experiment [16, 17]. We should have explored this aspect in the current study, as our focus was on a comprehensive observation of the proportions of all cell types in the COPD group compared to the normal control group. This suggests that further studies are needed to fully understand the role of neutrophils in chronic respiratory diseases.

Several studies have used scRNA-seq to profile gene expression in patients with COPD. Xiuying et al. [18] reported that scRNA-seq is useful for identifying transcriptional changes and individual protein levels that may contribute to the development of emphysema in a cell type-specific manner. They performed scRNA-seq on peripheral lung parenchymal tissue to confirm cell-specific transcriptomic changes between normal and COPD groups, finding more differentially expressed genes in monocytes, macrophages, and ciliated epithelial cells in the COPD group compared to the control group. Notably, QKI and IGFBP5 protein levels were especially high in the ciliated epithelial cells of the COPD group. Thus, in the future, scRNA-seq profiling of gene expression could potentially be used as an additional tool for COPD diagnosis and severity assessment in the future. However, further investigations are needed before scRNA-seq profiling of gene expression can be established as a clinical biomarker for diagnosing COPD.

In this study, we found an increased ratio of monocyte/macrophages in the COPD group compared to the normal control group. The increased ratio of monocytes/macrophages in the COPD group compared to the control group highlights the importance of the inflammatory processes, including the role of monocytes/macrophages, in the pathogenesis of COPD. This is because monocytes/macrophages are key players in immune responses and inflammation [19]. Monocytes have been implicated in several diseases. For instance, a recent study suggested that peripheral blood monocyte count in patients with idiopathic pulmonary fibrosis might be a valuable predictor of poor prognosis [20]. Similarly, activated circulating monocytes have been observed in patients with systemic sclerosis, indicating their potential involvement in disease processes [21]. Based on these findings, we hypothesized that the increased presence of monocyte/macrophages in COPD would contribute to its pathogenesis and that the upregulated or downregulated genes in these cells could serve as significant biomarkers for the early diagnosis of COPD in the future.

Among the DEGs in monocyte/macrophages, we found 15 upregulated genes (EGR1, NR4A1, CCL3, CXCL8, PTGS2, CD83, BCL2A1, SGK1, IL1B, BTG2, NFKBIZ, DUSP2, MAFB, PLAUR and CCL3L1) and 7 downregulated genes (FOLR3, RPS4Y1, HLA-DRB5, NAMPT, CD52, TMEM176A and TMEM176B) in the COPD group compared to the normal control group.

The EGR1 (early growth response 1) gene encodes a member of the EGR family of C2H2-type zinc-finger proteins, which are nuclear proteins that function as transcriptional regulators. These proteins are essential for differentiation and mitogenesis and have been suggested to act as cancer suppressors in other studies. EGR1 plays a critical role in promoting autophagy and apoptosis in response to cigarette smoke exposure [22]. In pulmonary epithelial cells, EGR-1 is induced by cigarette smoke and contributes to proinflammatory mechanisms that are likely involved in the pathogenesis of COPD in the lungs of smokers [23]. These data are particularly interesting given the important role of chronic inflammation in COPD.

The CXCL8 (C-X-C motif chemokine ligand 8) gene encodes a protein that is a member of the CKC chemokine family and serves as a major mediator of the inflammatory response. Especially, CXCL8 is considered a key regulator of local bronchial inflammatory response in COPD [24]. Additionally, evidence suggests that CXCL8 is involved in peripheral muscle damage in myositis [25] and in the relationship between peripheral muscle strength and systemic levels of inflammatory markers in patients hospitalized for acute exacerbation of COPD [26]. This underscores the relevance of the inflammatory response in the pathogenesis of COPD.

In conclusion, the 15 upregulated genes (EGR1, NR4A1, CCL3, CXCL8, PTGS2, CD83, BCL2A1, SGK1, IL1B, BTG2, NFKBIZ, DUSP2, MAFB, PLAUR and CCL3L1) found in monocyte/macrophages in the COPD group are associated with proinflammatory functions related to COPD. These findings highlight the critical role of inflammatory responses in COPD and underscore the potential of these genes as biomarkers for diagnosing and assessing the severity of COPD.

Among the 7 downregulated genes (FOLR3, RPS4Y1, HLA-DRB5, NAMPT, CD52, TMEM176A and TMEM176B) in the COPD group compared to the normal control group, HLA-DQB1 stands out as particularly notable. Its function is described below.

The HLA-DRB5 (major histocompatibility complex, class II, DR beta 5) is part of the human leukocyte antigen (HLA) class II beta-chain paralog. It plays a central role in the immune system by presenting peptides derived from extracellular proteins. HLA-DRB5 is an expression quantitative trait loci (eQTL) expressed in neutrophils and may be involved in immune recognition and regulation. It has also been implicated in mitochondrial function and associated with diseases such as inflammatory bowel disease, ulcerative colitis, Type 1 diabetes, multiple sclerosis, and systemic lupus erythematosus [27]. In another study, HLA-DRB5 was identified as hypermethylated in patients with smoking-related COPD [28]. Future studies are required to correlate DNA methylation status with the transcriptomic profiles of the genes identified in this study and to elucidate their relationship with COPD.

Functional analysis using GO and KEGG pathway analyses revealed that the upregulated DEGs in monocyte/macrophages of the COPD group were enriched in ribosome biogenesis. This finding is consistent with another study that also identified DEGs related to COPD as being enriched in ribosomal biogenesis through GO enrichment analysis [29]. Additionally, both GO and KEGG pathway analyses in our study showed that the downregulated DEGs in monocyte/macrophages of the COPD group were enriched in oxidative phosphorylation.

Finally, this study utilized peripheral blood samples, which are more easily obtained than lung tissue or BAL samples, highlighting a key advantage of our approach. Unlike studies that primarily use tissue samples for scRNA-seq, the use of blood samples offers several benefits: it is less invasive, more cost-effective, and has the potential for developing useful blood biomarkers.

This study had several limitations. First, the number of participants was limited, with only eight patients with COPD and seven healthy controls included in the analysis. Further studies with larger sample sizes could provide more robust results and identify reliable biomarkers for COPD development. Nonetheless, this study made significant contributions by analyzing cell populations and gene expression at the single-cell level in the peripheral blood of COPD patients, revealing important biological pathways. Second, our study found that the monocyte/macrophage group displayed the most significant difference between the control and COPD groups. Since macrophages are known to arise from circulating blood monocytes, these two cell types share many markers in common. Nevertheless, the monocyte/macrophage population can be re-clustered, but may display distinct proportion of subtypes within the population depending on the stage of COPD [30]. Thus, large numbers of COPD tissue samples from the different stages of COPD patients are needed to obtain an exact trend of subtypes after monocyte/macrophage re-clustering. Third, two key factors influence the appearance of genes during DEG analysis. Firstly, when the number of cells in the two comparison groups is small, comparing genes poses a significant challenge, resulting in fewer genes appearing accurately. Secondly, in the case of low-content RNA, i.e., cells with low gene expression values, the number of genes may vary due to their identification as non-cells in the initial filtering stage and subsequent removal. This study did not incorporate cell number correction; instead, it relied on the average gene expression value of all identified cells in the corresponding cell type. In future studies, we aim to determine the number of cells and the cell inclusion ratio by cell type to enhance the accuracy of our analyses. Additionally, scRNA-seq only recovers a few thousand unique transcripts from a single cell, which is far less than would be recovered from a complete transcript profile, thus limiting its range of applications. However, this study revealed meaningful cell-to-cell gene expression variability.

In conclusion, we found differences in cell type distribution, especially an increase in monocyte/macrophages, in COPD group compared to normal control group, using scRNA-seq on peripheral blood samples. These findings suggest that monocyte/macrophages may be associated with the mechanisms underlying COPD development. Further research is needed to elucidate the precise role of monocyte/macrophages in the development of COPD.

The five upregulated and seven downregulated genes identified in monocyte/macrophages of the COPD group have potential as biomarkers for evaluating COPD diagnosis and severity. Further studies are needed to fully understand the roles of the these genes in COPD and validate their utility in clinical settings.