Alveolar macrophage transcriptomic profiling in COPD shows major lipid metabolism changes.

BACKGROUND: Immune cells play a major role in the pathogenesis of COPD. Changes in the distribution and cellular functions of major immune cells, such as alveolar macrophages (AMs) and neutrophils are well known; however, their transcriptional reprogramming and contribution to the pathophysiology of COPD are still not fully understood.
METHOD: To determine changes in transcriptional reprogramming and lipid metabolism in the major immune cell type within bronchoalveolar lavage fluid, we analysed whole transcriptomes and lipidomes of sorted CD45+Lin-HLA-DR+CD66b-Autofluorescencehi AMs from controls and COPD patients.
RESULTS: We observed global transcriptional reprogramming featuring a spectrum of activation states, including pro- and anti-inflammatory signatures. We further detected significant changes between COPD patients and controls in genes involved in lipid metabolism, such as fatty acid biosynthesis in GOLD2 patients. Based on these findings, assessment of a total of 202 lipid species in sorted AMs revealed changes of cholesteryl esters, monoacylglycerols and phospholipids in a disease grade-dependent manner.
CONCLUSIONS: Transcriptome and lipidome profiling of COPD AMs revealed GOLD grade-dependent changes, such as in cholesterol metabolism and interferon-α and γ responses.

Introduction

COPD is characterised by progressive airflow obstruction, inflammation in the airways and systemic comorbidities, such as cardiovascular diseases and diabetes [1]. It is a leading cause of morbidity and mortality worldwide [2], and induces a substantial and increasing economic and social burden [3, 4]. Current guidelines for COPD patient classification use the Global Initiative for Chronic Obstructive Lung Disease (GOLD) criteria on the basis of spirometry-estimated severity of airflow limitation and range from mild (GOLD grade 1) to very severe (GOLD grade 4) disease [5]. Cigarette smoking is the leading environmental risk factor for COPD, yet even for heavy smokers, fewer than 50% develop COPD during their lifetime, indicating that the disease initiates from a complex interaction between environment and genome. Some genetic factors, e.g. deficiency of α-1 antitrypsin [6] or single-nucleotide polymorphisms in MMP12 [7] have been associated with the risk to develop COPD.

Innate immune cells are key players in protecting the lung from airway infections and their impairment plays a major role in COPD pathogenesis and exacerbation episodes [1, 8, 9]. Although changes in cellular distribution and functions within the alveolar space have been described in COPD before [10], these findings are far from being complete. Furthermore, the advent of omics technologies even extending to metabolomics [11], as well as a better understanding of cell type classification requires revisiting these aspects of innate immune cell biology of the alveolar space, both in healthy individuals and COPD patients. Alveolar macrophages (AMs) are described to be the most abundant cells in humans [12]. Under physiological conditions, AMs are the major cell type on the luminal surface of the alveolar space and are major regulators of the initiation and resolution of inflammation [13, 14] and surfactant metabolism [15]. AMs catabolise pulmonary surfactants, which are composed of phospholipids and proteins and humans lacking functional AMs develop pulmonary proteinosis [16].

Recent reports have shown that in respiratory diseases, such as asthma, COPD and pneumonia a lipidomic remodelling of the bronchoalveolar lavage fluid (BALF) is observed, with changes in the levels of cholesterol, sphingomyelins, phospholipids and fatty acids being the most prominent [17-19]. In particular, the role of sphingolipids and cholesterol metabolism are believed to play a role in the pathogenesis of COPD [20, 21]. Nevertheless, changes of the lipid metabolism within AMs in COPD patients have not been addressed so far.

In this study, we hypothesised that transcriptional changes leading to altered lipid metabolism occur in a disease grade-dependent manner. Using multi-colour flow cytometry (MCFC), we defined AMs as the major cell type in COPD patients for which we provide clear evidence of global transcriptional reprogramming, which differed between GOLD grade 2 and GOLD grade 3/4 COPD patients. Among the major alterations in the AM transcriptome of COPD patients, we identified genes associated with lipid metabolism, which was linked to significant changes in several lipid species by lipidomics analysis.

Methods

Human specimens

Human studies were approved by the ethics committees of the University of Bonn and University Hospital Bonn (local ethics vote 076/16). All patients provided written informed consent according to the Declaration of Helsinki before specimens were collected. Patients with COPD were diagnosed and stratified according to the guidelines of GOLD [5]. Eligible patients were aged 18 years or older and were either current, past or nonsmokers (table E1). Age-matched individuals suffering from chronic idiopathic cough, demonstrating an exquisitely sensitive cough reflex without underlying pathology [22], served as control donors. Patients with other pulmonary diseases (termed “other” (figure 1a)) were diagnosed as asthma, asthma–COPD overlap, bronchiectasis, cancer, fibrosis, pneumonia and sarcoidosis (table E2), but were excluded from further evaluation within this study.

