Transcriptomics Analysis Identifies the Presence of Upregulated Ribosomal Housekeeping Genes in the Alveolar Macrophages of Patients with Smoking-Induced Chronic Obstructive Pulmonary Disease

Introduction

The prevalence of chronic obstructive pulmonary disease (COPD) is increasing worldwide, resulting in a higher mortality in patients with this disease. The global prevalence of COPD has been found to be 10.1%. In 2017 alone, this represented 3.2 million deaths.¹ It has become the third largest cause of death, and the fifth largest global medical burden.^1–3 Long-term smoking, from the consumption of tobacco, and biomass smoke, is the most common and major cause of COPD. As a heterogeneous and complex airway disease causing chronic inflammation, the pathogenesis of COPD is characterized by the irreversible limitation of airflow, accompanied with hypersecretion, remodeling and an impaired immune response of the mucous membrane. Proteolytic enzymes and reactive oxygen species contribute to tissue damage, and results in other debilitating conditions such as emphysema. The available data shows that only 15–20% of smokers develop severe COPD, but the mechanisms behind the development of COPD are yet to be elucidated.^3–5

Alveolar macrophages (AM) participate in the phagocytosis and removal of foreign particles or pathogens and plays a vital role in maintaining alveolar structure and function, as well as the non-specific immunity in lung tissue. The impact of smoking on AM has been reported in several studies.^6–9 The transcriptomic alterations of AM were, however, only investigated in a few small sample studies.^5,10,11 Although COPD has been well known as chronic inflammation, the exact pathways or cytokines implied in the process of COPD are still not well understood. Thus, using bioinformatics we retrieved and sorted the datasets of similar studies available in the GEO database in the last 15 years, explored the transcriptomic influence of cigarette smoking exposure on COPD patients and healthy people, seeking the hub genes. Providing the greatest correlation in the relevant gene module or interaction network, hub genes are often key players in a number of clinical traits. The results will provide a greater insight and increased understanding of the pathogenesis and biological features associated with smoking-induced COPD at the transcription level.

Materials and Methods

