Cell-type-specific gene expression changes under different conditions
We performed scRNA-seq on lung tissue samples from 3 COPD patients (62 ± 11.53) and 6 controls with varied histories of smoking (4 active smokers and 5 never-smokers) (Fig. 1, Additional file 3: Dataset 1, and see Additional file 1: Supplementary Methods). Among controls, 3 were age-matched to COPD patients (72 ± 2.00), while the rest were from young individuals (30 ± 4.36). Hence, by deconvolving molecular complexity with a linear mixed model, we were able to characterize the molecular pathogenesis with the roles of aging and smoking assessed independently (see Additional file 1: Supplementary Methods). Data of scRNA-seq are notorious for false-positive and trivial findings as a result of the large amount of data. Therefore, to identify the most prominent alteration not only at the level of individual cell types but also in the whole system, we adopted the following analytic strategy. Starting with cell-type prioritization analysis, we first identified the most affected cell types under each condition, including COPD, aging, and smoking. Subsequently, we explored causality between the prioritized cell types and others with an array of advanced bioinformatic tools, such as CCIs analysis, integrating cell-type-level alterations into a system-level malfunction. The system-level malfunction occurring under the COPD condition was defined as the “COPD core” pathology. Ultimately, we combined findings under all conditions and established a new comprehensive COPD pathological model in which how aging and smoking facilitate COPD development was also clarified (Fig. 1).
Upon quality control (Additional file 2: Fig. S1B–S1D) and data processing (see Additional file 1: Supplementary Methods), we classified cells into 15 cell types, including 14 known cell types (macrophages, dendritic cells (DCs), monocytes, mast cells, neutrophils, natural killer cells (NKCs), T cells, B cells, alveolar type 1 cells (AT1s), alveolar type 2 cells (AT2s), club cells, ciliated cells, stromal cells, and endothelial cells) and a distinct unknown cluster highly positive for cell proliferation markers named proliferating cells, which was also detected in a recently published scRNA-seq dataset of human lung tissues [23] (Fig. 2A, Additional file 2: Fig. S2A and S2B, and see Additional file 1: Supplementary Methods). In general, cell compositions were not biased towards any level of conditions (Additional file 2: Fig. S2C–S2E).
To independently assess molecular alterations associated with COPD, aging, and smoking in a cell-type-specific manner, we adopted a linear mixed model to identify differentially expressed genes (DEGs) under each condition (see Additional file 1: Supplementary Methods, and Additional file 4: Dataset 2). Expectedly, aging-associated DEGs were enriched with molecular hallmarks of aging [10], reflecting the general status of aged lungs from elderly subjects; the biological ontologies associated with smoking appeared to be the most heterogeneous, yet they still highlighted cigarette-smoking-induced cellular stress and changes in lipid metabolism and vesicle transport (Fig. 2B and Additional file 2: Fig. S2F). Notably, regulation of response to cytokine stimulus, cytokine production, and innate immune response were among the top dysregulated Gene Ontology (GO) terms associated with COPD, emphasizing the involvement of cytokines and the innate immune system in COPD pathogenesis. To further reveal genes responsible for COPD in particular, we showed that the top COPD-associated DEGs among immune cells were from monocytes, including nuclear factor-kappa B (NF-κB) signaling genes such as REL, IL1B, KLF4 and FLNA; among non-immune cells, the top DEGs were from AT2s and club cells, including AZGP1, SNRPG, HNRNPM, EIF3E, TIE1, and IL7R (Fig. 2C). We also integrated genome-wide association study (GWAS) data [24] with COPD-associated genes to determine cell types contributing to COPD heritability and demonstrated that the GWAS COPD risk genes were mostly overrepresented among DEGs in monocytes, followed by AT2s and endothelial cells (Fig. 2D and see Additional file 1: Supplementary Methods).
To determine the prioritized cell types under COPD, aging, and smoking, we downsampled the data to compare the same number of cells across all cell types (see Additional file 1: Supplementary Methods) and found that monocytes had the largest number of DEGs in COPD, while club cells and macrophages were the most affected types under aging and smoking, respectively (Fig. 2E and Additional file 2: Fig. S2G). This finding was validated by AUGUR, which prioritizes cell types by quantifying the separability of perturbed and unperturbed cells within a high-dimensional space (Fig. 2F and see Additional file 1: Supplementary Methods). Note that both analyses revealed that monocytes were among the top-prioritized cell types under all conditions; therefore, they were the focus of our following analyses.
Condition-specific status of monocytes presents snapshots of system-level malfunctions
Given the shared functions and transcriptomic similarity among cell types in the mononuclear phagocyte system (MPS), we first ruled out the possibility of clustering errors by comparing our data with 2 recently published lung scRNA-seq datasets [25, 26]. Our in-depth classification of mononuclear phagocytes was in high agreement with previous annotations (Additional file 2: Fig. S3A and see Additional file 1: Supplementary Methods), including 5 subsets of DCs (myeloid DC1s (mDC1s), myeloid DC2s (mDC2s), plasmacytoid DCs (pDCs), IGSF2 + DCs, and TREM2 + DCs), 2 major clusters of macrophages (FABP4− macrophages and FABP4+ macrophages) and particularly 3 subtypes of monocytes (CD14+ classical monocytes, CD14+ intermediate monocytes, and CD16 + non-classical monocytes) (Additional file 2: Fig. S3B-3D and Additional file 5: Dataset 3). We also noted that the non-classical monocytes were predominantly from young controls (Additional file 2: Fig. S3E); hence, there was a lack of statistical power to detect expression changes in this subtype (Additional file 2: Fig. S3F), so we restricted our analyses to two CD14+ subtypes of monocytes.
Next, we performed GSEA analysis to uncover perturbed functions associated with COPD, aging, and smoking (see Additional file 1: Supplementary Methods). Intriguingly, monocytes appeared to play very distinct immunological roles under different conditions (Additional file 2: Fig. S3G). In COPD, inflammatory responses, mediated primarily by NF-κB as well as tumor necrosis factor (TNF) signaling, were aberrantly activated in monocytes (Fig. 3A), as well as some downstream processes, such as negative regulation of apoptosis and leukocyte migration (Fig. 3B), a clear indication that monocytes function as a component of innate immunity. (Additional file 2: Fig. S3H), suggesting that at least in a subset of COPD patients, it was monocytes that triggered innate immune-driven inflammation, rather than macrophages, a long-recognized orchestrator in COPD pathogenesis. With regard to aging, inflammation appeared to be low grade, an expected phenomenon termed “inflammaging” [11], as many upstream regulators of the NF-κB and TNF pathways (such as REL and MAP2K3) were marginally increased in aged monocytes, and many pro-inflammatory effectors (such as IL1B, CCL20, and ICAM1) even showed decreased expression (Additional file 2: Fig. S3I and Additional file 6: Dataset 4). Notably, genes involved in interferon gamma (IFN-γ) production, including the typical IFN-γ stimulator IL18, were significantly up-regulated in aged monocytes (Fig. 3C, D), suggesting that interleukin-18 (IL-18)-mediated signaling is responsible for the maintenance of inflammaging. Recent studies have shown that telomere dysfunction can induce the activation of IL-18-mediated signaling via the YAP1 complex [27]. Consistently, our study showed that down-regulated DEGs in aged monocytes were enriched for DNA damage response and telomere maintenance (Fig. 3C), and components of the YAP1 complex were significantly activated, along with 2 deubiquitinases vital to both inflammasome activation and YAP1 stabilization, namely, USP7 and USP47 (Fig. 3D). Regretfully, the expression of YAP1 was too low to detect (Additional file 2: Fig. S3J). Altogether, these findings provide a plausible mechanism underlying the initiation of inflammaging. Other aging-compromised functions might also contribute to inflammaging, such as translation, vesicle transport and secretion (Additional file 2: Fig. S3K), as their malfunctions could lead to retarded immune responses. Collectively, given the well-established role of IL-18 in the development of optimal T helper type 1 (Th1) cell responses via IFN-γ [28], our study revealed that monocytes link innate and adaptive immunity in aged lungs. In striking contrast to COPD and aged monocytes, monocytes in lungs from active smokers mainly functioned as antigen-presenting cells (APCs) in that genes involved in MHC class II antigen presentation were overrepresented in up-regulated DEGs when active smokers were compared with never-smokers (Fig. 3E and Additional file 2: Fig. S3L). Note that our study and others consistently defined CD14 + intermediate monocytes as APCs in the population of monocytes [29] (Additional file 2: Fig. S3C). Concordantly, we observed a large increase in the proportion of the intermediate subtype in active smokers (Fig. 3F). To rule out the possibility that the deviation in cell composition was simply a result of random sampling, we performed trajectory analysis to identify potential driver genes responsible for monocyte differentiation from the classical to the intermediate subtype (Additional file 2: Fig. S3M, S3N, Additional file 7: Dataset 5, and see Additional file 1: Supplementary Methods) and found that smoking-associated DEGs were indeed enriched among driver genes, whereas COPD- or aging-associated DEGs were not (Fig. 3G and Additional file 2: Fig. S3O). Accompanying antigen processing and presentation via MHC class II was T-cell activation, leukocyte chemotaxis, and particularly the dysregulation of the biosynthesis of sphingolipids (Additional file 2: Fig. S3G), which have gained growing attention for their increasingly recognized immunoregulatory role in inflammation, antigen presentation and thereby the pathophysiology of COPD [30,31,32]. We summarized the key points regarding the distinct immunological roles of monocytes in a schematic illustration in Fig. 3H. We also validated these findings by leveraging recently published COPD scRNA-seq datasets [17, 22] and showed that the combined analysis findings were in good agreement with our findings (Additional file 2: Fig. S3P, S3Q and see Additional file 1: Supplementary Methods).
With the rationale that the prioritized cell types likely responded or contributed most to the disease phenotype and therefore their status could to some extent be a proxy to uncover system-level malfunctions, we examined whether other cell types underwent similar or biologically relevant changes conditionally by using the GSEA enrichment score as a metric to measure the degree of the activation or repression of any biological process of interest. In line with our expectations, inflammatory responses and cytokine production were more prevalent in the COPD-associated environment, and antigen presentation via MHC class II was more active in smoking-associated APCs (Fig. 3I). In particular, we observed widespread activation of antigen presentation via MHC class I in aged structural cells, corresponding to the induction of Th1 cell responses mediated by aged monocytes (Fig. 3I). In summary, our analyses demonstrated that the top-prioritized cell type, herein monocytes, was indeed a crucial orchestrator in COPD pathogenesis, especially when aging and smoking were heavily involved.
“COPD core” pathologies: altered intercellular communications and the pro-apoptotic effect of monocytes on AT2s
The activation of pro-inflammatory signaling pathways and cytokine production in COPD-associated monocytes suggested that this prioritized cell type likely drives transcriptional alterations in other cell types via intercellular communication. We therefore performed connectome analyses to identify differences in CCIs between control and COPD lungs. More specifically, we quantified changes in CCIs by measuring differences in the numbers of CellPhoneDB-predicted ligand-receptor pairs among the 15 cell types (see Additional file 1: Supplementary Methods). Overall, CCIs between innate immune cells and non-immune structural cells were increased, which was particularly obvious for monocytes (Fig. 4A and Additional file 2: Fig. S4A), strongly supporting their vital role in modulating microenvironmental signals in COPD lungs. We then focused on CCIs with a > 10% increase in predicted ligand-receptor interactions from control to COPD lungs (Fig. 4B). Interestingly, both types of alveolar epithelial cells (AT1s and AT2s) exhibited the largest enhanced interactions with monocytes (Fig. 4B). A Circos illustration of ligand-receptor pairs responsible for elevated CCIs between monocytes and alveolar epithelial cells as a result of significant expression changes highlighted several pro-inflammatory cytokines up-regulated in COPD-associated monocytes, including IL-1β, TNF, and EREG (Fig. 4C and see Additional file 1: Supplementary Methods). However, the Circos illustration of ligand-receptor pairs accounting for enhanced CCIs between monocytes and other immune cells revealed that the most extensively involved ligands and receptors were all from COPD-associated monocytes, including chemokines (CCL20 and CCL4), adhesion molecules (ICAM1), and chemokine receptors (CCRL2 and CXCR4) (Additional file 2: Fig. S4B). These findings implied that in addition to being reactive to surrounding cells, monocytes within COPD lung tissues dominated their enhanced interactions with both structural cells and other immune cells. Such frequent intercellular communications relying on pro-inflammatory ligand-receptor pairs fit well with the role of monocytes in orchestrating the inflammation core to COPD pathogenesis. To further explore the downstream effects of the interactions between monocytes and alveolar epithelial cells, we performed NicheNet analysis to prioritize ligands that could best predict the DEGs in COPD-associated alveolar epithelial cells (see Additional file 1: Supplementary Methods). Among the top-prioritized ligands, IL-1β and TNF were predicted to be influential upstream ligands that could regulate many genes involved in the NF-κB and mitogen-activated protein kinases (MAPK) pathways in both alveolar epithelial cell types (Fig. 4D and Additional file 2: Fig. S4C). In accordance with this finding, the upstream analysis by Metascape predicted that several transcription factors in NF-κB, IL-1, TNF, and MAPK signaling pathways could potentially have regulated the identified DEGs in AT1s and AT2s (Fig. 4E and see Additional file 1: Supplementary Methods). Collectively, our analyses suggested that elevated CCIs due to activation of pro-inflammatory signaling pathways in COPD-associated monocytes, especially IL-1 and TNF signaling, might be major drivers of alveolar epithelial inflammation in COPD.
Next, to verify the predicted pro-inflammatory effect of monocytes and further clarify the biological consequences of such interactions on alveolar epithelial cells, we performed an in-depth analysis of the transcriptomic changes in alveolar epithelial cells, with a specific focus on AT2s, which ranked second only to monocytes in their contribution to COPD heritability (Fig. 2D) and showed more than 10% enhanced interaction with monocytes uniquely among all cell types in COPD lung tissues (Fig. 4B). To better characterize signaling topologies, we performed WGCNA analysis to find COPD-associated gene clusters with highly correlated expression levels and found a module, referred to as the Blue module, of 838 genes whose overall expression was significantly up-regulated in COPD (Fig. 4F, G, Additional file 8: Dataset 6, and see Additional file 1: Supplementary Methods). We used an enrichment map to reveal functional topologies among major biological processes of the Blue module (Fig. 4H), including external stimuli (cellular response to chemical stimulus), cellular damage (apoptotic process and programmed cell death), and signaling cascades in accordance with those predicted in Fig. 4E (MAPK family signaling cascades), strongly attributing the cellular damage in AT2s to the external immune stimulus from monocytes mediated by pro-inflammatory ligands such as IL-1β and TNF. Subsequently, we constructed a functional interaction network by STRINGDB for the proteins encoded by Blue module genes that were enriched in more than 2 biological processes in Fig. 4H to better demonstrate the key biological processes in COPD-associated AT2s as well as their interrelationships (Fig. 4I and see Additional file 1: Supplementary Methods). Notably, the STRINGDB network not only emphasized the biological processes (such as cellular response to external stimulus, response to unfolded protein, mitochondrion organization, and apoptotic processes) most associated with COPD in AT2s but also contained the predicted signaling pathways (including NF-κB, IL-1, TNF, and MAPK signaling pathways) in AT2s influenced by ligands from monocytes (Fig. 4E). More importantly, this network clearly showed that external stimuli from monocytes, especially represented by IL-1β and TNF, led to the apoptosis of AT2s via endoplasmic reticulum (ER) stress, mitochondrial dysfunction, and the MAPK signaling pathway, wherein several hub proteins were linked by NF-κB signaling (Fig. 4I).
When exploring the possibility that cellular damage in AT2s might compromise their capacities to differentiate into AT1s [33] by slingshot (Additional file 2: Fig. S4D), inspired by the observation that AT2s from COPD lung tissues were more clustered in the early stage of the differentiation path (Additional file 2: Fig. S4E), we further identified potential driver genes responsible for AT2 differentiation and indeed observed an overrepresentation of COPD-associated DEGs in monocytes among them, as well as aging-associated DEGs (Additional file 2: Fig. S4F and see Additional file 1: Supplementary Methods). Notably, a number of pro-inflammatory and pro-apoptotic genes were identified as driver genes, along with a recently discovered AT2 surface marker, CD74, which is helpful for discriminating AT2s from AT1s [33] (Additional file 2: Fig. S4G and Additional file 9: Dataset 7), implying that differentiation suppression of COPD-associated AT2 might be partly accounted for by inflammation and apoptosis.
In summary, our analyses established a signaling pathway underlying the “COPD core” pathogenesis, which started from pro-inflammatory factors secreted by monocytes, transmitted signals by CCIs, and ultimately resulted in the apoptosis of AT2s as well as restrained differentiation towards AT1s and hence the potential failure of alveolar epithelial repair after injury.
Potential mechanism of aging facilitating COPD development: an autoimmune airway niche modulated by aged club cells underlying the injury of airway epithelia
Since club cells were the most prioritized aging-associated cell type (Fig. 2E and F), we focused our attention on club cells to illustrate the potential mechanism by which aging might facilitate COPD development. Prevalent overexpression of MHC class I molecules hinted at increasing susceptibility to autoimmunity during aging (Fig. 3I). Hence, to better characterize the cellular status of each single cell, we performed GSVA analysis to assess single-cell functional enrichment (see Additional file 1: Supplementary Methods). The heatmap listing the top altered biological terms between club cells from lungs of the elderly and young individuals revealed a distinct subset of aged club cells characterized by hallmarks of aging, such as a loss of proteostasis (regulation of IRE1 mediated unfolded protein response) and cellular senescence (negative regulation of mitotic cell cycle), which were mainly derived from COPD lungs (Fig. 5A). Similarly, other airway epithelial cells, including AT1s, AT2s, and ciliated cells, showed signs of these hallmarks with age (Additional file 2: Fig. S5A–S5C). Note that single-cell functional enrichment of the MHC class I antigen presentation pathway was highly correlated with those of negative regulation of the cell cycle and IFN-γ-mediated signaling pathway among all 4 types of airway epithelial cells (Fig. 5B and Fig. S5D–S5F), strongly suggesting a potential association between autoimmune susceptibility as well as IFN reactivity and cellular senescence of epithelial cells and COPD pathology.
Genes encoding IFN receptors (including IFNGR1, IFNGR2, IFNAR1, and IFNAR2) were significantly increased in airway epithelial cells with age (Fig. 5C and Additional file 2: Fig. S5G–S5I), and the up-regulation of these IFN receptors together with their ligands was also observed in the aged nervous system [34,35,36], which reportedly resulted in an increased level of MHC class I molecules in aged neurons [34]. Correspondingly, vital components of the MHC I system, especially β2-microglobulin (B2M), showed a general up-regulation in aged airway epithelial cells (Fig. 5C and Additional file 2: Fig. S5J-5L), which could injure lung epithelia by MHC class I-mediated cytotoxicity and directly induce epithelial cell senescence, as previously reported [37]. Particularly, in club cells, we also observed negative correlations between MHC class I antigen presentation and cell fate commitment as well as the fibroblast growth factor (FGF) receptor signaling pathway, which is indispensable to club cell proliferation and differentiation [38] (Fig. 5B). The expression of FGFR2 indeed declined in aged club cells (Fig. 5C). Given that club cells are multifunctional bronchoalveolar epithelial progenitor cells that are of great importance to airway epithelium maintenance [39], such negative correlations raised the possibility that the autoimmunity of aged club cells may compromise their stemness through the modulation of FGF signaling, as was reported in other types of progenitor cells and stem cells [35]. Concordantly, we further identified an autoimmune-prone sub-cluster unique to aged club cells (Fig. 5D, Additional file 2: Fig. S5M, S5N, and Additional file 10: Dataset 8) and enriched in COPD lungs (Fig. 5A), for which antigen presentation via MHC I and IFN-γ-mediated signaling pathway were remarkably activated (Fig. 5E), whereas cell proliferation and the FGF receptor signaling pathway were suppressed (Additional file 2: Fig. S5O and S5P) compared to other aged club cells. Next, after dissecting the inferred differentiation trajectories from club cells to other epithelial cells (Fig. 5F, Additional file 2: Fig. S5Q, Additional file 11: Dataset 9, and see Additional file 1: Supplementary Methods), we found that aged club cells were more clustered in the earlier phase in all three differentiation modes and even more so for the autoimmune-prone sub-cluster of aged club cells (Fig. 5G). In addition, we observed an overrepresentation of genes regulating proliferation among differentiation driver genes, as well as aging-associated DEGs related to proliferation, such as down-regulated pro-proliferative DEGs (including SOX2, TMEM45A, CRNDE, and CD55) and up-regulated anti-proliferative DEGs (including IGFBP7 and TIMP2) in aged club cells, the degree of which was even greater in the autoimmune-prone sub-cluster (Additional file 2: Fig. S5R). Collectively, we revealed an autoimmunity-stimulating IFN/MHC I axis generally activated in aged lung epithelial cells, which both directly led to epithelial cell injury and specifically dampened the stemness of club cells by inhibiting FGF signaling.
Since our data demonstrated cellular senescence and an autoimmune response in non-immune structural cells from aged lungs, we further investigated whether and how other cell types were involved in these processes by analysing the CCIs. Intriguingly, aged club cells displayed the strongest enhanced communications with other cell types (Fig. 5H), suggesting a non-negligible role of this aging-associated prioritized cell type in modulating the microenvironment that increases susceptibility to COPD. Hence, we took a deeper look into the ligand-receptor pairs responsible for enhanced CCIs between club cells and other cell types (Fig. 5I). In general, cytotoxic responses to senescent cells were augmented in aged immune cells, particularly monocytes, DCs, NKCs and T cells, supported by significantly increased expression of ligands that promote CD8+ T-cell and NKC cytotoxic responses, such as CRTAM, HAVCR2, TNFRSF1A, BST2, and FAS. It is noteworthy that aged monocytes displayed the largest increase in CCIs with aged club cells (Fig. 5H), suggesting that the possible monocyte recruitment of aged club cells may further amplify the epithelial damage attributed to autoimmune responses and that monocytes may simultaneously elevate the autoimmunity of aged club cells in other direct ways in addition to IL-18. Aged club cells also presented enhanced NOTCH2-mediated interactions with other epithelial cells (Fig. 5I), implying that other types of aged epithelial cells might prevent the replenishment of themselves from club cells by changing the differentiation direction of club cells into goblet cells [40], a cell type not detected in our data. With regard to CCIs between club cells and endothelial cells, we proposed that the increased secretion of CD47 [41] and carcinoembryonic antigen-related cell adhesion molecule 1 (CEACAM1) [42] from aged club cells was likely to promote vascular damage and aging processes. In summary, these findings illustrated an autoimmune airway niche during aging responsible for the cytotoxic damage of airway epithelia, which was triggered by the autoimmunity-inducing IFN/MHC I axis, further augmented by monocytes, and modulated by aged club cells.
Potential mechanism of smoking facilitating COPD development: alveolar impairment due to detrimental intercellular communication between macrophages and endothelial cells
GSEA highlighted the activation of antigen presentation via MHC class II among innate immune cells (mainly monocytes, macrophages, DCs, and neutrophils) in the lungs from active smokers (hereafter referred to as “smoking lungs”) (Fig. 3I). A deeper examination of the functional enrichment of these cells under smoking by GSVA confirmed this observation, as many processes indispensable to MHC class II trafficking and processing were also enriched, such as endocytosis and vesicle transportation (Additional file 2: Fig. S6A–S6D). We also noted that macrophages, the top-prioritized cell type associated with smoking, appeared to be much more heterogenous than other cell types (Additional file 2: Fig. S6A), partly because this cell type was most abundantly sampled in our dataset. Hence, differential expression-based analysis was not adequate to explore alterations in macrophages. We therefore addressed heterogeneity by k-means clustering based on the top 400 highly variable GSVA terms and segregated macrophages into 3 subtypes (A, B, and C) (Fig. 6A and see Additional file 1: Supplementary Methods). The GSVA enrichment of subtypes A and B fit well with a role in the innate immune response and antigen presentation via MHC class II, respectively, thereby largely corresponding to the homeostatic function of macrophages (Fig. 6B). Subtype C, however, was distinguished by the expression of genes involved in mitochondrial electron transport (Fig. 6C), suggesting a distinctive cellular metabolism of this subtype. Intriguingly, subtype C comprised a disproportionate number of cells from active smokers (odds ratios (ORs) > 3.5, p-value < 1.1e−172, one-sided Fisher’s exact test, alternative = “greater”; Fig. 6A). Moreover, subtype C was nearly absent in never-smoking COPD patients (ORs < 0.51, p-value < 7.8e−09, one-sided Fisher’s exact test, alternative = “less”; Fig. 6A), strongly suggesting that subtype C might be an important smoking-induced subtype of macrophages that contributes to COPD pathogenesis.
Along with mitochondrial gene expression, we also observed an overexpression of several genes encoding enzymes implicated in sphingolipid metabolism in subtype C macrophages (Additional file 2: Fig. S6E). Inspired by the pleiotropic effect of bioactive sphingolipids on multiple cell types in the lung [43], we sought to determine whether and how sphingolipid metabolism in macrophages plays a role in smoking lungs. In addition to macrophages, altered sphingolipid metabolism may affect endothelial cells, which was implicated by the enrichment of smoking-associated DEGs from endothelial cells in “cellular response to lipid” (Fig. 2B). Indeed, we found evidence for enhanced pro-inflammatory sphingolipid biosynthesis (including ceramide and sphingosine) along with the diminished biosynthesis of anti-inflammatory sphingolipids (including sphingosine-1-phosphate and ceramide-1-phosphate) in macrophages from smoking lungs, as demonstrated by increased macrophage expression of SPTLC1, SPTLC2, CERS2, CERS5, ASAH1, and SGPP1 as well as decreased expression of CERK and SMPD4 in smoking lungs (Fig. 6D). Additionally, such expression changes in these genes were observed in cigarette smoke extract (CSE)-treated bone-marrow-derived macrophages (BMDMs) from C57BL/6 mice (Fig. 6E and see Additional file 1: Supplementary Methods), further supporting pro-inflammatory sphingolipid biosynthesis in macrophages from smoking lungs.
Next, given the vesicle-transportation-dependent paracrine activities of ceramide and sphingosine reported in other tissues [44] and their detrimental role in lungs exposed to cigarette smoke [43], we hypothesized that excessive ceramide and sphingosine biosynthesized by macrophages could also have an effect on endothelial cells in smoking lungs. However, such intercellular communications based upon sphingolipid metabolites could not be assessed by the analytic strategy relying on known ligand-receptor pairs. Hence, we designed a novel approach to infer the potential contribution of activated sphingolipid biosynthesis in smoking macrophages to transcriptomic alterations in other cell types, which takes ligand-receptor pair-mediated CCIs into account at the same time due to the participation of the other 2 subtypes of macrophages in innate and adaptive immune responses (Fig. 6B). Briefly, we first summarized cell-type-specific gene expression at the individual level, and then we examined whether the expression of significantly dysregulated sphingolipid metabolic enzymes as well as up-regulated ligands in smoking macrophages changed concordantly with smoking-associated DEGs in other cell types across individuals (Fig. 6F and see Additional file 1: Supplementary Methods). We found the highest concordance of expression change between smoking-associated DEGs in endothelial cells and the abovementioned dysregulated sphingolipid metabolic enzymes and ligands in macrophages (Fig. 6G). Note that the sphingolipid biosynthesis process was not activated in endothelial cells, which ruled out the autocrine effect of sphingolipids secreted by endothelial cells (Additional file 2: Fig. S6F–S6H). To gain insight into the mechanism underpinning such metabolic communication, we sought to identify correlated genes across individuals with the 8 dysregulated sphingolipid metabolic enzymes in macrophages not only in endothelial cells to infer cell-nonautonomous effects but also in macrophages to reflect cooperative cell-autonomous alterations (Additional file 12: Dataset 10). As seen from the Sankey plot in Fig. 6H, the enriched biological processes among the correlated genes to these enzymes in macrophages included autophagy, vesicle transport, and endomembrane system organization, while those in endothelial cells included membrane trafficking, cellular response to lipid, and cytoskeleton organization, indicating that pro-inflammatory sphingolipid-based metabolic communication is likely to be transmitted in the form of vesicles from macrophages to endothelial cells and may have a damaging effect on both cell types in smoking lungs. Several enzymes and their highly correlated genes from both macrophages and endothelial cells are shown in Fig. 6I. On the other hand, we applied the same strategy to prioritize up-regulated ligands in smoking macrophages that could best account for expression changes in endothelial cells (see Additional file 1: Supplementary Methods). Specifically, we first compiled lists of highly correlated genes in endothelial cells with each smoking-associated up-regulated ligand in macrophages across individuals, and then we examined the functional overlap of each list with the smoking-associated DEGs of endothelial cells (Additional file 2: Fig. S6I). Two ligands in macrophages, namely, PDGFC and TWSG1, showed remarkable functional overlap with the smoking-associated DEGs of endothelial cells (Additional file 2: Fig. S6I), underscoring the importance of these two ligands in mediating the effect of macrophages on endothelial cells under smoking. Similarly, PDGFC, TWSG1 and their representative correlated genes from endothelial cells are shown in Additional file 2: Fig. S6J.
Overall, these findings represented a scenario in which alveolar impairment by smoking was mediated by APCs due to their enhanced capacity for antigen presentation and heightened intercellular communication between macrophages and endothelial cells in both ligand- and sphingolipid-metabolite-dependent manners, which may exacerbate COPD.
Validation of findings under the three pathological conditions by a meta-analysis
Given the high clinical heterogeneity of COPD and the single-centre source of the lung tissue samples in our study, we sought to confirm our findings, especially those of system-level malfunctions centred on prioritized cell types under three conditions (COPD, aging, and smoking), in publicly available COPD datasets (see Additional file 1: Supplementary Methods). To this end, we first combined the scRNA-seq dataset of 9 samples from both healthy and COPD individuals with different ages, races, and smoking histories (Additional file 13: Dataset 11) recently published by Adams et al. [22] with our scRNA-seq dataset, forming a mixed dataset (hereafter referred to as the “9 + 9 dataset”). We then repeated the downsampling analysis on the 9 + 9 dataset. Once again, monocytes and club cells were prioritized in this combined dataset under COPD and aging conditions, respectively (Fig. 7A). However, the result under smoking conditions in the 9 + 9 dataset differed slightly from that in our dataset (Fig. 7A), which was largely due to the heterogeneity of macrophages in both datasets. This result demonstrates the rationality and necessity of our analysis strategy being based on the GSVA enrichment score for subtype classification of macrophages, which can be extended to investigate the influence of smoking on macrophages in the 9 + 9 dataset in the near future. Similarly, we also measured the activation or repression states of the same biological processes in Fig. 3I across all the cell types under the three conditions in the 9 + 9 dataset by using GSEA enrichment scores (Fig. 7B). Although the 9 + 9 dataset did not fully recapitulate the status of every biological process in our dataset, these datasets showed substantial overlap in heightened inflammatory responses together with the cytokine production prevalent in monocytes and AT2s in COPD lung tissues and shared the general activation of antigen presentation via MHC class I along with DNA and telomere dysregulation among structural cells in aged lung tissues. The states of biological processes were more heterogeneous under the condition of smoking than under aging and COPD due to the high complexity in the impact of smoking, although sphingolipid biogenesis in smoking-associated macrophages as well as MHC class II-mediated antigen presentation among smoking-associated APCs were still highlighted in the 9 + 9 dataset.
Overall, this integrated analysis of published scRNA-seq data indicated that the conditional prioritized cell types and system-level malfunctions reflected by the dysregulated biological processes we found were present in people with different genetic backgrounds and clinical contexts, confirming the reliability and reproducibility of our results.