FIGURE 1

Identification of myeloid immune populations in the bronchoalveolar lavage fluid (BALF) of COPD and control patients. Schema of the pipeline for collection and processing of human BALF samples. Of 177 BALF samples, 69 satisfied the quality criteria (>30% recovery rate and no blood and/or mucus contamination) and were included in the study for multiple colour flow cytometry, transcriptomics and/or lipidomics. a) 61 patients suffering from asthma, asthma–COPD overlap (ACO), bronchiectasis, cancer, fibrosis, pneumonia or sarcoidosis were excluded from the study. b) Representative multicolour flow cytometry (MCFC) analysis of the myeloid compartment of a control patient. c) Absolute numbers of myeloid immune cells in the BALF of control, Global Initiative for Chronic Obstructive Lung Disease (GOLD) stage 2 and GOLD3/4 COPD patients calculated with traditional gating. Data are from 8–29 patients per group except for mast cells GOLD3/4 (n=2) and are represented as mean±sd. Data are represented as mean±sem and were analysed with an unpaired two-tailed t-test. d) Giemsa staining of cytospins for sorted BALF myeloid cell types. RR: recovery rate; AM: alveolar macrophage; DC: dendritic cell; . Scale bars=10 µm. *: p<0.05.

Bronchoscopy procedure

Bronchoscopy was performed as a part of the diagnostic workup by two bronchoscopists through oral access and with light conscious sedation and was performed in the middle lobe or, if not accessible, the lingular lobe. Warmed saline (six syringes of 20 mL each) was instilled into the airways to enable BALF recovery.

BALF processing

Human BALF was obtained from all patients included in the study (control, COPD, other) through bronchoscopy. BALF specimens were washed with PBS, suspended with 0.02% EDTA-2Na and washed again for final re-suspension with 2% fetal calf serum (FCS)/1 mM EDTA. CD45+Lin−HLA-DR+CD66−Autofluorescence+ AMs were sorted using a FACS Aria III cell sorter (BD Biosciences, USA).

Cell counting

Total cell counts were determined with (1:5) Trypan Blue exclusion (Sigma-Aldrich) under an optical microscope. BALF cells were diluted 1:10 in the Trypan Blue solution and counted in a Neubauer haemocytometer.

Flow cytometry/FACS

Single-cell suspensions were stained with Live/Dead yellow fluorescent dye (Thermo Fisher Scientific, USA) for 15 min at room temperature and were washed with PBS at 300×g for 5 min at 4°C. They were then re-suspended in 100 µL PBS and blocked with 5 µL human FcR blocking reagent (Miltenyi, Germany) for 15 min on ice and were subsequently stained with the listed anti-human antibodies (table E4) in buffer containing PBS, 2% FCS, 1 mM EDTA for 30 min on ice. Cells were spun at 300×g for 5 min at 4 °C and re-suspended in buffer containing PBS, 2% FCS, 1 mM EDTA for analysis. Data acquisition was performed on a FACS Aria III cell sorter. Analysis was performed with FlowJo v.10 software (Tree Star, USA).

Cytospin preparation

Cytospins were obtained by centrifuging 2×105 cells in 200 µL PBS on microscope slides at 20% power for 5 min. Excess buffer was carefully discarded and slides were air dried for 3 min followed by fixation with 100% methanol for 5 min at 4°C. The slides were subsequently washed with PBS and stained with 1:20 Giemsa solution (Sigma, USA) for 25 min at room temperature. A final rinsing step with H2O and air drying before mounting followed. Cell morphology was examined by microscopic evaluation of stained cells using an Axio Lab A1 microscope (Zeiss, Germany).

RNA extraction and library preparation

Total RNA was isolated from human AMs with the miRNeasy Micro kit (QIAGEN, Germany) according to the manufacturer's protocol. cDNA libraries were prepared from 5 ng total RNA with the SMART-seq2 protocol [23] and were tagmented with the Nextera XT kit (Illumina, USA). Library size selection was carried out with AMPure beads (Beckman-Coulter, USA). Libraries were sequenced for SR 75 cycles on a NextSeq500 system (Illumina) using High Output v2 chemistry. Base call files were converted to fastq format and demultiplexed using bcl2fastq v.2.20.

Data pre-processing and RNA-sequencing analysis