Keywords such as “COPD”, “smoking”, “cigarette”, “tobacco”, “macrophage” were used as terms for searching the GEO database (https://www.ncbi.nlm.nih.gov/geo/), and the species was limited to “homo sapiens”, preferring to select the same popular GPL platform, in accordance with case-control trials, until August 2020. A total of 4 datasets: GSE2125 (n=45), GSE8823 (n=24), GSE13896 (n=88), GSE130928 (n=70), representing a total of 203 cases of GSM, and all on an Affymetrix HG-U133 Plus 2.0 (GPL570 microarray), were obtained (Figure 1).

Figure 1 Flow Chart of Data Processing for 4 GSEs.

The total of 188 cases obtained (for smoking capacity calculation) were comprised of 63 non-smokers (22 female, 41 male) with an average age of 40.35±7.95 years, 91 non-COPD smokers (23 female, 68 male) with an average age of 43.68±8.56 years, and 34 COPD smokers (6 female, 28 male) with an average age of 54.25±8.35 years. A difference in smoking history existed within the COPD and non-COPD smokers, where non-COPD smokers consumed 26.91±17.20 pack-years (n=76), which was less than that of the COPD smokers, who consumed 45.40±27.87 pack-years (n=34, P<0.01), which is determined using the t-test, with the exception of 15 smokers in GSE2125 that had no information concerning their smoking history. AM were collected by bronchoalveolar lavage (BAL), as previously described.¹²

One hundred and twenty-five cases of GSM were collated as the smoking population, and the healthy smokers (non-COPD, n=91) were compared with the COPD smokers (n=34). After data cleaning using R 3.6.1 software, irrelevant or duplicate data were excluded (some data within the GSE13896 set came from GSE8823, and GSE2125 included some asthmatic data) and 54,675 probes were obtained in total. The combined expression matrix, grouping information and GPL annotation information were also extracted. Following the process of the batch effect (Figure 2), analysis and visualization of the differentially expressed genes (DEGs), including a PCA (principal component analysis) chart, heat-map and volcano plot were performed using “means” of the corresponding R package “limma”; The above data then underwent enrichment analysis by Gene Ontology (GO) and Gene Set Enrichment Analysis (GSEA). The interacting genes were retrieved by STRING platform (www.string-db.org/) based on DEGs and the protein–protein interaction (PPI), which was constructed to screen for hub genes using Cytoscape 3.8.1 software. A gene possessing high degrees in a PPI network would be defined as a hub gene in the module. The hub genes of interest were also identified and analyzed in comparison using the t-test. In addition, 155 cases of GSM were collated as the non-COPD population, the healthy control (non-smokers, n=64) versus the healthy smokers (non-COPD, n=91). The dataset was processed using the same analysis process described above (Figure 1).

Figure 2 PCA of 4 GSEs before and after Batch Effect Adjustment. The 4 datasets unprocessed (against Batch Effect) are shown divergently by means of PCA (A); After the adjustment allowing for the Batch Effect, the datasets are integrated by the COPD smokers and the non-COPD smokers and the inherent difference can be sharpened (B). The colored dots represent the probes corresponding to relevant genes symbolically and quite a number of genes are differently distributed between the groups (B). The identical process of PCA was applied to the non-COPD groups with the unprocessed data shown (C). As the normalized data is displayed, there is only a limited but definite difference of gene expression between the non-smokers and the non-COPD smokers (D).

R package “limma”, which is usually used for the analysis of a gene expression microarray or sequencing data using a generalized linear model and empirical Bayes method, can fit a linear model to the expression data for each gene and “retrain” the false positives, making the overall analysis stable.¹³ GO enrichment analysis is designed to calculate the P value of each pathway (predefined according to GO) with a differential gene expression, which is part of an over-representation analysis based on hypergeometric distribution.¹⁴ As a computational method, GESA is used to determine whether a priori defined set of genes (obtained from KEGG) shows statistically significant, concordant differences between two biological states. GSEA first ranks all the genes according to their various expression by means of functional class scoring (based on a Kolmogorov–Smirnov test). The ranked list is then used to assess how the genes within each gene set are distributed across the entire ranked list.¹⁵ Adjusted P<0.05 and |log FC| (fold change) > 1.5 were considered to be statistically significant.

Results

Principal Component Analysis (PCA)

A PCA highlighted that there was a significant difference between the genes expressed by the COPD smokers and non-COPD smokers (Figure 2B), whereas there was only a limited difference between the genes expressed by the non-smokers and non-COPD smokers (Figure 2D).

Differential Gene Expression

When comparing the COPD smokers and non-COPD smokers, the top ten significant gene differences (sorted by adj.P<0.05) were: NADK, DRAP1, DEDD, NONO, KLHL12, PRKAR1A, ITGAL, GLE1, SLC8A1 and SVIL (Table 1). Most of these are down-regulated genes, with KPNA2 being the only up-regulated genes (in the top 20, not shown). When comparing the non-COPD smokers and non-smokers, the top ten genes with significant differences (sorted by adj.P<0.05) were: VGLL3, SLC26A11, NOL4L, ATP6V0D2, ITGAE, CADPS2, LPAR3, PDE4B, IFI27 and SEMA4D (Table 1). The five most distinct down-regulated genes were VGLL3, NOL4L, CADPS2, LPAR3 and PDE4B, while the three most significant up-regulated genes were SLC26A11, ATP6V0D2 and ITGAE. With the default logFC (0.585) and adj.P (<0.05), the differently expressed genes have been visualized in the heatmap and volcano plot presented in Figure 3.

Table 1 The Top 10 Genes with Most Significant Difference of Transcription Between Relevant Groups

Figure 3 The Visualization of Differently Expressed Genes Between Relevant Groups. According to logFC, the Heatmap (A) displays the top 40 of differential genes expressed (both upregulated and downregulated) between the COPD smokers and the non-COPD Smokers. With the default logFC (0.585) and P value (0.05), the volcano plot (B) shows the top 20 differential expression of relevant genes. Likewise, the Heatmap (C) displays the top 40 differential genes upregulated or downregulated between the non-COPD smokers and the non-Smokers. With the same default of logFC and P value, the volcano plot (D) shows the top 20 differential expression of relevant genes.

The Enrichment Analysis

For the COPD and non-COPD smokers, a GO enrichment analysis shows that there is a remarkable enrichment of cell components, biological processes and molecular functions. Cellular components include intercellular and cytoplasmic adhesion components, lysosomal membranes and other structures; Biological processes include neutrophil activation, migration, degranulation and related immune processes. Molecular functions consist mainly of binding entities such as cell adhesion molecules (CAM) including cadherin and others (Table 2, Figure 4A–C).

Table 2 The Top 10 of Most Significant Pathways Between Relevant Groups

Figure 4 GO Enrichment Analysis of AM Between Relevant Groups. (A–C) GO comparison of the COPD Smokers and the non-COPD Smokers; (A) for Cell Components (Top 10), (B) for Biological Process (Top 10), (C) for Molecular Function (Top 10). (D–F) GO comparison of the non-COPD Smokers and the non-Smokers; (D) for Cell Components (Top 10), (E) for Biological Process (Top 10), (F) for Molecular Function (Top 10).

It was also found that more than 120 enriched pathways were acquired by GSEA (P<0.05, minGeneset Size: 20, maxGeneset size: 200). Ribosome and RNA transport are the only two up-regulated pathways (red-colored in Figure 5A) in contrast with the majority of down-regulated pathways. The latter includes processes such as the inflammation signaling pathways of various viral, bacterial and even protozoal infections, pinocytosis processes, cytoskeleton and connection structures, cellular adhesion molecules (CAMs) and chemokine signaling pathways, pattern recognition receptor (PRR) signaling including NLR and TLR, NF-κB pathways; PI3K-Akt pathways, mTOR signaling; phospholipase D signaling, TNF signaling; lysosome and phagosome, cell senescence and necroptosis, osteoclast differentiation; a variety of carcinomas and related processes including relevant misregulating transcription and microRNAs, viral carcinogenesis, PD-L1 and PD-1 pathways, HER2 signaling, tumor proteoglycans; EGFR and VEGF pathways; insulin-related signaling including relaxin and FoxO signaling as well as RAGE-AGE signaling, Type II diabetes processes and carbohydrate metabolism; adrenergic signaling, hypertrophic and dilated cardiomyopathy; JAK-STAT signaling and hematopoietic lineage; T and NK lymphocyte activation and differentiation, BCR and FcεRI signaling; leukocyte migration, platelet activation; atherosclerosis process; adipocyte metabolism; neurotrophin signaling and cholinergic synapses; prolactin, growth hormone, parathyroid and thyroxine processes; Hippo signaling, sphingolipid signaling, HIF-1 signaling, Notch signaling, Hedgehog signaling, etc. As in Figure 5A, only pathways with P<0.0005 are displayed for the figure size.

Figure 5 GSEA of AM Between Relevant Groups. (A) GESA for the non-COPD Smokers vs COPD Smokers. The blue-colored part on the left represents the relatively down-regulated pathways according to the enrichment analysis by GSEA, while a few of the red-colored parts represent the relatively up-regulated pathways. The default cutoff of P value is 0.05, with a total of more than 120 pathways being achieved, thus only pathways with P value <0.005 are shown here. (B) GSEA for the non-Smokers vs non-COPD Smokers. The blue-colored part on the left represents the relatively down-regulated pathways according to the enrichment analysis by GSEA, while a few of the red-colored part represent the relatively up-regulated pathways. The default cutoff of P value is 0.05, with a total of 50 pathways being achieved. The trend of up- or down- regulated is based on the algorithm of functional class scoring (FCS).

For the non-COPD smokers and non-smokers (ie, the healthy control), the GO enrichment analysis highlighted evident enrichment on molecular function (MF) and biological processes (BP). The cytokine activity and its receptor binding pathways were found to be the most significant processes for molecular functions. As for biological processes, leukocyte migration, inflammatory cells releasing chemokines, and non-specific immune responses to pathogen components, such as LPS are the most prevalent (Table 2, Figure 4D–F). Likewise, using GSEA, 50 enrichment pathways were acquired (P<0.05, minGeneset Size: 20, maxGeneset size: 200) with quite a proportion of them being analogous to the above, include the relevant pathways or mechanisms of viral or bacterial or protozoal infections, T lymphocyte differentiation and activation pathways, and other various steps associated with inflammatory processes (Figure 5B).

Protein–Protein Interaction (PPI) Analysis

The PPI analysis of genetic differential expression between the COPD and non-COPD smokers highlighted that ribosomal protein-coding genes, such as RPS11, RPL14, RPL27A, RPS27L, RPL37A and RPL38, are the hub genes, which are also related to PPP2R1A, EIF4B, TUBB, SMURF2, ACTB, and are closely connected to each in a network (Figure 6A). All of the ribosomal genes (above) in COPD smokers have been remarkably lifted, as shown in Figure 6C, with the most significant being RPL38. The result is in accordance with the up-regulated tendency in GSEA.

Figure 6 Protein–Protein Interaction Network and Hub gene Expression. All of the nodes could be regarded as hub-genes of the relevant network in (A and B). (A) PPI comparison of the Cytoscape between the COPD and non-COPD Smokers. The Hub genes have been displayed as the color nodes. The most deep-colored are RPL27A, RPL14, RPL37A, RPL38, RPS11 and RPS27L, signifying the most concentrated genes of the network. (B) PPI comparison of the non-COPD smokers and non-smokers. The most deep-colored are CCL5, CCR5, CXCL9 and CXCL11, signifying the most concentrated genes of the network. (C) The expression of ribosome-encoding genes such as RPS11, RPL14, RPL27A, RPS27L, RPL37A and RPL38 between the COPD and non-COPD Smokers. All of these genes in COPD Smokers have been remarkably up-regulated, with RPL38 being most significant, with a P value of 3.2×10−7.

As for the non-COPD smokers and the non-smokers, the PPI analysis highlighted that the hub genes present were dominantly CCL5, CCR5, CXCL9 and CXCL11 and other chemokines (or receptors) of the network, which are related to ITGAE, LPAR3, EMR1, GBP2/5 and other G protein coupling or related components, as well as CD69/80 and other T lymphocyte activation receptors and TNF superfamily members. Out of these, the most prominent were the aforementioned chemokines (Figure 6B).

Discussion

Although it has been well-recognized that cigarette smoke is the most common cause of COPD, the specific mechanisms responsible have not yet been elucidated.^16,17 This study analyzed the last 15 years of datasets (until August 2020) available from the GEO database, and found that there was an extensive and significant degree of gene differential expression of AM between the COPD smoking and non-COPD smoking group. It was found that people with COPD have a number of significantly differently expressed genes related to energy metabolism or oxidative stress than those without COPD. For example, NADK catalyzes the phosphorylation of nicotinamide adenine dinucleotide (NAD) into NADP. The latter is widely considered as an electron donor in the process of biosynthesis, especially the most basic oxidative phosphorylation of cells and oxidative stress. SLC8A1 encoding protein NCX can affect the process of oxidative phosphorylation by regulating the Ca²⁺ concentration in mitochondria.¹⁸

GSEA manifests a definite trend of up-regulation of the ribosomal pathway as well as a substantial down-regulated trend of most other pathways, which is coincident with the transcriptional down-regulation of massive genes. A PPI analysis indicated that quite a few hub genes are encoding proteins of ribosomal subunits, also known as housekeeping genes. It was also found that the expression of ribosomal genes in COPD smokers were all significantly elevated (Figure 6C, sorted by adj.P<0.05), which suggests that the process of protein synthesis has been affected. Nevertheless, the results obtained using the GO enrichment analysis were found to be dissimilar to the results obtained using GSEA. GO enrichment analysis is based on differential gene expression, whereas GSEA is designed to be not limited in the DEGs and could reveal whether the accumulation of subtle changes in a priori defined set of genes could be statistically significant and reflect the difference between the two biological states or phenotypes.

Studies have found that the AMs in smokers lose a certain degree of their pro-inflammatory properties, such as phagocytosis, their release of inflammatory mediators, and reduction of NO, which highlights the imbalance between the so-called M1 pro-inflammatory phenotype and the variant M2 anti-inflammatory phenotype.^7,8 This imbalance is in accordance with the results of this study, which highlights the presence of vast numbers of down-regulated pathways. This imbalance not only raises the instances of mortality through certain pathogenic infections, but also enhances the occurrence and development of tumors.⁴ Several transcriptomics studies of AM have focused on COPD.^5,11,19 An earlier study found that the expression of osteopontin (encoded by the gene SPP1) was up-regulated (4 times, P=0.006) in the AMs of healthy smokers, with significant differences in CCR5, CD69 and ITIH5. This also suggest that AM is activated by exposure to smoking, however the sample size was limited.¹⁰ Another study found that MERTK, the receptor responsible for recognizing apoptotic bodies on AM surfaces of healthy asymptomatic smokers, was significantly augmented, even more in patients with mild or moderate COPD.¹¹ Using an RNA sequencing method, Morrow et al showed the enrichment of amyloid and telomere-related pathways across the three tissues (AM, peripheral blood and airway epithelium) in a study comparing 11 COPD smokers and 6 non-COPD smokers. In this study, the SCGB1A1 (secretoglobin family 1A member 1), which encodes for CC16 (Club Cell Secretory Protein), a blood biomarker of COPD, was differently expressed,²⁰ which as well implied somewhat systemic effect of COPD.

Shaykhiev et al obtained a different result, whereby M1 inflammation-related chemokines, such as CXCL9, CXCL10, CXCL11 and CCL5 (obtained by analysis of the AMs of healthy asymptomatic smokers), were relatively down-regulated, while M2-related species such as MMP2, MMP9 and ADORA3 were up-regulated.⁵ A later study, which adopted a multi-factor mixed method to stimulate AM (multiple cytokines and different combinations), not only confirmed the existence of the M1 or M2 polarized inflammation phenotypes, but illustrated there were even up to 9 different expression patterns of the inflammatory response existing, highlighting the complexity of chronic inflammation.²¹ Not only an increased number of AM of patients with COPD, but the plasticity of AM to change phenotypes, have been already discovered. Different phenotypes and the distribution of AM largely depend on their microenvironment, including various inflammatory mediators or cytokines and even the loci of small airways.²² In another research, it was found that AMs in patients exposed to smoking could simultaneously express both markers of M1 (iNOS) and M2 (arginase), while dual positive (both M1 and M2 marker-positive) macrophages were more frequently distributed in lung tissues of subjects with severe COPD. These results allow the postulation that AM populations with M1 and M2 signatures do not exclude each other, but can often coexist.²³

Our study further compared the transcriptomic profile of AM between the non-smokers (the healthy control) and the non-COPD smokers (the healthy smokers). A PPI analysis indicated that CCL5, CCR5, CXCL9 and CXCL11 were the hub genes. For COPD patients with a smoking history, the levels of chemokines, such as IL-6, CCL26 of blood and IL8, CXCL1, CCL5 of AM had increased significantly, which could also be correlated with the degree of airflow limitation.^24–26 The results obtained through the PPI and enrichment analysis performed in this study suggests that related chemokines could play a crucial role, which is consistent with previous studies and can therefore be regarded as a verification of the bioinformatic analysis contained within this paper.

The main challenge with this work is that the relevant differential genes or hub genes determined as being significant in this study have not yet been verified in the laboratory. The sample size of the smoking COPD patients was relatively small compared to that of the non-smoking COPD patients (34:91) and this may have lead to a statistical bias, however the results obtained are consistent with a number of other similar studies. Though ICS (inhaled corticosteroid) usage has an impact to alter AM phenotype and function,²⁷ there is no relevant information provided within the data set.

In conclusion, it is clear that continuous exposure to cigarette smoking could trigger the onset of COPD. This study has revealed the comprehensive trend of repressed pathways and elevated ribosome proteins, each associated with the development of the disease. A large number of hub genes are up-regulated ribosomal genes, suggesting that exposure to cigarettes could adversely affect the synthesis of relevant proteins, and hence increase the subsequent pathogenic risks.