The 75 bp single-end reads were aligned to the human reference transcriptome hg38 from UCSC by kallisto v.0.44.0 using default parameters. Data were imported into DESeq2 (v.1.10.1; [24]) using the TXimport (v.1.2.0, [25]) package. DESeq2 was used for the calculation of normalised counts for each transcript using default parameters. All normalised transcripts with a maximum over all group means <10 were excluded resulting in 33 032 present transcripts. Unwanted or hidden sources of variation, such as batch, sex and smoking status were removed using the sva package [26]. The normalised rlog-transformed expression values were adjusted according to the five surrogate variables identified by sva using the function removeBatchEffect from the limma package [27]. Differentially expressed (DE) genes were defined by a p-value cut-off of 0.05 and an adjusted p-value (independent hypothesis weighting) <0.5.

Gene set enrichment analysis

Gene set enrichment analysis (GSEA) [28] was performed on all present genes of the dataset using the gene ontology set of biological processes. Information of gene ontology was obtained from the biological process gene set “c5.bp.v7.0.symbols.gmt”, downloaded from the Molecular Signatures Database (MSigDB). All present genes were used as background (universe).

Linear support vector regression

Linear support vector regression [29] was employed to characterise the relative contribution of 28 different activation signatures derived from Xue et al. [30] to the control and COPD patients (1000 permutations).

Filtering for transcription factors, epigenome, surfaceome and secretome

All present transcripts were filtered and sorted by their variance in the dataset. The 20 most variable genes of each category were selected and visualised in heat maps. Transcription factor (TF) lists were extracted from Fulton et al. [31], the epigenome gene list was derived from the literature, surface and secretome markers were extracted from the Human Protein Atlas [32, 33].

Construction of coexpressed network analysis: automated

Coexpressed network analysis (CoCena2) was performed to elucidate similarities and differences within the gene expression patterns of the three different patient groups. Pearson correlation was calculated on the 6000 most variable genes within the dataset using the R package Hmisc (v.4.3-0; [34]). Data were filtered for significant (p-value <0.05, Bonferroni correction p<0.05) and positive (r-value >0) correlation values. The group fold change (GFC) was calculated for each gene and each condition on the inverse logarithmic count data using the R package gtools (v.3.8.2; Unbiased clustering was performed using the R package igraph (v.1.2.4.1; [35]). The clustering algorithm “cluster_louvain” was selected as it achieves the highest modularity score. Clusters with fewer than 35 genes are not shown. Network generation was performed with the R package igraph. The network information was imported to and exported from Cytoscape using the R package RCy3 (v.2.6.2; [36]).

Biological function-related bioinformatic analysis of network modules

GSEA was performed on the patient group-related modules identified by CoCena2 using the R package ClusterProfiler (v.3.12.0; [37]). Information of hallmark genes was obtained from the hallmark gene set “h.all.v6.1.symbols.gmt”, downloaded from MSigDB. All genes present in the network were used as background (universe).

Lipidomics

Sorted AMs (5×104) pelleted and frozen at −80°C until analysis. Extraction mix (chloroform 1:5 methanol-containing internal standards) was spiked to the pellets before sonication. Samples were treated with chloroform and 1% acetic acid, the lower phase was transferred after centrifugation, and let evaporate in the vacuum concentrator (45°C for 10 min). After addition of spray buffer and sonication, samples were analysed separately with a Thermo Q Exactive Plus spectrometer equipped with the HESI II ion source for shotgun lipidomics.

Lipidomics analysis

LipidXplorer software was used for analysis using custom mfql files to identify sample lipids and internal standards. Absolute amounts were calculated using the internal standard intensities followed by normalisation on the sum of all measured lipid species per sample. %mol values were averaged for each patient group, log2-transformed and then used for fold change calculations. CoCena2 was used to find patient group-specific coregulated lipid species using Pearson correlation and clustering using the “cluster_louvain” algorithm.

Statistics

A two-tailed Welch's unpaired t-test was used to analyse data from two groups. Equality of population variance was assessed with the F-test statistic for two independent groups. A nonparametric Wilcoxon test was used to perform a pairwise comparison between patient groups for all enriched macrophage activation signatures. For more than two groups, normality and homoscedasticity were first assessed using the Shapiro–Wilk and Levene tests in R (v.3.6.1). A nonparametric Kruskal–Wallis test with Dunn's multiple correction post hoc was used in figure E3 because the data did not follow a normal distribution. Statistical significance was inferred when p<0.05.

For a full list of the methods, see the supplementary material.

Results

Macrophages are the most abundant cell type in the alveolar space in COPD

We first evaluated the cell type distribution in the BALF to ensure that AMs are the most abundant cell type in control and COPD samples, the latter of which an enrichment of neutrophils had been postulated [38]. We performed MCFC integrating markers to determine all major immune cell types expected within the alveolar space [39]. Individuals with chronic coughing served as controls (table E1). Only BALF samples of the highest quality following quality criteria established previously for the processing of BALF [40], such as recovery rate >30% and absence or minimal blood/mucus contamination were included in the study (figure 1a). Out of 177 screened BALF samples in the clinic, 72 were derived from COPD patients, 44 from control individuals, while patients suffering from asthma, ACO, bronchiectasis, cancer, fibrosis, pneumonia or sarcoidosis were excluded (table E2). 69 patients (36 COPD and 33 controls) finally qualified for further analysis (table E3).

We stained for markers of the myeloid cell compartment (figure 1b and table E5), whereas lymphocytes were labelled with antibodies against CD3, CD19 and CD56 (figure 1b). Our data show that AM numbers undergo a nonsignificant increase in COPD patients, whereas we could detect a statistically significant upregulation (p<0.05) in the neutrophil numbers of GOLD3/4 grade patients compared to controls (figure 1c). Eosinophil and monocyte numbers varied widely within groups and reached no significance, whereas mast cells were found to be significantly higher (p<0.05) in GOLD3/4 compared to GOLD2 COPD patients (figure 1c).

Lymphoid cells were of much lower frequency and using a second panel for lymphoid cell markers (figure S1a and table E5), we observed that only innate lymphoid cells (ILCs) were statistically significantly (p<0.05) changed in all COPD samples (figure S1b and c), with a similar trend when taking GOLD grades into account (figure S1b). Within ILCs, ILC1s were the only significantly increased (p<0.05) ILC population in COPD patients (figure S1d). Cell type classification was also confirmed on sorted cells by Giemsa staining of cytospins (figure 1d). Collectively, AMs are the most abundant cell type in BALF, which was then our focus for subsequent transcriptome and lipidome analyses.

Transcriptional changes in AMs within the alveolar space of COPD patients

To elucidate molecular changes of AMs as the most prevalent cells within the alveolar space, we sorted CD45+Lin−HLA-DR+CD66b−Autofluoreschencehi AMs from the BALF of nine COPD and six control patients and generated whole transcriptomes by RNA sequencing (figure 2a). Principal-component analysis (PCA) of all present genes (figure 2b) and hierarchical clustering of the 25% most variable genes within the dataset (figure S2a) revealed a clear disease grade-dependent structure within the data suggesting significant differences between COPD GOLD2, GOLD3/4 and control samples.

FIGURE 2

Whole transcriptome analysis of alveolar macrophages (AMs) reveals lipidomic dysregulation in COPD patients. a) Schematic representation of the bioinformatics workflow and methods used for the transcriptome analysis of 15 RNA sequencing (RNA-Seq) samples from Global Initiative for Chronic Obstructive Lung Disease (GOLD) stage 2, GOLD3/4 COPD and control patients. b) Principal component (PC) analysis of all 33 032 present transcripts. c) Distribution and amount of differentially expressed (DE) genes. d) Volcano plots of the comparisons GOLD2 versus control, GOLD3/4 versus control and GOLD2 versus GOLD3/4. Upregulated genes are marked in red, downregulated genes in blue. Fold change cut-off was set to 1.5, adjusted p-value cut-off was set to 0.05. TF: transcription factor; Epi.: epigenetic; GSEA: gene set enrichment analysis.

Next, we calculated DE genes between GOLD2 or GOLD3/4 COPD and control patients, respectively. There were 563 DE genes (229 up, 334 down, p<0.05, fold change (FC) >1.5) between GOLD2 and control, 1456 (618 up, 838 down, p<0.05, FC>1.5) between GOLD3/4 and control and 1147 (572 up, 575 down, p<0.05, FC>1.5) between GOLD2 and GOLD3/4 COPD patients (figure 2c), which we visualised in Volcano plots for all pairwise comparisons (figure 2d).

We first assessed overall transcriptional changes. We provide information for the top 20 altered TFs, epigenetic modulators, cell surface molecules (surfaceome) and soluble mediators as micro-environmental interactors (figure S2b–e). Upregulated TFs included ETS2 in COPD patients consistent with its function as a biomarker for loss of lung function [41] (figure S2b). Consistent with the reported dysregulation of histone deacetylase (HDAC) and sirtuin epigenetic signalling in AMs from COPD patients [42-45], the expression of histone-modifying enzymes HDAC9 and SETD1A was downregulated in GOLD3/4 COPD (Figure S2c). Regarding the surfaceome of the AMs, we measured decreased antigen presentation molecules (HLA-DQA2 and CD1B) [46, 47] and CCR2 expression in GOLD3/4 COPD patients (figure S2d). Notably, expression of the ligand chemokine to CCR2, CCL2, was elevated in the secretome-associated gene set of GOLD3/4 COPD patients, similar to other members of the CC family of chemokines, including CCL8 and CCL20 (figure S2d), highlighting the role of chemokines in COPD pathophysiology [48]. Finally, we detected the downregulation of MMP7 and MMP9 in GOLD3/4 COPD patients, with MMP12 being expressed significantly higher only in GOLD2 COPD patients (figure S2e), in agreement with the role of metalloproteinases in the progression of COPD [49-52].

Activation state analysis indicates changes in lipid metabolism

We have previously reported that macrophage activation follows a multi-dimensional model of cell activation [30]. However, for AMs in COPD, it has been proposed that AMs in COPD follow a dual polarisation model with deviation towards an anti-inflammatory phenotype [47, 49, 53, 54] with no hints on alterations in lipid metabolism. To analyse how the observed transcriptional changes in COPD AMs related to our multi-dimensional model of activation states, we performed linear support vector regression utilising our previous human macrophages activation spectrum model based on 29 conditions [30] (figure 3a). Our activation spectrum model could be grouped into nine major activation programs (c1–9) and were used here as the bait to deconvolute the AM transcriptomes from COPD and control patients. The COPD and control transcriptomes were enriched in a spectrum of activation signatures ranging from pro-inflammatory to anti-inflammatory, as well as those extracted from fatty acid-stimulated macrophages (c5 module signature), which show deviations in lipid metabolism. AMs from GOLD2 COPD patients exhibited significantly lower enrichment (p<0.05) of the c4 module signature (P3C/PGE2, P3C and PGE2 stimulation) compared to GOLD3/4 COPD and control patients and c8 module signatures (TNF-α/P3C, TNF-α and P3C stimulation) in comparison with GOLD3/4 COPD patients (figure 3a). In addition, COPD GOLD2 AMs showed a higher enrichment (p<0.05) of the c5 module signatures (lauric acid (LA), linoleic acid (LiA), oleic acid (OA), palmitic acid (PA) and stearic acid (SA) stimulation) compared to AMs from control patients (figure 3a).

FIGURE 3

Overlay of human macrophage activation states in alveolar macrophages (AMs) from COPD patients show lipid term enrichment. Overall, 29 macrophage activation signatures from [30] were used as input to identify the relative fraction of these activation signatures in COPD and control patients. Grouping of activation signatures in nine clusters was used as proposed in [30]. a) Data were statistically analysed with the nonparametric Wilcoxon test. Gene set enrichment analysis (GSEA) using the gene ontology (GO) database filtered for “fat” and “lipid” as reference signatures. b) Normalised enrichment scores (NES) and enrichment p-values for the five most enriched gene ontology (GO) terms for the comparison Global Initiative for Chronic Obstructive Lung Disease (GOLD)2 COPD versus control and c) GOLD3/4 COPD versus control patients. Heat map of the mean expression of the 20 most variable genes filtered by d) the GO term fatty acid metabolism and e) fatty acid catabolic process. IFN: interferon; HDL: high-density lipoprotein; IL: interleukin; GC: glucocorticoid; upLPS: ultra pure LPS; IC: poly(I:C); P3C: Pam3CSK4; PG: prostaglandin; LA: lauric acid; LiA: linoleic acid; OA: oleic acid; PA: palmitic acid; SA: stearic acid; TNF: tumour necrosis factor; sLPS: soluble LPS; TPP: TFNF/PGE2/P3CXXX;. *: p<0.05; **: p<0.01.

Based on these results indicating potential changes in lipid metabolism-associated-genes, we next performed GSEA for either GOLD2 or GOLD3/4 COPD patients compared to controls (figure 3b and c). Lipid-related gene sets, such as fatty acid catabolic process, fatty acid oxidation and regulation of cholesterol biosynthetic process were enriched in AMs from GOLD2 COPD patients (figure 3b). AMs from GOLD3/4 COPD patients presented with enrichment in cholesterol and lipid storage compared to controls (figure 3c).

To highlight the genes that are involved in the lipid-associated phenotype of COPD AMs, we visualised the 20 most variable genes for the gene ontology term “fatty acid metabolism” of the c5 module from figure 3a (figure 3d) and “fatty acid catabolic process” from figure 3b (figure 3e). Several genes that are involved in cholesterol homeostasis were dysregulated in COPD patients. For example, enzymes involved in fatty acid metabolism, such as ACAT2, FAAH, MLCYD and ACOX2 or enzymes that synthesise diacylglycerols (LPIN3) were expressed higher in GOLD2 COPD patients. On the contrary, the acyltransferase CRAT and the triglyceride synthesis enzyme LPIN2 were overexpressed in GOLD3/4 COPD patients. Fatty acid synthase FASN and the anti-inflammatory TF ATF3 [55, 56] were both downregulated in COPD patients irrespective of GOLD grade (figure 3d and e). These findings strongly support that AMs do not follow a simple polarisation model, but rather display a multi-activation phenotype with a shift towards lipid activation signatures.

Coexpression analysis reveals changes in cholesterol homeostasis in COPD

Coexpression network analysis is an alternative to statistical methods describing changes in transcriptome data, such as the calculation of DE genes based on a set fold change and p-value [57]. To further investigate active gene programs in AMs, we performed gene coexpression analysis using CoCena2 (figure 4a). The clustering of the genes within the constructed network led to the identification of groups of transcriptionally similarly regulated genes within the dataset. The mean GFC of these gene modules in GOLD2, GOLD3/4 COPD and control patients were visualised in a heat map (figure 4b). We detected 11 modules, three of which demonstrated GOLD grade-associated average expression; for instance, modules 7 and 8 had a higher GFC in GOLD2 COPD patients, modules 3 and 4 were GOLD3/4-related, while modules 9–11 had a higher GFC in AMs from control patients.

FIGURE 4

Coexpression network analysis uncovers lipid metabolism-associated functions for alveolar macrophages (AMs) in COPD. a) Schematic representation of the bioinformatics workflow for coexpression network analysis of the 6000 most variable genes in the dataset by CoCena2. b) CoCena2 cluster-condition heat map. c) Significant enrichment (q-value <0.1) of hallmark genes within patient groups. Patient group modules consist of the respective cluster-condition heatmap modules. d) CoCena2 network with Global Initiative for Chronic Obstructive Lung Disease (GOLD) stage 2-associated genes within the hallmark term “Cholesterol homeostasis” marked by black edging and coloured according to their cluster name. Each node represents one gene in the network and each edge represents a coexpression. The expression levels are presented in more detail in the following heat map. e) CoCena2 network with GOLD3/4-associated genes within the hallmark term “interferon-α response” marked by black edging and coloured according to their cluster name. Each node represents one gene in the network and each edge represents a coexpression. The expression levels are presented in more detail in the following heat map. f) Genes enriched within the hallmark term “G2M checkpoint” marked by black edging and coloured according to their cluster name. Each node represents one gene in the network and each edge represents a coexpression. The expression levels are presented in more detail in the following heat map. GFC: group fold change; GSEA: gene set enrichment analysis; DN: downregulated; IL: interleukin; JAK: Janus kinase; STAT: signal transducer and activation of transcription.

Next, GSEA was applied to concatenated modules per disease group (control: 1, 2, 6, 9, 10, 11, GOLD2: 5, 7, 8, GOLD3/4: 3, 4) to link the coexpressed genes to biological processes (figure 4c). While genes of control patients were enriched in apical junctions, genes from the modules enriched in AMs derived from GOLD2 patients suggested changes in cholesterol homeostasis in addition to the findings based on statistical analysis (figure 3). However, GSEA of modules enriched in AMs derived from GOLD3/4 patients were not revealing any enrichment of terms related to lipid metabolism, but rather terms associated with cell cycle and the secretion of pro-inflammatory cytokines.

The coexpression of genes associated with cholesterol homeostasis were not previously described to be altered in AMs of COPD patients. Therefore, we investigated the GOLD2-associated modules in more detail and identified genes involved in sterol biosynthesis (DHCR7, SQLE), cholesterol transportation to cells (LDLR) and fatty acid metabolism (FADS2) (figure 4d). In contrast, genes of the interferon response (CXCL10, IFIT3, IFITM1, IFITM3, OAS1, RSAD2, and ISG15) (figure 4e) and genes involved in the G2M checkpoint including the proliferation marker MKI67, cyclin-dependent kinases (CDK1) and DNA replication machinery (TOP2A, CENPA) (figure 4f) were highly expressed in GOLD3/4 AMs.

In summary, AMs derived from the BALF of COPD patients show numerous transcriptional changes that suggest broad reprogramming including genes involved in lipid metabolism with different functionalities that follow GOLD grade-specific patterns.

Lipidome analysis of AMs reveals GOLD-specific changes in COPD

The AMs are the major cell type in the alveolar space being involved in lipid metabolism, e.g. of the surfactant [58]. Performing a transcriptome analysis of AMs derived from COPD and control patients using two different mathematical approaches identified reprogramming of lipid metabolism. To elucidate whether the changes we observed in these cells’ transcriptomes in COPD is associated with changes in their lipidome, CD45+Lin−HLA-DR+CD66b−Autofluoreschencehi AMs were sorted from the BALF of eight COPD and seven control patients and were used for mass spectrometry-based lipidomics (figure 5a). Seventeen different lipid classes including a total of 202 lipid species were quantified. Glycerophospholipids, especially the phosphatidylcholines (PCs), alongside phosphatidylethanolamines (PEs), phosphatidylinositols (PIs) and phosphatidylserines (PSs), the sphingolipids (ceramides (Cer) and sphingomyelins (SMs)) and the monoacylglycerols (MAGs) constituted the major lipid mass of AMs (figure 5b). Comparison of the AM lipidome of COPD with that of control patients indicated most prominent changes in the cholesterol ester (CE) and MAG classes (figure 5b) and quantification of the mean GFC revealed disease severity-specific changes in the analysed lipid classes (figure 5c). In accordance with the predicted change in cholesterol metabolism, AMs from GOLD2 COPD patients exhibited an upregulation of CE, but also downregulation of hexosylceramides (HexCer). In contrast, AMs from GOLD3/4 COPD patients showed an increase in MAG and diacylglycerols (DAGs) and a decrease in the PC class. GOLD grade-specific patterns were reflected as well on the single lipid species level, revealing the diverse changes within each lipid class (figure 5d and e), all in all prompting us to investigate this in greater detail.

FIGURE 5

Quantitative lipidomics of alveolar macrophages (AMs) of COPD patients describe disease severity. a) Schematic representation of the analysis workflow for analysis of the 202 lipid species assessed by lipidomics of AMs from bronchoalveolar lavage fluid of 11 COPD and four control patients. b) Mean mol% per patient group for lipid class sums of lipid species normalised to total identified lipid content. c) Hierarchical clustering of average log2 fold change (FC) of mol% of lipid class sums. d) Average log2 FC of mol% for single lipid species per lipid class for Global Initiative for Chronic Obstructive Lung Disease (GOLD) stage 2 versus control or e) GOLD3/4 versus control. Bubble size indicates p-value. f) CoCena2 cluster-condition heat map. g) Hierarchical clustering of the average lipid abundance for lipid species encompassed in the GOLD2-specific modules 2 and 3 and h) GOLD3/4-specific module 8. GFC: group fold change; CE: cholesteryl ester; Cer: ceramide; CL: cardiolipin; DAG: diacylglycerol; DiHexCer: dihexosylceramide; HexCer: hexosylceramide; LPC-O: ether-lysophosphatidylcholine; LPC: lysophosphatidylcholine; MAG: monoacylglycerol; PC-O: ether-phosphatidylcholine; PC: phosphatidylcholine; PE-O: ether-phosphatidylethanolamine; PE: phosphatidylethanolamine; PI: phosphatidylinositol.

For this purpose, we determined lipid-lipid level correlations using the CoCena2 pipeline. We identified nine different modules encompassing similarly regulated lipid species across the disease groups (figure 5f). Modules 1–3 showed enriched lipid species in GOLD2, whereas module 8 encompassed lipid species specifically enriched in GOLD3/4 COPD patients (figure 5f).

A closer look at the GOLD2-specific modules revealed enrichment for many unsaturated PC (e.g. PC(32:2, 36:3, 38:4, 38:5, 40:4) and CE (CE(16:1, 18:0, 18:1, 18:2, 18:3, 20:4)), as well as a few SM (e.g. SM(18:0, 22:1, 24:0,24:2)), lysophosphatidylcholines (LPCs) (e.g. LPC(20:4, 22:5, 22:6)) and Cer, such as Cer(d18:1/24:0) or (d18:1/24:1) (figure 5g). On the other hand, distinct lipid species were correlated in the module specific for GOLD3/4 AMs, which showed an accumulation of neutral lipids, such as saturated MAGs (MAG(16:0, 18:0)) and DAGs (DAG(32:0, 34:0, 36:0)) (figure 5h). Taken together, AMs from COPD patients display changes in their lipidome, which, in addition to their transcriptome, are GOLD grade-associated and could constitute a hallmark in the pathogenesis of COPD.

Discussion

Chronic inflammatory diseases such as COPD are associated with functional and transcriptional alterations of immune cells involved in the pathophysiological processes and disturbances of organ-related homeostatic functions, e.g. surfactant biology in the lung. Here we focused on the most abundant immune cell type in the alveolar space, namely AMs to determine transcriptional changes particularly in those genes associated with lipid metabolism. By defining the AM transcriptomes, we observed numerous GOLD grade-specific changes in COPD patients. Evident changes in gene expression related to lipid homeostasis were further evaluated by lipidome profiling of AMs, which demonstrated that prominent lipid changes correlate with disease severity.

We focused on AMs, the most prevalent immune cell in the alveolar space of COPD patients and investigated their transcriptomes. Previously, AMs were shown to convert to an anti-inflammatory phenotype in COPD smokers based on the simple paradigm of M1, M2 polarisation [47, 49, 54]. Additionally, a recent article found that some of those cells may carry a dual M1, M2 pattern in severe COPD [53]. However, we have previously demonstrated that macrophages follow a multi-dimensional model of activation [30]. Comparison of these various activation signatures with the transcriptome of AMs from COPD patients revealed that these cells can be better described by various activation signatures (including PGE2 (c4, c8) and lipid (c5) signatures). To further explore the overall transcriptional changes in AMs derived from COPD patients, we provide the data for further easy inspection at FASTGenomics.org. Here we focused on the changes related to lipid metabolism.

AMs maintain lung function by supporting the air–liquid interface on the alveolar surface, which, when compromised, results in augmented surface tension and pulmonary proteinosis [16]. It has been shown that lipid homeostasis in COPD is perturbed in both humans [21] and cigarette smoke-induced mouse models [59]. In addition, decreased BALF levels of cholesterol, SM and phospholipids correlate with loss in lung function in COPD patients [17]. Our analysis shows that AMs from GOLD2 COPD patients are enriched in PC and cholesterol esters, and express higher levels of lysosomal lipase LIPA, which suggests that in mild severity stages, AMs may contribute to the decrease in surfactant lipids in the alveolar space by increasing their surfactant catabolism rate. In contrast, in more severe stages, AMs are characterised by the enrichment in saturated MAGs and DAGs, which indicates a degradation defect is very likely due to the accumulation of glycerolipids.

Recently, lipidome analysis of perturbed macrophages showed that changes in lipid composition correlate with functional responses [60]. AMs exhibit a foam-like phenotype with augmented transforming growth factor-β secretion in response to silica dust and cholesterol uptake [61]. Other findings were changes in cholesterol biology. Here one could postulate that the balance between cholesterol synthesis and the liver X receptor (LXR) pathway acts as a rheostat that controls AM inflammatory responses. In GOLD2 patients, cholesterol synthesis was upregulated and LXR signalling was repressed. In contrast, in GOLD3/4 COPD patients, cholesterol synthesis was downregulated, which would be consistent with the upregulation of the LXR pathway and the induction of lipogenesis [62]. This switch could have implications in AM biology in COPD, as LXR activation confers antiviral protection [63] and inflammation resolution [64].

In addition, the observed heightened levels of CE and SM in AMs from GOLD2 COPD patients may be responsible for the dysregulation of AM functionalities, such as efferocytosis of apoptotic cells as shown before in rats exposed to cigarette smoke [65]. Dampened PI levels curtail viral infections in mice via Toll-like receptor pathway inhibition [66]. Our data provide the enriched interferon responses as additional mechanistic insight in AMs from GOLD3/4 patients.

As with every exploratory study, some of the limitations of our study need to be mentioned. First, the control cohort consisted of patients suffering from chronic idiopathic cough, which might be regarded as a potential confounder given the presence of lipid-laden macrophages in human BALF samples [67]. However, the assessment of these patients demonstrated that they suffered from an exquisitely sensitive cough reflex without any underlying lung pathology. On average, the control group was also younger than the COPD group and too small to thoroughly investigate differences between smokers and nonsmokers in COPD and control patients. To alleviate any confounding age, smoking status effects in our transcriptome analysis, we computationally removed technical and unknown variance before downstream analysis. We suggest that any effect of age or smoking on high-dimensional data, including lipidomics might be indeed interesting questions for future studies.

Taken together, although the numbers of AMs are not significantly changed in the alveolar space in COPD, we conclude that they undergo transcriptional reprogramming that leads to GOLD grade-dependent changes in lipid metabolism. Elucidation of the concomitant alterations in the transcriptome and lipidome of AMs aids the understanding of their role in COPD and provides druggable molecular pathways, such as cholesterol metabolism and interferon-α and γ responses.

Please note: supplementary material is not edited by the Editorial Office, and is uploaded as it has been supplied by the author.

Supplementary material 00915-2020.SUPPLEMENT

Figure S1 00915-2020.FIGURES1

Figure S2 00915-2020.FIGURES2

Table E1 00915-2020.TABLEE1

Table E2 00915-2020.TABLEE2

Table E4 00915-2020.TABLEE4

Table E5 00915-2020.TABLEE5

Table E6 00915-2020.TABLEE6

Table E7 00915-2020.TABLEE7