Alveolar macrophages in early stage COPD show functional deviations with properties of impaired immune activation

Despite its high prevalence, the cellular and molecular mechanisms of chronic obstructive pulmonary disease (COPD) are far from being understood. Here, we determine disease-related changes in cellular and molecular compositions within the alveolar space and peripheral blood of a cohort of COPD patients and controls. Myeloid cells were the largest cellular compartment in the alveolar space with invading monocytes and proliferating macrophages elevated in COPD. Modeling cell-to-cell communication, signaling pathway usage, and transcription factor binding predicts TGF-β1 to be a major upstream regulator of transcriptional changes in alveolar macrophages of COPD patients. Functionally, macrophages in COPD showed reduced antigen presentation capacity, accumulation of cholesteryl ester, reduced cellular chemotaxis, and mitochondrial dysfunction, reminiscent of impaired immune activation.


Introduction
Worldwide, chronic obstructive pulmonary disease (COPD) is the third leading cause of death (1,2). Due to smoking and increasing air pollution, the current prevalence of 10.1% is estimated to further increase in the next decades (2). Considering the enormous medical and financial burden of COPD, there is a need to develop efficient biomarker-based diagnostics, and molecularly guided therapies. It is now accepted that COPD is a heterogeneous disease manifesting as a clinical syndrome with structural pulmonary abnormalities, lung function impairment, chronic respiratory symptoms, or any combination of these. Consequently, the pathogenesis of the disease is complex with numerous co-existing mechanisms with inflammation being one of the most prominent and important mechanisms (3). Lung inflammation in COPD is characterized by alterations in the number and function of immune cells. Alveolar macrophages (AMs) are considered to be one of the major orchestrators (4). Yet, little is known about the heterogeneity of AMs in COPD as well as the underlying molecular mechanisms leading to AM alterations, particularly during earlier disease stages.
To characterize molecular and functional alterations in the myeloid compartment in COPD, we here applied single-cell transcriptomics combined with extended data analytics, as well as phenotypic and functional assays to characterize the molecular changes in myeloid cells derived from bronchoalveolar lavage fluid (BALF) and peripheral blood obtained from patients with early-stage COPD (Global Initiative for Chronic Obstructive Lung Disease (GOLD) stage 2).

Heterogeneous cellular states of macrophages in the human alveolar space
We obtained freshly isolated BALF material and peripheral blood ( Figure 1A) from COPD patients and donors with chronic cough, but without any signs for pathophysiological alterations of the lung (hereafter referred to as 'control') (Table S1). We conducted a pilot experiment, in which we obtained single-cell RNA-sequencing (scRNA-seq) data using the most widely used droplet-based solution [Chromium from 10x Genomics (5)] and a well-based method [Seq-Well (6)]. After identification of cell-types based on marker gene expression of defined clusters ( Figure S1A), we compared the cell populations between the two technologies. As ground truth, we characterized the cellular compartment in the alveolar space using multi-color flow cytometry (MCFC) ( Table S2, see methods). All three approaches identified macrophages as the predominant cell type in the alveolar space ( Figure S1B). When determining the cell type distribution for the droplet-and well-based scRNA-seq methods, granulocytes (neutrophils, eosinophils) were almost undetectable in the droplet-based method ( Figure S1B).
Since the population structure in Seq-Well was more closely related to MCFC data, we continued with the wellbased scRNA-seq method and generated 60,925 single-cell transcriptomes from BALF derived from 9 patients with early-stage COPD (GOLD stage 2) and 6 controls, as well as 54,569 single-cell transcriptomes from peripheral blood of 6 B C D E F A FIGURE 1 scRNA-seq data of BALF samples obtained from COPD patients and healthy controls. (A) Schematic workflow of the present study. Bronchoalveolar lavage fluid (BALF) and peripheral blood was obtained from control donors and COPD patients (GOLD stage 2). After enrichment for immune cells (CD45+ cells), single-cell RNA-seq was performed. (B) UMAP representation of integrated BALF data obtained from all COPD patients and control donors. Coloring and numbering according to identified main clusters. (C) Heatmap of the calculated marker genes per main cluster with a bar chart representation of the relative cell type proportions at the top. The marker gene expression per cluster is represented as a z-transformed value (across all clusters). Rows of the heatmap are clustered hierarchically. At the bottom of the plot, the main cell type is displayed, which is contained in the respective main cluster. (D) Schematic workflow of the four-step annotation approach, including machine learning-based cell type annotation, clustering, assignment and subsequent confirmation of a cluster to a cell type according to the machine learning-based cell type annotation, and identification of 'contaminating' cells (referred to as 'mixed cells'). (E) Final cell type annotation of integrated BALF data according to the four-step annotation approach. (F) Volcano plot visualization of log2 fold changes and negative log10 p-values (Wilcoxon rank sum test) of changes in cell type occurrence in BALF of samples obtained from COPD patients and controls. BALF, bronchoalveolar lavage fluid; alv., alveolar; MDM, monocyte-derived macrophage; DC, dendritic cell; n, number; MФ, macrophage.
COPD patients and 4 controls ( Figure 1A; Table S1). Starting with BALF cells, we first used a classical clustering approach and visualized the data via UMAP in 17 clusters ( Figure 1B). By marker gene identification on the majority of cells in each of the clusters, we identified the major cell types present in BALF ( Figure 1C). A more detailed inspection of individual clusters revealed further cluster substructures. To better describe the cellular compartment in BALF we developed and applied a four-step cell type annotation procedure (Figure 1D, S1C-G) (for details see methods section 'four-step cell type annotation'). Macrophages were not only the most prevalent, but also the most heterogeneous class of immune cells in the alveolar space ( Figure 1E, Table S3), but we also identified dendritic cells (DCs), monocytes, neutrophils, eosinophils, mast cells and T cells in BALF which is in line with recent reports (7). Determination of relative frequencies between COPD and control revealed one of the macrophage states (MФ8) to be elevated in COPD, while the majority of cell types and states did not significantly differ between the COPD and control group ( Figure 1F). Collectively, single-cell transcriptomics reveals a heterogeneous landscape of myeloid cells in BALF with slight shifts in cell state distributions between early stage COPD (GOLD stage 2) and controls.

Proliferating and monocyte-like macrophage states are elevated in COPD
To further characterize the most prevalent and heterogeneous cell types in the alveolar space, we subclustered macrophages and monocytes excluding non-immune cells, neutrophils, basophils, eosinophils, mast cells, DCs, and T cells, which resulted in a total of 13 clusters (Figures 2A, B). Except for cluster 10 (monocytes), all other clusters expressed macrophage cell lineage markers (MSR1, MRC1, MARCO). BALF-derived macrophages disp layed remarkable transcriptional plasticity. The MФ8 macrophage state, elevated in COPD, was characterized by proliferation-associated genes (MKI67, TOP2A, and NUSAP1), as well as increased expression of histone genes (HIST1H4C and HIST1H1D) and most of the MФ8 cells were computationally assigned to the G2/M cell cycle phase ( Figure S2A), strongly supporting these cells representing proliferating macrophages. MФ6 and MФ9 macrophage states were highly enriched for major histocompatibility class (MHC) II expression (HLA-DQ and HLA-DR respectively), while the MФ12 cell state carried hemoglobin genes (HBA2, HBA1, and HBB) either due to engulfed erythrocytes or induction of hemoglobin genes in macrophages. FIGURE 2 Exploration of the macrophage and monocyte cell types and states in human BALF. (A) UMAP representation and clustering of cells annotated as monocytes or macrophages by the four-step annotation approach (according to Figure 1). (B) Heatmap of marker genes per macrophage/ monocyte cluster (referred to as 'macrophage states'; according to Figure 2A). The marker gene expression per macrophage state is represented as a z-transformed value across all macrophage states. On the left side of the heatmap, conserved macrophage markers are depicted. Columns and rows of the heatmap are sorted by hierarchical clustering. (C) Box visualization plot (with marked median values) of most significant differences in population sizes within the identified macrophage states between COPD and control (error bars indicating the standard deviation; statistics based on Wilcoxon rank sum test). n, number; MФ, macrophage.

B C A
Except for the macrophage states MФ12, 2 and 11, we did not identify any donor effect ( Figure S2B) with the latter being characterized by CCL5 expression (CCL5+ macrophage state). The MФ5 state exhibited relatively strong expression of the monocyte-associated genes VCAN and S100A8 together with the monocyte attractant CCL2 and the late monocyte-tomacrophage differentiation marker CHIT1 and was therefore designated as 'monocyte-like' (mono-like) macrophages. Furthermore, these cells, together with the proliferating macrophages, exhibited the largest relative increase in population size in COPD ( Figure 2C). The monocyte-like macrophages also shared some markers with the MФ7 cell state, which was additionally high in interferon-response genes (IFIT1 and IFIT2), and MФ3 cells characterized by increased expression of complement components (C1QA-C) and alpha-1antitrypsin (SERPINA1).
Next, we predicted the functions of each macrophage state by gene set variation analysis (GSVA) ( Figure S2C, D), illustrating shared, but also cluster-specific functions. Among the shared terms, we found enrichment of 'antigen presentation', 'endocytosis', 'oxidative phosphorylation' and 'b-oxidation', which reflect some of the basic cellular processes of macrophages in the alveolar space. Intriguingly, the MФ4 cell state revealed a specific enrichment of the mTOR signaling pathway, which was described to be associated with cellular senescence in non-immune cells from the lung (8). This was further corroborated by enrichment analysis of gene sets associated with cellular senescence, namely genes associated with cell aging and mitochondrial functions ( Figure S2E). Furthermore, a senescent molecular phenotype of the MФ4 cell state was supported by downregulation of the genes also downregulated in the recently described IMM-age signature derived from aged immune cells (9 ) (Figure S2E). Collectively, macrophages in BALF exist in numerous different molecular and functional states with proliferating and 'monocyte-like' macrophage states being elevated in COPD.

Altered lipid metabolism and stressed macrophage phenotypes in COPD
To determine overall functional differences between control and COPD based on macrophage state information, we developed 'GO-shuffling' as a GO enrichment approach ( Figure S3A, see methods section 'Gene set distance analysis of annotated cell types' for more detail). This enrichment analysis showed that mainly metabolism-associated terms contributed to the separation of COPD patients from control donors ( Figure 3A). To examine potential COPD-associated changes in metabolism, we applied the Compass algorithm (10) to comprehensively model the metabolic differences between COPD and control macrophage states. The largest differences were found in amino acid and lipid metabolism ( Figure 3B), with an overall higher predicted metabolic activity in COPD samples ( Figure 3C). Among the differential lipid-associated metabolites and reactions, phosphorylation of inositol was most prominent, but we also found altered metabolites and reactions, indicating increased transport (monoacylglycerol), synthesis (phospholipids and cholesterol) and degradation (boxidation) of lipids in COPD macrophages. In concordance with the increase in lipid metabolism in COPD patients (predicted by Compass), we observed an overall higher expression of genes found in lipid-associated gene sets that were contained in the top 1% of the functional gene sets ( Figure  S3B). Among these genes, we found several receptors for cholesterol uptake (CD36, LDLR, MSR1, and TREM2) and genes of cholesterol storage mediated by cholesteryl ester synthesis (ACAT1/2 and SOAT1), but also genes associated with cholesteryl ester hydrolases (LIPA, CES1, and NCEH1) ( Figure S3B, C). Next, we validated the in silico prediction of altered lipid metabolism in COPD by performing lipidomics analyses of 229 lipid species in macrophages obtained either from COPD GOLD 2 patients or control donors. We observed the greatest difference in the lipid class of cholesteryl esters, which was significantly higher in COPD macrophages than in controls ( Figure 3D, E). These findings indicate that the macrophages in COPD patients show a pulmonary foam celllike response, which has been reported for other lung diseases, such as pulmonary alveolar proteinosis (11). This cellular phenotype is characterized by the cells being predominantly cholesterol-laden. The accumulation of cholesterol in macrophages of pulmonary alveolar proteinosis patients has been associated with downregulated expression of the cholesterol transporter ABCG1 (12). Furthermore, the accumulation of cholesteryl ester has been described in microglia as a consequence of deficient TREM2 signaling (13). Surprisingly, we found an upregulation of ABCG1 and TREM2 expression in COPD macrophages ( Figure S3B). However, NOTCH signaling was also predicted as a strong separator of COPD and control macrophages ( Figure 3A) and, as a consequence, might result in perturbed TREM2 signaling. Investigation of this gene set revealed increased expression levels of the metalloprotease-disintegrins ADAM10 and ADAM17 and the g-secretase component APH1A ( Figure  S3D). These enzymes can cleave TREM2 from the surface and thus interfere with the downstream signal transmission (14). It is possible that elevated TREM2 expression in macrophages is a consequence of COPD-mediated tissue damage and thus increased cellular stress.
Since both increased metabolic activity ( Figure 3C) and putative cell stress demand high amounts of energy, we hypothesized that energy turnover might be increased in m ac r o p h a g e s fr om C O P D p a t i e n t s a n d t h e r e f o r e investigated the mitochondrial function of AMs. In three COPD patients and 2 control donors, we were able to isolate Frontiers in Immunology frontiersin.org sufficient numbers of viable cells to measure mitochondrial function. Indeed, we observed an increased baseline respiration rate in macrophages derived from COPD patients ( Figure 3F, G, Figure S3E), which reflects an elevated energy demand. In line with previous reports (15), we found a significant increase in proton leakage in COPD macrophages, despite similar levels of ATP production, which is indicative for mitochondrial dysfunction and increased ROS production in COPD (16). Reduction of chemotaxis was also predicted for COPD macrophages ( Figure 3A). While CCL3 was elevated in BALF from COPD patients ( Figure 3H), the chemotaxis of COPD macrophages towards CCL3 was reduced ( Figure 3I), indicating that single-cell transcriptomes indeed correctly predicted macrophage function, while elevated chemokine levels in BALF did not serve as a surrogate for cellular function. Taken together, the heterogeneous landscape of BALF-derived macrophages is linked to numerous molecular and cellular alterations in COPD, of which we highlight metabolic and chemotactic changes together with evidence of pronounced cellular stress.

COPD leads to downregulation of MHC expression
We next intended to determine differential gene expression across different macrophage states between COPD and controls. Here, we applied an approach, which includes patient information by testing all possible pairs of patients and controls followed by utilizing the median Wilcoxon score of the pairwise tests as a test statistic ( Figure  S4A-D, for more detail see methods section 'Distribution-free DE analysis across patient groups'). Visualization of the DE genes per macrophage state shows that the majority of the observed transcriptional differences are macrophage state-specific ( Figure 4A), albeit trends for differential expression in the same direction were often seen for other macrophage states as well ( Figure 4B). Interestingly, transcriptional differences are mainly attributable to increased expression in COPD.
In accordance with the Compass analysis ( Figure 3), lipid metabolism-associated genes (e.g. CD36, COLEC12, SOAT1, and PPARG) were identified to be upregulated in COPD ( Figure 4B). Further, metalloprotease-disintegrins ADAM9, ADAM10 and ADAM17, as well as the surface molecule CD163 were elevated across many macrophage states in COPD, which corroborates earlier findings for CD163 by immunohistochemistry (18) ( Figure 4B). Gene set enrichment analysis (GSEA) revealed terms associated with focal adhesion and antigen processing and presentation ( Figure 4C).
When plotting the expression of the top expressed MHC class I-encoding genes (HLA-A, HLA-B, HLA-C, HLA-E) ( Figure 4D) and MHC II-encoding genes (HLA-DRA, HLA-DRB1, HLA-DRB5, HLA-DPA1, HLA-DPB1, and HLA-DQB1) ( Figure S4E), we found these genes largely to be downregulated in COPD. We identified similar downregulation of MHCencoding gene expression in bulk transcriptome data (17) comparing BALF-derived macrophages from healthy donors, healthy smokers, and COPD patients ( Figures 4E, F). Downregulation of MHC molecules was most pronounced in COPD and thus not solely due to smoking. Next, we isolated BALF macrophages from additional patients and measured surface protein levels of MHC class I (HLA-A/-B/-C) ( Figures 4G, H) and class II (HLA-DR) ( Figure S4F, G). MHC class I was significantly reduced on macrophages derived from COPD patients, while MHC class II molecules only showed a trend towards lower expression. In summary, DEG expression analysis revealed significant transcriptional changes in macrophages, including the downregulation of MHC Iencoding gene expressions, which was also apparent on protein level. Modeling of the metabolic landscape and alterations in macrophages. (A) Word cloud of the most common words in the top predicted terms of the GOshuffling approach across all macrophage clusters. (B) Compass results of the modeled metabolic landscape in macrophages. The pie chart summarizes and categorizes the predicted metabolites and pathways that are significantly different between COPD and control. (C) Heatmap showing the metabolites and pathways that were predicted by Compass as altered in COPD and that were associated with lipid metabolism. Metabolites are shown in black and reactions in red. Columns and rows of the heat map are sorted by hierarchical clustering. (D) Volcano plot visualization of log2 fold changes and negative log10 p-values (Wilcoxon rank sum test) of lipid class levels between COPD and control macrophages obtained by lipidomics analysis. (E) Box plot with marked median values of cholesteryl ester proportions with the representation of individual donors. (F) Evaluation of mitochondrial function via the timedependent course of the oxygen consumption rate (OCR) in macrophages using baseline-corrected values. Error bars indicate the standard deviations (control n = 2, COPD n = 3). Dashed arrows represent the injection of various compounds (shown at the top of the plot) used to assess different aspects of mitochondrial function (according to Figure S3E). (G) Bar plots showing quantifications of different aspects of mitochondrial function inferred from the OCR measurement in Figure 3F (according to Figure S3E; error bars indicating the standard deviation; statistics based on t-test). (H) Heatmap representation of proteins detected in BALF with a p-value < 0.1 according to the Wilcoxon rank sum test between COPD patients and control donors (control n = 11, COPD n = 12). The mean protein expression (identified by Olink Proteomics) per donor is represented as a z-transformed value (across all donors). Columns of the heatmap are sorted by hierarchical clustering. (I) Quantification of the migratory capability of macrophages towards CCL3 displayed in a box plot with marked median values and the representation of individual donors (control n = 4, COPD n = 4; error bars indicating the standard deviation; statistics based on t-test). BALF, bronchoalveolar lavage fluid; OCR, oxygen consumption rate.

Cell-to-cell communication via TGF-b signaling explains changes in macrophage states
To define potential upstream regulators for changes observed in COPD, we focused on those macrophage states with a minimum of 30 DE genes between COPD and control ( Table S5). Representation of predicted transcriptional regulators in an UpSet plot showed that YY1, which is an important modulator of TGF-b1 and NOTCH signaling, was the only predicted transcription factor (TF) shared by all macrophage states included in the analysis ( Figure 5A). Elevated TGF-b signaling was further supported by the identification of the TFs TFE3 and MYOD1 with co-regulation being present in monocyte-like macrophages (cluster MФ5) and C1Q + macrophages (cluster MФ3) which was similarly true for the NOTCH signaling related TFs HES1 and HEY1. Other predicted signaling cascades included WNT signaling (e.g. TCF3/4, MYC and NFATC1/3) and TNF/NF-kB signaling (e.g. CEBPB and REL). These major pathways suggested that signals from the microenvironment are important drivers for transcriptional alterations in macrophages. We next applied CellPhoneDB, which models cell-to-cell communications based on known receptor-ligand interactions (19). Network construction of cell-to-cell interactions within control samples revealed monocyte-like and C1Q + macrophages to be the major network hubs ( Figure 5B). In COPD, cell-to-cell communication was increased, which was particularly obvious for C1Q + and monocyte-like macrophages ( Figure 5B). Among the predicted monocyte-like macrophage interactions, which showed the strongest difference between COPD and the control, we identified several receptor-ligand pairs associated with the TNF superfamily ( Figure S5A). Furthermore, we found an increased likelihood of interaction between the ligand TGF-b1 and the receptor TGFBR1 in COPD.
To corroborate this model, we applied NicheNet (20) to monocyte-like and C1Q + macrophages exhibiting the most cellto-cell-interactions ( Figure 5B) and most DE genes in COPD ( Figure 4B). Ligand activity analysis allowed selection of the top 3 ligands that best predicted DE genes in one of the two macrophage states ( Figure 5C, Figure S5B, C). TGFB1, PTPRM and PSAP were predicted to regulate monocyte-like cells, while C1Q + macrophages were influenced mainly by INHBA and to a lesser extent ADAM12. As INHBA is part of the TGF-b superfamily and shares the same signaling cascade via SMAD2/3/4/7, there might be more commonality in ligand activity within these two macrophage states. Most of the genes, for which expression is predicted to be regulated by the aforementioned ligands, showed a clear DE pattern between COPD patients in both C1Q + and monocyte-like macrophages ( Figure 5C), but only weak expression in the other macrophage states ( Figure S5D). In contrast, visualization of the expression of the predicted ligands across the different immune cell types from BALF revealed no clear differences between COPD versus control cells for INHBA, PSAP and ADAM12 ( Figure 5D). If these genes play a role in COPD, the major sources might be cells not present in BALF. For instance, Activin-A, whose subunit is encoded by INHBA, is known to be upregulated on lung epithelial cells from COPD patients (21). However, for the ligands CSF1, PTPRM, and TGFB1, we found a direct link between their ligand activity and the gene expression in BALF cells from COPD patients ( Figure 5D). Since TGFB1 was both predicted as a signaling pathway of transcriptional regulation ( Figure 5A) and identified as a cell-to-cell-interaction partner for monocyte-like macrophages by CellPhoneDB ( Figure S5A), we focused further analysis on this ligand. TGFB1 is upregulated in COPD patients in eosinophils, C1Q + macrophages, monocytelike macrophages, neutrophils and mast cells ( Figure 5D). To assess whether the increase in TGFB1 expression is translated into elevated protein levels, we examined the BALF of COPD patients and control donors for the latency-associated peptide TGF-b1 (LAP TGF-b1), which serves as a surrogate for TGF-b1 protein levels. This analysis showed a tendency towards increased LAP TGF-b1 levels in COPD ( Figure 5E), which is further supported by reports on elevated TGFB1 levels in peripheral lung tissue from COPD patients (22).
In addition to elevated TGFB1 expression in COPD, the receptors with the highest predicted interaction potential score for TGFB1 (TGFBR1 and TGFBR2) exhibited also higher expression in monocyte-like macrophages from COPD patients ( Figure 5F). Further, we visualized NicheNetpredicted signaling and transcriptional regulation events between TGFB1 and its putative target genes shown to be DE in COPD ( Figure 5F). The nodes in the constructed path were colored according to the expression fold change between COPD and control. Among the transcriptional regulators were the classical TGF-b signaling mediators SMAD3 and SMAD4, with SMAD4 showing increased expression in COPD ( Figure 5F).
Finally, further support for the importance of TGF-b signaling in COPD came from elevated expression of genes within the TGF-b signaling cascade in COPD patients but not smokers when compared to healthy non-smokers, as assessed in the dataset from Shaykhiev et al. (17) (Figure 5G). In summary, we predicted TGF-b signaling to be a prominent regulator of gene expression in BALF-derived macrophages in the context of COPD.

The macrophage pool is supplied by blood monocytes in COPD
We predicted TGF-b1 as an important regulator of monocyte-like macrophages ( Figure 5C), and this cytokine has recently been identified as a crucial cytokine in macrophage differentiation (23). Tissue macrophage replenishment is linked to the local proliferation of tissue-resident cells, but also influx and subsequent differentiation of monocyte-derived cells from the circulation (24). The monocyte-like macrophages had transcriptional similarities to monocytes ( Figure 4B) and, at the same time, their expression profile was regulated by typical signaling pathways of cell differentiation ( Figure 5A, C), suggesting that they may be derived from monocytes. To investigate whether the monocyte-like macrophage state represents an early stage of monocyte-to-macrophage differentiation, we used a gene signature of murine monocytederived macrophages (MDM) from the lungs of smoke-exposed mice (Wohnhaas et al., unpublished data) and assessed the enrichment of orthologous genes in the human macrophage states ( Figure 6A). The strongest enrichment of the MDM signature was found in monocyte-like (cluster MФ5) and C1Q + macrophages (cluster MФ3). Utilizing orthologous gene signatures derived from murine lipid-associated macrophages (LAMs), which were shown to be monocyte-derived by lineage tracing (25), also revealed the strongest enrichment in monocyte-like and C1Q + macrophages ( Figure 6A). These enrichment analyses supported the hypothesis that monocytelike, but also C1Q + macrophages, were derived from monocytes.
To establish a direct link from circulating monocytes to the monocyte-related macrophages, we performed scRNA-seq of blood immune cells ( Figure 1, Figure 6B) from the same donors from whom the scRNA-seq data of alveolar space immune cells were obtained. Application of the four-step cell type annotation approach ( Figure 1B) identified the three known blood monocyte populations comprising classical monocytes (CD14 + monocytes), intermediate monocytes (CD14 + CD16 + monocytes) and non-classical monocytes (CD16 + monocytes) along with a small monocyte population that expressed high numbers of interferon-associated genes (IFIT + monocytes) ( Figure 6B). We next described the relationship between blood-derived monocytes and alveolar space-derived monocytes and macrophages by building a model to determine which of the monocyte subtypes in the blood would most likely give rise to the monocyte-like macrophage state. For this purpose, we combined the blood and BALF data while considering donor batches. While this approach enabled the combination of the blood and alveolar space data, we observed a reduced resolution of the defined macrophage states and therefore continued with a simplified annotation for the analysis of the embedded data ( Figure 6C, Figure S6A). Projection of RNA velocity vectors calculated by the scVelo method (26) in a batch-corrected manner onto the embedded data ( Figure S6B) and inference of the main average vector flow visualized by velocity streamlines ( Figure 6C) revealed a clear motion of blood monocytes towards the macrophages, further supporting circulating monocytes to be precursors of macrophages in the alveolar space. Since RNA velocity visualization on the UMAP did not reveal a clear link between individual macrophage states and blood monocyte subsets, we calculated a higher-order representation using partition-based graph abstraction analysis (PAGA) (27) ( Figure 6D). The strongest connection was derived between blood monocytes and monocytes identified in the alveolar space. To evaluate the connectivity of the PAGA network more precisely, we used the connectivity matrix as a test statistic to define the highest likelihood for each of the blood monocyte subtypes to be related to the different macrophage states in the alveolar space ( Figure 6D). The monocytes within the alveolar space served as positive controls indicating very high relationships. However, among the macrophage states, we could establish the strongest connections between the CD16 + monocyte subtype in blood and the monocyte-like macrophages in the alveolar space, further supporting that the monocyte-like macrophages are most likely an early functional state of macrophages after circulating monocytes enter this tissue compartment.
Lastly, we investigated whether the DE genes in macrophages from COPD patients ( Figure 4B, C) were already altered in blood monocytes. For this purpose, we used the DE genes as signatures of up-and downregulated genes. Clearly, these signatures were altered in the different blood monocyte subtypes derived from COPD patients with CD14 + CD16 + and CD16 + monocyte subtypes showing the strongest enrichment of macrophage DE genes upregulated in COPD ( Figure 6E). Of particular interest, MHC class I and II genes were found to be expressed at lower levels in COPD-derived monocytes, supporting a systemic component of COPD leading to transcriptional changes in circulating monocytes ( Figure 6F).
In summary, we provide evidence that blood monocytes contribute to the macrophage pool, with monocyte-like macrophages providing a link between blood and lung. The monocyte-like macrophages are elevated in the alveolar space of COPD patients ( Figure 2C), suggesting an increased infiltration of blood monocytes. In addition, blood monocytes already show transcriptional changes reminiscent of those observed in cells from the alveolar space strongly arguing for a systemic component in COPD.

Discussion
COPD is an inflammatory lung disease with a high global burden, increasing incidence, prevalence, morbidity and mortality, mainly due to rising air pollution and high smoking rates worldwide (2). Yet, the cellular and molecular mechanisms of this heterogeneous disease are far from being fully understood. Not surprisingly, the diagnosis of COPD is solely based on clinical parameters due to the lack of molecularly defined biomarkers and, as a consequence, causal therapies are lacking because of an incomplete understanding of the complex pathophysiology. Assessing the relationship between blood monocytes and BALF macrophages. (A) Violin plots (with marked median values) displaying enrichment of human orthologues of murine monocyte-derived macrophage signature genes across macrophage states in COPD and control based on Area Under the Curve (AUC). (B) Integrated scRNA-seq data of blood immune cells annotated according to the four-step annotation approach (according to Figure 1B). (C) UMAP of embedded macrophages/monocytes from BALF and blood monocytes. Inferred main average vector flow is indicated by velocity streamlines that are projected as vectors. Locations of the main cell types (acc. to the combined labels from Figure S6A) in the UMAP are indicated by the heat maps at the bottom. (D) PAGA graph derived from embedded BALF and blood data (according to Figure 6C). The weight of an edge, which reflects a statistical measure of connectivity, is represented as the edge width. The table below summarizes the results of the PAGA connectivity calculation, where a value of 1 indicates a strong connection and 0 indicates a weak connection between two cell types. Here, we characterized COPD-associated changes in immune cells from BALF and blood using scRNA-seq in combination with the application of advanced computational approaches. Focusing on alveolar macrophages, the most prevalent cell compartment in BALF, we found specific alterations in lipid metabolism, reduced expression of MHC class I molecules, and identified TGF-b1 as a major factor responsible for transcriptional reprogramming in COPD. Overall, our results indicate stressed and dysfunctional macrophages in COPD. Changes of the molecular phenotype were further supported by functional analysis, illustrating mitochondrial leakage and reduced chemotaxis. In addition, proliferating and monocyte-like macrophages were elevated in COPD, with evidence that the latter were derived from blood monocytes ( Figure 7).
Recently, it has been hypothesized that reprogramming of disease-related cells as a potential therapeutic option might only be possible at earlier stages (28). We therefore focused on patients diagnosed with early clinical stage disease (GOLD stage 2). Single-cell transcriptomes from BALF showed many different cellular states within the myeloid compartment, both in COPD patients and controls. We identified numerous alterations, both cell state-specific but also myeloid compartment-wide changes between COPD and controls. Of particular interest is the identification of reduced expression of MHC class I molecules across macrophage states in COPD. This finding is in accordance with previous studies linking downregulation of surface MHC class I in COPD with impaired immunoproteasome activity (29). Macrophages expressing low-level MHC class I are less efficient in inducing antiviral immune responses, which may explain the high susceptibility of COPD patients to viral infections, one of the main reasons for disease exacerbations.
To understand the regulation of the DE genes, we performed transcription factor binding prediction, receptor-ligand interaction modeling, and downstream transcriptional signature prediction, which indicated TGF-b signaling followed by NOTCH-, WNT-, and TNF-signaling to be elevated in COPD. The predicted pathways might also be involved in immunosenescence in COPD. For example, TGF-b1 can signal via the mTOR pathway, which was recently associated with cellular senescence in lung cells (30). As COPD develops preferentially in elderly people who often suffer from several comorbidities, cellular aging has been suggested as a hallmark of the disease (31). Features of cellular senescence comprise an increase in the number of mitochondria and mitochondrial dysfunction, which is reflected by increased proton leakage and an associated increase in reactive oxygen FIGURE 7 Schematic representation of the key findings of the present study In healthy lungs, alveolar macrophages survey the alveoli and remove pathogens and debris to enable proper gas exchange. In the alveoli of COPD patients, the alveolar macrophages accumulate cholesteryl esters. In addition, blood monocytes invade the alveoli and differentiate into alveolar macrophages. The transcriptome of COPD alveolar macrophages indicate TGF-b1-associated cell signaling especially in the early stages of monocyte-to-macrophage differentiation. The alveolar macrophages in COPD show a reduced ability to migrate towards chemokine. Furthermore, they express fewer MHC molecules; especially MHC class I. Together with the reduced phagocytosis of alveolar macrophages in COPD, the ability of these cells for immune surveillance is severely limited during the disease. In addition, their mitochondria are leaking (e.g. to protons) and therefore produce high amounts of reactive oxygen species. Taken together, the guardians of normal lung function (alveolar macrophages) are severely altered in COPD, preventing them from fulfilling their important physiological functions properly. Furthermore, the observation of reduced MHC expression in blood monocytes indicates that the manifestation of COPD has a strong systemic component. species (ROS) production (16). Oxidative stress due to increased ROS production is a feature of COPD and there is evidence that this is partly due to mitochondrial dysfunction (32). In line with cellular senescence, we found increased proton leakage in mitochondria of macrophages from COPD patients. Additionally, the reduced chemotactic capacity of macrophages in COPD might also be a result of aged immune cells (33). Reduced migratory capacity of macrophages can have deleterious consequences for the lung, as it reduces the efficient removal of pollutants from the alveolar space, which can lead to cell death and the induction of inflammation. Moreover, the clearance of the alveolar space is further deteriorated due to decreased phagocytosis in macrophages in COPD (34).
TGF-b signaling can induce downregulation of MHC expression. This effect has been associated with signaling via SMAD4 (35), which gives a direct link between the predicted intercellular signaling pathways and the DE genes observed between COPD and control patients. Moreover, TGF-b1 is a known inducer of ADAM10 and ADAM17 expression (36, 37) and is described to be essential for macrophage homeostasis and the differentiation of monocytes into macrophages (23). Following up on this idea, we performed PAGA and RNA velocity analysis of cells from the peripheral blood and the alveolar space. This model suggested that a proportion of the macrophage pool is replenished from the systemic monocyte pool circulating in peripheral blood. A recently proposed model (24) suggested that under homeostatic conditions survival of tissue-resident macrophages is supported by self-renewal within the local microenvironment while monocyte recruitment is rather limited. During inflammation, tissue-resident macrophages retain the ability to self-renew, but at the same time many blood-derived monocytes are recruited (38). In COPD, we found evidence for both, local proliferation of some macrophages and recruitment of blood monocytes. Indeed, monocyte-like and C1Q + macrophages exhibited strong enrichment of monocyte-derived macrophage signatures (25). Further, RNA velocity analysis supported a differentiation process from blood monocytes, particularly the CD16 + subset towards the monocyte-like cell state within the alveolar macrophage compartment, which is in line with previous findings demonstrating that the murine counterpart of human CD16 + monocytes can differentiate into lung macrophages (39).

Limitations of this study
The analysis of high quality-biosamples is an important prerequisite for high-resolution analysis such as single cell transcriptomes. While we screened many more patients, only a subfraction of BALF samples was of sufficiently high quality for further analyses. In addition, since the beginning of the pandemic we were not able to obtain further BALF samples due to hospital restrictions. While the study was comparable in size to many single cell transcriptome studies prior the pandemic, the last two years have seen an explosion of larger studies, in particular related to COVID-19. As a consequence, the size of the study now appears rather small. Yet, we have identified several important biological findings that characterize early-stage COPD. We anticipate that our study will provide a framework for further functional studies on the immune compartment in COPD, with a particular emphasis on metabolism, but also to further understand the patient heterogeneity we observed for some of the functional outcomes. In this context, it might be of particular interest that the blood compartment also showed already alterations that might be more easily assessed in future studies due to the easier access to this tissue compartment.
We clearly could show that COPD-related signatures derived from BALF-derived macrophages were already enriched in the peripheral blood monocyte pool, particularly in CD14 -CD16 + and CD14 + CD16 + subsets. These findings indicate that the pathophysiology of COPD is not restricted to the lung. More specifically, reduced MHC expression was also observed on circulating blood monocytes, which further underlines a systemic component of COPD (40). Importantly, elevated levels of TGFB1 have been described in plasma of COPD patients (41) that could explain the low MHC expression in blood monocytes. Finally, as we provide all single-cell transcriptome data and analyses in an integrated fashion on https://www.fastgenomics.org/ ( Figure S7) our data are easily accessible for further analysis.

Contact for reagent and resource sharing
Further information and requests for resources and reagents should be directed to and will be fulfilled by the lead contact, Prof. Dr. Joachim L. Schultze (joachim.schultze@dzne.de). A detailed list of the used reagents and resources is provided in Table S7.

Subject and method details 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. Each individual included in this study was diagnosed and the disease stage was stratified according to the recommendations of the global initiative for chronic obstructive lung disease (GOLD) (COPD recommendations, 2020), with a ratio of post-bronchodilator (salbutamol 400 μg) forced expiratory volume in 1 s (FEV1) to forced vital capacity (FVC) of less than 0.7, and moderate airflow limitation (50% <= FEV1 < 80%). For scRNA-seq, the eligible patients were aged 40 years or older and were either current or ex-smokers. Since COPD has recently been suggested to be a clinical syndrome rather than a single disease (42), we anticipated that despite the focus on GOLD 2 patients, the current study should include a spectrum of COPD patients (Table S1). For example, the generated dataset comprised COPD GOLD 2 patients with different emphysema proportions, exacerbation histories and even a patient suffering from combined pulmonary fibrosis and emphysema (CPFE). The latter patient was admitted based on an external diagnosis of COPD that was later diagnosed as CPFE. This disease type was first described by Cottin et al. (43) and is defined radiologically by the presence of classical features of emphysema in the upper lobes and pulmonary fibrosis in the lower lobes and subnormal lung volumes and severe reduction of CO transfer. Irrespective of the expected heterogeneity within the COPD GOLD 2 patient cohort, stringent exclusion criteria for the current study were a primary diagnosis of asthma with a physician-judged need for oral corticosteroid therapy, clinically significant cardiovascular disorders or laboratory abnormalities and unstable concurrent disease (e.g. exacerbation of disease) that could have affected safety (as judged by the investigator). Individuals suffering from chronic cough without any signs of severe lung pathophysiology or subnormal lung functions served as control donors.

Isolation of cells from bronchoalveolar lavage fluid
Human BALF was obtained from patients with or without COPD via bronchoscopy (at the University hospital Bonn). BALF was performed according to the official American Thoracic Society guideline for interstitial lung disease patients to ensure highest quality of biospecimen material (44). According to these guidelines, we excluded more than half of the clinical samples from further analyses because either the volume of saline solution recovered compared to the amount previously injected into the lungs during bronchoscopy was too low, or blood contamination or increased upper respiratory secretion was present. Each of these factors has an influence on the differential cell count of BALF samples and would have therefore had a negative effect on the analysis results. BALF samples fulfilling the quality criteria were once washed with PBS supplemented with 1 mM EDTA followed by washing with PBS supplemented with 2% fetal calf serum (FCS) and 1 mM EDTA. Throughout the isolation process, the samples were kept at 4°C and centrifugation steps performed at 300 g for 10 min. To exclude any macroscopic non-cellular particles and nonimmune cells from further analyses, immune cells were enriched with MACS columns by using CD45 microbeads according to manufacturer's instructions.

Isolation of peripheral blood mononuclear cells and granulocytes
For the assessment of relationship analysis of the myeloid cell compartment in BALF with cells from the systemic circulation, we obtained venipuncture blood on the day of bronchoscopy. PBMC were obtained by Pancoll density centrifugation (at 20°C and 700 g for 25min with centrifugation break was turned off) of the peripheral blood. After harvesting PBMC from the interphase, all further steps were conducted at 4°C. Granulocytes were recovered from the granulocyte/erythrocyte fraction using cold ACK (ammonium chloride potassium) lysing buffer (1.5M NH 4 Cl, 0.1M KHCO 3 and 1mM EDTA in H 2 O with pH 7.4 at 8°C) to lyse erythrocytes, followed by a washing step with PBS supplemented with 2% FCS and 1 mM EDTA. All centrifugation steps required for granulocyte isolation were performed with max. 300 g for 10 min. To assess the granulocyte fraction in further analyses (particularly in scRNA-seq experiments, Table S1), it was mixed with the PBMC fraction in the ratio PBMC:granulocytes = 2:1. Finally, the PBMC/granulocyte mix was stained with CD45 microbeads for 15 min in order to use a magnetic field in the cell loading of Seq-Well arrays (see below). This artificial ratio allowed to assess the granulocytes in addition to the PBMCs without sequencing the majority of blood immune cells being granulocytes allowing sufficient granularity in the PBMC fraction.

Flow cytometric data generation
Cells were resuspended in PBS supplemented with 2% FCS and 1 mM EDTA for surface marker staining (Table S2). To distinguish live from dead cells, the cells were incubated with LIVE/DEAD Fixable Yellow Dead Cell Stain Kit (1:1000) at room temperature for 15 min protected from light. After washing, human FcR blocking reagent was included to reduce unspecific staining (incubation on ice for 15 min). Next, surface antibodies were added and after 30 min incubation at 4°C in the dark, cells were washed and analyzed either on BD FACSAria III (Becton Dickinson; 3 lasers: violet, blue, and red) for acquisition and sorting or on BD FACSCanto II (Becton Dickinson; 2 lasers: blue and red) for acquisition only (Table S2). Fluorescence-minus-one (FMO) controls were prepared for non-lineage markers.
We performed differential marker intensity measurements across individuals based on the Cohen's d definition of effect size as follows: with FMO = fluorescence minus one and SD = standard deviation. This procedure was followed since we observed strong variability in autofluorescence intensities of macrophages among donors, despite strictest standard operating procedure (SOP) compliance and the use of SOPs for application settings (50) during flow cytometry to minimize potential biases that can occur during sample-to-sample flow cytometry comparisons.

MitoStress assay on seahorse
For the analysis of the metabolic state of donor-derived alveolar macrophages, freshly obtained BALF was centrifuged for 10 min at 300 g. Cell pellet was then washed carefully in PBS (supplemented with 0.02% EDTA) and finally resuspended in MACS buffer. Cell suspension was then stained for 15 min with CD66b microbeads and depleted from granulocytes according to manufacturer instructions. Granulocyte-depleted cell suspension was counted and seeded in Seahorse XF RPMI medium (supplemented with 2 mM L-glutamine, 1 mM sodium pyruvate, 10 mM glucose, adjusted to pH 7.4 prior to the assay) at a concentration of 200,000 cell per well; for each sample, 2 to 4 technical replicates were performed. Cells were then incubated for 30 min in a 37°C incubator, washed two times with pre-warmed Seahorse XF RPMI medium to remove all nonadherent cells and loaded onto the Seahorse XFe96 Analyzer (Agilent). After 3 cycles of baseline measurement, whereby one cycle is defined as 3 min of initial mixing and 3 min measurement, the cells were subsequently injected with Oligomycin (1:1000), FCCP (1:500) and finally a combination of Antimycin A and Rotenone (both 1:2000). Following each injection, oxygen consumption rate (OCR) was measured for 3 cycles.
After the assay, the relative cellular number was determined via crystal violet staining. Shortly, cells were fixed with 4% PFA for 5 min at room temperature and stained for 30 min with crystal violet (0.05% in H 2 O). After two washes with H 2 O the staining was air dried and the formed crystals were dissolved in 200 μL of methanol. Absorbance at 590 nm was measured and used to normalize the Seahorse assay within the Wave software (Agilent). The normalized data were finally exported, further analyzed and visualized in R, with values adjusted to the measured baseline (baseline-corrected). Basal respiration was calculated as baseline OCR at the beginning of the measurement -(OCR after addition of rotenone + antimycin A), maximal respiration as OCR after addition of FCCP -(OCR after addition of rotenone + antimycin A), and proton leak as OCR after addition of oligomycin -(OCR after addition of rotenone + antimycin A) according to schema in Figure S3E.

Migration assay
Migration was analyzed in 24-well transwell plate containing a 8 μm polycarbonate membrane. Macrophages were first purified by FACS according to the expression of CD45, CD66b, HLA-DR and the absence of CD3, CD19 and CD56. Cells were also selected according to the strong autofluorescence signal (51). Macrophages were cultured in 300 μL starvation medium (RPMI 1640 medium supplemented with 0.5% FCS and 1% penicillin/streptomycin) and 50,000 macrophages were seeded in each upper well, while the lower chamber was filled with 700 μL starvation medium only. After an incubation of 1 h in a 37°C incubator, the medium in the upper chamber was exchanged with 300 μL fresh starvation medium and the medium in the lower chamber with 700 μL starvation medium supplemented with 100 ng/mL recombinant human CCL3. The seeded macrophages were incubated at 37°C overnight. Next, cells on the upper filter surface were removed with a cotton swab. Transmigrated cells on lower filter surface were incubated with 2 μM CFSE in 700 mL PBS for 10 min in a 37°C incubator. The transwell inserts were then transferred into wells containing 700 μL RPMI 1640 medium supplemented with 10% FCS and 1% penicillin/streptomycin and incubated for 10 min in a 37°C incubator. Finally, transwell inserts were washed with PBS and imaging of cells was performed using an inverted fluorescent microscope (Nikon) with a 10-fold objective and GFP filter. The number of migrated cells was quantified using ImageJ [version 2 (52)].

Measurement of proteins in BALF
After isolation of cells (see above), the supernatant of BALF samples of both COPD patients and controls were collected and frozen at −80°C before proteomics measurement. Protein levels from cell-free BALF samples were determined using the INFLAMMATION panel from Olink Proteomics, a commercial multiplex immunoassay for high-throughput detection of 92 inflammation-related protein biomarkers. The obtained normalized results (Table S4) were further analyzed in R, whereby proteins were kept for visualization that showed a statistically significant difference (Wilcoxon rank sum test-based p-value < 0.1) between COPD and control samples. Raw files were converted to mzml files and imported into and analyzed by LipidXplorer (version 1.2.8 (53)) software using custom mfql files to identify sample lipids and internal standards. For further data processing, absolute amounts were calculated using the internal standard intensities followed by normalization of the identified lipids on total lipid content. Lipid class sums were calculated for each donor and log 2 -transformed. Differential lipid classes were calculated between COPD GOLD 2 vs control samples using the 'limma' package [version 3.42.2 (54)] under consideration of 'date of sampling'.

Nanodroplet-based scRNA-seq
For comparison of nanodroplet-based scRNA-seq with array-based scRNA-seq (Seq-Well technology, see below), cell preparations derived from three blood and three BALF donors were split in half to be further processed with the two different scRNA-seq technologies by two teams simultaneously. For each donor, 10,000 BALF or blood-derived cells were loaded onto the Chromium ™ Controller instrument (10x Genomics) using the Chromium ™ Single Cell A Chip Kit together with the Chromium ™ Gel Bead Kit v2 following the manufacturer's recommendations. Libraries were prepared using Chromium ™ Single Cell 3' Library Kit v2 according to manufacturer's recommendations and sequenced paired-end as followed: Read 1 26 cycles, i7 index 8 cycles and Read 2 56 cycles on a NextSeq500 instrument (Illumina) using High Output v2.1 chemistry. Single-cell data was demultiplexed and converted into fastq format using bcl2fastq2 (v2.20).

Preparation of Seq-Well arrays
Seq-Well arrays were prepared as described by Gierahn et al. (6). Briefly, Sylgard base and crosslinker were mixed at 10:1 ratio for 10 min, placed under vacuum pressure for 15 min to remove air bubbles and were next poured for a 2 h incubation at 70°C into a wafer with a mounted 86,000 well pattern-holding microscope slide. The arrays were then removed from the molds, excess silicone was cut off with a blade and were prepared for the functionalization process. This protocol adds chemical moieties to the surface of the arrays which facilitate the sealing of a semi-permeable polycarbonate membrane and the interchange of lysis and RNA hybridization buffers. Arrays were rinsed with EtOH, plasma treated for 10 min and successively submerged in APTES (0.05% APTES in 95% EtOH), acetone and PDITC buffers (0.2% PDITC, 10% pyridine, 90% DMF). Upon further washes with acetone, the arrays were spun and dried at 70°C for 2 h. Among the most critical steps in the protocol was the incubation of the arrays with 0.2% chitosan solution (pH=6.3) at 37°C for 1.5 h, after which an overnight incubation in PGA buffer (20 μg/mL polyglutamic acid, 2 M NaCl, 100 mM sodium carbonate (pH=10)) at room temperature under vacuum pressure followed. Finally, the arrays were removed from the vacuum and were rotated for 3 h at room temperature and subsequently moved to 4°C for at least 24 h before use.

Preparation of Seq-Well libraries and sequencing
Seq-Well libraries were generated as recently described by Gierahn et al. (6). After loading of the functionalized arrays with mRNA capture beads, 20,000 CD45 + cells were applied that were previously coated with CD45 + magnetic beads (see above) and suspended in RPMI 1640 medium supplemented with 10% FCS. During the incubation time of 10 min, the loaded arrays were placed on a strong magnetic plate to support the settling of the cells via a magnetic field. After repetitive washing with PBS and soaking with RPMI 1640 medium, the arrays were sealed using polycarbonate membranes that were 7 min treated with air plasma under mild vacuum (Diener electronic). Following a 30 min incubation time in a 37°C cell culture incubator, the arrays were incubated in lysis buffer (5M guanidine thiocyanate, 1mM EDTA, 0.5% Sarkosyl and 1% b-mercaptoethanol in H 2 O) for 20 min and then placed in hybridization buffer (2M NaCl, 3mM MgCl 2 and 0.5% Tween-20 in PBS) for 40 min. Next, the mRNA capture beads were washed from the arrays and collected using washing buffer (2M NaCl, 3mM MgCl 2 and 20mM Tris-HCl pH 8.0 in H 2 O). The reverse transcription was performed on the bead pellet using a Maxima Reverse Transcriptase reaction (Maxima RT buffer, 4% Ficoll PM-400, 1mM dNTPs, 1U/μL RNase inhibitor, 2.5 μM template switch oligonucleotide (TSO) primer and 10U/μL Maxima Reverse Transcriptase in H 2 O) for 30 min at room temperature followed by 90 min incubation at 52°C with end-over-end rotation. The reaction was stopped by washing the beads with TE buffer (10mM Tris-HCl pH 8.0 and 1mM EDTA in H 2 O) supplemented with 0.1% Tween-20 (TE-TW) and TE buffer supplemented with 0.5% SDS (TE-SDS). After a washing step in 10mM TrisHCl pH 8.0, excess primers were digested in an exonuclease reaction (ExoI buffer and 1U/μL ExoI in H 2 O) for 50 min at 37°C with end-over-end rotation and washed in TE-TW and TE-SDS. Beads were resuspended in 500 μL H 2 O and counted with a Fuchs-Rosenthal cytometer in bead counting solution (10% PEG, 2.5 M NaCl). Pools of 5,000 beads (10 μL) were then added to 40 μL PCR reactions (2X KAPA HiFi Hotstart Readymix and 25 μM SMART PCR primer in H 2 O) for the amplification of reverse transcribed cDNA libraries (95°C for 3 min, 4 cycles of 98°C for 20 s, 65°C for 45 s, 72°C for 3 min, 12 cycles of 98°C for 20 s, 67°C for 20 s, 72°C for 3 min and final extension of 72°C for 5 min). After PCR, 16,000-20,000 beads were combined (thereafter referred to as 'pools') and further processed. The pools were cleaned with 0.6x volumetric ratio AMPure XP beads (5 min incubation with beads, followed by 3 min on the magnet, two washes with 80% EtOH, 5 min dry-out, elution with 13 μL H 2 O for 3 min, followed by 2 min on the magnet for collection of the eluent) and the library integrity was assessed using a High Sensitivity D5000 assay for the Tapestation 4200 (Agilent).
To reduce library costs, we produced homemade Tn5 transposase according to Picelli et al. (55). Briefly, the Tn5 coding sequence (tnpA gene from Escherichia coli, Uniprot accession number: Q46731, residues 1-476) was purchased as a synthesized gene containing the mutations E54K and L372P for hyperactivation of the enzyme. Overhangs with the restriction sites XbaI and SpeI were used for cloning into pTXB1 vector, generating a Tn5-Intein-CBD fusion construct. The Tn5 coding sequence was validated by Sanger sequencing. Next, the pTXB1-Tn5-Mxe-CBD plasmid was transformed into the E.coli strain BL21. Cells were grown in LB media supplemented with ampicillin at 37°C to an OD 600 0.8. The temperature was then lowered to 10°C and protein expression was induced by addition of 0.25 mM IPTG. After incubation at 23°C for 4 h cells were harvested by centrifugation at 15,000 rpm on a JA 25.50 rotor (Beckman) for 20 min at 10°C. The cell pellet was resuspended in running buffer (20 mM Hepes-KOH, 0.8 M NaCl, 1 mM EDTA, 10% glycerol, 0.2% Triton-X 100) supplemented with 1 mM PMSF and disrupted by sonication. After centrifugation of cell debris at 15,000 rpm on a JA 25.50 rotor (Beckman) for 30 min at 10°C, residual nucleic acid contaminations from E.coli were precipitated by dropwise addition of polyethyleneimine pH 7.5 to a final concentration of 0.3%. The lysate was cleared by centrifugation at 12,000 rpm on a JA 25.50 rotor (Beckman) for 10 min at 4°C. Chitin resin (10 mL) was equilibrated with running buffer and then incubated with the prepared lysate for 1 h at 4°C. Beads were washed with 10 column volumes of running buffer. For elution by self-cleavage via the intein-tag, the Tn5-loaded resin was incubated overnight at 4°C in 3 mL elution buffer (20 mM Hepes-KOH, 0.8 M NaCl, 1 mM EDTA, 10% glycerol, 0.2% Triton-X 100, 100 mM DTT), followed by dialysis at 4°C overnight in dialysis buffer (100 mM Hepes-KOH, 0.2 M NaCl, 0.2 mM EDTA, 2 mM DTT, 0.2% Triton-X 100, 20% glycerol). The protein concentration was determined using Bradford Assay. Glycerol was added to a final concentration of 50% to the protein sample.
The cDNA libraries (1 ng) were tagmented with the prepared single-loaded Tn5 transposase in TAPS-DMF buffer (50mM TAPS-NaOH (pH 8.5), 25mM MgCl 2 , 50% DMF in H 2 O) for 10 min at 55°C and the tagmented products were cleaned with the MinElute PCR kit following the manufacturer's instructions. Finally, a master mix was prepared (2X NEBNext High Fidelity PCR Master Mix, 2.5 μM barcoded index primer, 2.5 μM P5-SMART-PCR primer) and added to the samples to attach the Illumina indices to the tagmented products in a PCR reaction (72°C for 5 min, 98°C for 30 s, 15 cycles of 98°C for 10 s, 63°C for 30 s, 72°C for 1 min). The pools were cleaned with 0.8 x volumetric ratio AMPure XP beads, were run with a High Sensitivity DNA5000 assay on a Tapestation 4200 (Agilent), and quantified using the Qubit high-sensitivity dsDNA assay. Seq-Well libraries were equimolarly pooled and clustered at 1.4pM concentration with 10% PhiX using High Output v2.1 chemistry on a NextSeq500 system. Sequencing was performed paired-end as followed: custom Drop-Seq Read 1 primer for 21 cycles, 8 cycles for the i7 index and 61 cycles for Read 2. Single-cell data were demultiplexed using bcl2fastq2 (v2.20).

Processing of scRNA-seq raw data
For preprocessing, the generated fastq files from both Chromium ™ and Seq-Well were loaded into a data preprocessing pipeline (version 0.31, available at https://github. com/Hoohm/dropSeqPipe) that relies on Drop-seq tools provided by the McCarroll lab (56). STAR alignment within the pipeline was performed using the human GENCODE reference genome and transcriptome hg38 release 27 (57). The resulting datasets were imported into R for further analyses.
For datasets for which TSO primers were used based on the Smart-Seq2 protocol, sequences starting with either the sequence 5'-GGG-3', 5'-ATGGG-3' or cell barcodes with a Hamming distance of 1 to 5'-ATGGG-3' were excluded to avoid overlapping cell barcodes that are increased with this TSO primer. All other datasets were generated with the TSO primers as described in the original Seq-Well protocol. Next, datasets were examined for content of mitochondrial ribosomal transcripts. For further downstream analyses, the highly abundant mitochondrial transcripts MT-RNR1 and MT-RNR2 were excluded. The resulting datasets were then imported into the R package 'Seurat' [version 3.0.0 (58)] for downstream analyses.

Quality control of scRNA-seq data
We defined cells and genes to be included for further analyses by the following criteria for each donor separately (1): Only genes that were found in at least 3 cells were kept (2); To retain granulocytes that contain only very limited number of transcripts, a relatively low threshold of 100 expressed genes was used to keep cells for further analyses (3); With regard to the rate of endogenous-to-mitochondrial counts per cell, blood cells with a rate > 5% and lavage cells with a rate >10% were excluded. For the comparison of scRNA-seq methods for clinical applications, these quality control filters resulted in a Chromium  (Table S3).

Comparison of different single cell transcriptome technologies
We conducted a pilot experiment, in which we obtained single-cell RNA-sequencing (scRNA-seq) data using the most widely used droplet-based solution [Chromium from 10x Genomics (5)] and a well-based method [Seq-Well (6)]. After identification of cell-types based on marker gene expression of defined clusters ( Figure S1A), we compared the cell populations between the two technologies. As ground truth, we characterized the cellular compartment in the alveolar space using multi-color flow cytometry (MCFC) ( Table S2). All three approaches identified macrophages as the predominant cell type in the alveolar space ( Figure S1B). When determining the cell type distribution for the dropletand well-based scRNA-seq methods, granulocytes (neutrophils, eosinophils) were almost undetectable in the droplet-based method ( Figure S1B).

Dataset integration and dimensionality reduction of scRNA-seq data
If not stated otherwise, all following steps were conducted using the single-cell analysis pipeline Seurat. To account for variations in sequencing depth across cells, we applied a lognormalization strategy using CPM-normalization with a scale factor of 10,000. Next, the genes with the highest cell-to-cell variability in the dataset were determined by calculating the top 2,000 most variable genes by selecting the 'vst' method of the 'FindVariableFeatures' function in Seurat. For the comparison of scRNA-seq methods, the variable genes were determined separately for each technology, while for the integrated analysis of Seq-Well data from COPD GOLD 2 patients and control donors, variable genes were calculated separately for each donor.
To analyze the data without having any influence of batch effects resulting from either different donors or technologies, an integration approach based on 'anchors' across batches (59) was used to harmonize and integrate the different datasets by using the Seurat implementation with the default settings. After linear transformation of the remaining genes (scaling) to ensure homoscedasticity, the dimensionality of the data was reduced to 30 principal components (PCs) that was used as input for UMAP representation.
Next, doublet cells were identified utilizing the R package 'DoubletFinder' [version 2.0.2 (60)] by using the first 30 principal components of the non-integrated datasets, assuming a doublet formation rate of 10% and leaving all other parameters unaltered. The alleged duplicate cells were not removed from the dataset, but accumulations of these cells were highlighted and named accordingly. This procedure revealed, for example, that none of the identified macrophage states was defined by doublet cells (data not shown).

Clustering of the integrated scRNA-seq datasets
The cellular heterogeneity of the integrated datasets was determined using a shared nearest neighbor (SNN)-graph based clustering algorithm implemented in the Seurat pipeline. For both the BALF and the blood data, we used the first 30 principle components as input and set the resolution to 0.7 and 0.6, respectively. The default setting for number of neighbors were used (k=20).

Cell type annotation based on reference transcriptomic datasets
For cell annotation, we developed a slightly modified Python implementation of SingleR (61) (commit a4afed8, available at https://github.com/dviraran/SingleR) and an additional method called GenSigPro. We explicitly used two different methods with different reference datasets to capture variations in the annotation methods.
To compare and integrate these methods (with varying reference data), we first defined a common cell type standard to which all annotated cell types were matched. This standard, the mapping, and the actual reference data are available at F A S T G e n o m i c s ( h t t p s : / / b e t a . f a s t g e n o m i c s . o r g / p / bassler_scCOPD).
The SingleR method iteratively computes the bivariate correlations between the respective cluster expression vector and the multiple reference gene expression vectors for each cell type based on a set of differentially expressed (DE) genes. In each iteration, every cell type in the reference dataset is assigned a score based on these bivariate correlations with the different reference gene expression vectors of that cell type. The cell type with the lowest score is dropped and the DE genes among the remaining cell types are computed and, based on these genes, the bivariate correlations are computed again. This procedure thereby iteratively reduces the number of cell types until only one best fitting cell type is retained. We reimplemented the SingleR functionality to assign cell types per cluster in Python to use in our framework and in addition to the original algorithm, we included a threshold for the bivariate correlation score based on tests with randomized reference data. This made it possible to label cell clusters as "unknown" if the bivariate correlation score of the best fitting reference cell type was below 0.1 and thus no cell type could be assigned. As a reference for SingleR, we used data from both Blueprint+ENCODE (62, 63) and the Human Primary Cell Atlas (HPCA) (64). In addition to the implementation of the SingleR algorithm in Python, we also modified the reference datasets by reducing the reference to immune cells and lung tissue cells. Furthermore, based on the experimental setting of the reference dataset, we adapted some cell labels, e.g. the neutrophils were divided into mature, immature, and inflammatory neutrophils, whereas the original annotation had designated all these cells as neutrophils.
In order to capture potentially relevant variations in the annotation besides SingleR, we developed the similar but distinct statistical approach GenSigPro (Gene Signature Profiler) and incorporated a further reference dataset. To incorporate additional cell types, not included previously, we used manually curated reference data derived from the leukocyte expression dataset LM22 (65). This reference dataset encompasses one gene expression vector (signature) per cell type. While this made it incompatible to be used with SingleR, it allowed the reference dataset to be used in a multiple regression approach. The GenSigPro method fits a multiple linear regression for each cluster expression vector. The covariates in this regression are the reference expression vectors for each cell type that were obtained from the CIBERSORT algorithm (65).
The more similar the cluster expression vector is to one of the reference expression vectors, the higher the regression coefficient for the respective reference vector. If the highest regression vector is positive and above an uncorrected significance threshold of a = 0.05, the cluster is assigned the respective cell-type label of this reference cell type, otherwise, the cluster is labeled "Unassigned". To calculate the regression vectors, we used the Generalized Linear Model (GLM) with an added intercept from statsmodels [version 0.9.0 (66)] with the Gaussian family and left all parameters at their defaults.
GenSigPro does not alter the gene expression vectors during the annotation process, which is in contrast to SingleR where differentially expressed genes are calculated for each iteration. This is especially wanted for manually curated gene lists, like the used LM22 reference, where only high-confidence genes are included. Furthermore, whereas SingleR iteratively selects bivariate correlations, GenSigPro includes all reference gene vectors in a combined model. This allows us to assess the unique contributions of one reference gene vector over the others.
Using this different approach and manually curated reference data, also created more heterogeneous training data for the final consolidation by machine learning (see below). As reference data, we used a manually curated version of the leukocyte expression dataset LM22 (65), where the neutrophils were subdivided according to their activation state. We calculated the reference expression vectors by running CIBERSORT (version 1.06) on the modified LM22 dataset, leaving the default settings unchanged and setting the option "Filter non-hematopoietic genes from the signature matrix during construction". The obtained signature genes (derived from the calculated support vectors) were almost completely (>99%; data not shown) contained in the signature genes of the original CIBERSORT publication (65). As this reference data only provides a mean expression vector per cell type, it was not suitable to be used with the SingleR approach.
Although both SingleR and GenSigPro can be applied also to vectors of single-cell expressions, we applied it to the mean of expression vectors within a cluster for more robust results. Since both GenSigPro and the modified SingleR are Python implementations, we performed clustering using the Louvainclustering (67) function of Scanpy (68) by setting the number of neighbors to 24 and leaving the remaining parameters unaltered.
To assess the uncertainty of the annotation results, we added bootstrapping to GenSigPro and SingleR. The basic principle of bootstrapping is to create an artificial dataset by sampling subjects, in our case cells, with replacement such that in the resulting artificial dataset some cells will be excluded, whereas others will be included more than once. The analyses are then repeated on multiple of these artificial datasets, resulting in somewhat different results. For robust and certain patterns, different bootstrapped datasets generate similar results, while for random fluctuations different bootstraps result in highly different outcomes. Here, we conducted all cell typing analyses using 100 bootstrapped datasets.

Cell type annotation using machine learning
To aggregate and consolidate the initial cell type annotation, we trained a Gradient Boosting Classifier on the combined data of all datasets to classify each cell into a cell type. Gradient Boosting is a machine learning technique that combines multiple classification trees in order to assign an input to different classes. This method is highly flexible and robust in the classification task and has high predictive power. We used an implementation of the Gradient Boosting algorithm from scikit-learn [version 0.19.1 (69)], the leading machine learning library for Python. For training the model, we used the raw gene expression matrix of each cell as input feature for the classification. We additionally extracted features from the data such as the type of tissue, the number of genes per cell, counts per cell, and the percentage of mitochondrial gene expression per cell. The training target of this model were the three cell type labels from GenSigPro and SingleR (Blueprint+Encode and HPCA). For this, we triplicated the data such that each cell with its feature vector was included three times, each with one label of the three cell-type annotations. Our aim was to apply the classifier to all cells in our data. However, as no distinct training data were available, we conducted a 3-fold cross-validation. In this procedure, two random thirds of a data set were used as training data, and the model assigned cell type names to the remaining cells. Importantly, a cell with all three cell type labels was only assigned either to the test or the training dataset. A major advantage of this machine learning method is that the classifier learns the specific expression profile of cell types and can take any cell type annotation as input, independent of techniques, such as bulk RNA-seq or microarray used as initial cell type annotation reference. In addition, we were able to apply the classifier at the single-cell level instead of the cluster mean expression level and thus achieved a higher resolution to exploit the full potential of scRNA-seq. This also allowed us to detect cell types with very low frequency in individual patients. Normally, these cells might end up in larger clusters with a different cell type and are therefore not detected. For all these reasons, this machine learning-based cell type annotation is unbiased, reliable, reproducible and scalable.
Marker gene identification of scRNA-seq data DE genes between identified cell types/clusters (referred to as marker genes) were defined using a Wilcoxon rank sum test for differential gene expression implemented in Seurat. The significance threshold for marker genes were set to an adjusted p-value smaller than 0.001 and the logarithmic fold change cutoff to at least 0.4. In addition, the detected marker genes should have been expressed in at least 50% of the cells within the respective cell types/clusters. Visualization of the obtained marker genes were mainly done using Seurat functions, such as dot plot representation of cell type-/cluster-specific marker gene expression or heatmap representation of marker genes across single cells. A more global overview of the expression profiles was obtained by calculating the mean expression values of marker genes per clusters, followed by scaling and centering of these values and representing them in a heatmap graph using the R package 'pheatmap' (version 1.0.12, https://CRAN.R-project. org/package=pheatmap), in which the genes were clustered according to the 'ward.D' agglomeration method.
Similar to the clustering and marker gene identification of the complete BALF dataset, we performed the same steps also for the detailed characterization of the macrophage population. In addition, we assessed the reproducibility of the identified clusters in the Seq-Well dataset. For this purpose, we used the BALF cells of the Chromium ™ dataset. To consider possible influences of data integration on cell clustering, we used a different integration method, namely Harmony [version 1.0 (70)]. For determining the similarities between the Chromium ™ and Seq-Well clusters, we calculated marker genes and assessed the overlap of the genes per cluster using the matchScore2 package [version 0.1.0 (71)]. For the majority of clusters, we found strong concordance between the Chromium ™ and Seq-Well clusters (data not shown).

Four-step cell type annotation
For the final cell type annotation of the integrated 61K BALF and the 55K blood dataset we used a four-step strategy for cell annotation and for the identification and finally removal of cells of inferior quality. The steps of the strategy include 1) the machine learning-based classifier, 2) cell clustering, followed by 3) a manual classifier-to-cluster comparison and 4) clusterlevel marker gene analysis, including cleanup.
As the first step, the machine learning-based strategy is used to assign the most likely cell type to each cell in the dataset ( Figure S1E). To determine the validity of this approach, we needed a dataset, for which the ground truth of the cell type is known by a secondary method, e g. flow cytometry data. In order to test the validity of machine learning-based strategy we generated a benchmark dataset. The data were obtained by fluorescence-activated cell sorting of blood-derived immune cells using cell type-specific markers followed by SMART-seq2 single-cell sequencing, which gives flow cytometry (ground truth) and scRNA-seq data information for each cell. In this validation experiment for the computational method only cells were used for which RNA expression values of typical cell markers were available ( Figure S1F). Neither SingleR nor GenSigPro alone were able to correctly annotate all cells within the benchmark dataset due to incomplete cell type annotation within the reference ( Figure S1C) used in these approaches. In contrast, the machine learning-based cell classifier was successful in consolidating the annotation results and thus resolving the different cell types in the blood scRNAseq based benchmark dataset ( Figure S1F). Applying the machine learning-based cell annotation to the integrated BALF dataset revealed all major immune cell types and for some cell types a subset structure ( Figure S1G).
The second step consisted of clustering of the data in 18 main clusters ( Figure 1B), which agreed with the areas that were enriched for distinct cell types predicted by the classifier ( Figure  S1G). However, we also found some cells that were annotated, e.g. as dendritic cells (DCs) (Figure S1H), which scattered away from the other DCs.
In the third step of the cell type annotation procedure, we determined which cell type occurred most frequently per main cluster according to the machine learning-based annotation and compared it with the identified marker genes for each cluster ( Figure 1D) (step four). The application of this approach to all 18 major clusters in the BALF dataset led to a detailed resolution of the immune landscape in the alveolar space and this was similarly achieved for the blood dataset.  (77,78). In addition, we retrieved gene sets from WikiPathways (79). This search strategy resulted in a list of 12,755 gene sets, each containing a unique gene set term and a set of associated gene symbols.
As input, normalized scRNA-seq data was used, in which the cells were annotated according to the four-step cell-type annotation approach described above. Cell types containing at least 10 cells for each patient were retained and genes expressed in less than 5% of the cells in the respective cell type were excluded.
For each of the 12,755 gene sets, the "gene set distance" was calculated as follows for each cell type: Gene sets were taken into account that were present with a minimum of 3 genes. For each gene set, the Euclidean distance between all donors was calculated using the get_dist function from the R package 'factoextra' (version 1.0.5). Next, the mean distance of COPD patients, the mean distance of controls and the overall mean distance was calculated. The "gene set distance" was then defined as the overall mean distance divided by the mean distance of COPD patients plus the mean distance of control patients.
This metric allows to determine for which gene set the quotient takes a value close to or greater than 1, which means that the distance within the groups (COPD (dist COPD ) or control (dist CTRL )) is smaller than the overall distance (dist overall ) and consequently the distance is mainly defined by the difference between the groups. Since the Euclidean distance metric is prone to be affected by outliers in higher dimensions, we also tested this approach by using the Manhattan distance and got comparable results. For each cell type, we ranked the gene sets by their gene set distance. Visualization of the most frequent terms contained in the upper percentile of the predicted gene sets in the macrophage states was performed using the R package 'wordcloud' (version 2.6), in which filler and connective words were excluded. Alternatively, the gene sets in the upper percentile were filtered for association with 'NOTCH' or 'lipidomics' and the expression of the involved genes visualized in a heatmap.

Modeling of metabolic pathways based on scRNA-seq data
The metabolic landscape of macrophage states was modeled using the Compass method [version 0.9.5 (10, 80)] by leaving the standard settings unaltered (model: RECON2 (81); lambda: 0; media: media1, which represents a rich extracellular medium, as defined in the Compass manuscript). As input, we simplified the single-cell data of the macrophages by using the 'applyMicroClustering' function of the R package 'VISION' [version 2.1.0 (82)], resulting in approximately 20 microclusters per patient. Next, we applied Compass to the microclusters for each donor separately. The output tables representing Compass scores for single reactions and synthesis of single metabolites of the individual donors were imported into R. They were concatenated and finally transformed as described in the Compass manuscript, except for disabling the division into meta-reactions. In detail, the concatenated output table x was first negatively log-transformed (y = -log(1+x)), the global minimum value of table y was subtracted from the values (z = y -min(y)) and the resulting table z was then used for further analysis. To determine which reactions and metabolites are significantly different between control donors and COPD patients, with the differences being reproducible in the COPD population, we performed Wilcoxon rank sum tests on Compass scores. We first computed the Wilcoxon p-value for every patient separately against all controls, took the median of these p-values, and kept reactions/metabolites for which -log10 (median p-value) ≥ 2.5. We derived a second list of reactions and metabolites by similarly comparing control donors separately against all patients. The reactions and metabolites that have significant differences are the union of these two lists. Next, we excluded reactions with the lowest confidence score in the metabolic reconstruction (83), i.e., we discarded reactions with a confidence score of 1 and kept confidence scores of 2-4 (as well as 0 which is reserved for unannotated confidence). We also excluded metabolites that localize to cellular compartments other than the cytoplasm [c], extracellular space [e] or mitochondria [m]. Finally, the remaining reactions and metabolites were annotated using the Virtual Metabolic Human (VMH) database (84) and visualized in a heat map.

Cell cycle state analysis of scRNA-Seq data
To categorize the cells within the macrophage states into the respective cell cycle states, we applied the 'CellCycleScoring' function of Seurat and substantiated the results using the 'cyclone' function (85) implemented in the R package 'scran' (version 1.10.2 (86)).

Gene set variation analysis
To predict the functions of the macrophage states, we performed gene set variation analysis (GSVA) (87) by using the R package 'GSVA' (version 1.30.0) and defining 'Poisson' for the non-parametric estimation of the cumulative distribution function of expression levels across donors. For the GSVA input expression table, we calculated the sum of the expression of normalized scRNA-seq data for each patient in any macrophage state. As gene sets we used the gene set collection described in the section 'GO-shuffling' and additionally included the 'ImmuneSigDB' collection of MsigDB, whereby this collection was reduced to gene sets that had one of the following terms in the gene set description: 'Mono', 'Macro', 'MDC', 'MDM', 'Dend' and 'DC'. This resulted in 14,160 gene sets. Similar to GO-shuffling, we filtered this collection for gene sets that were present with a minimum of 3 genes in a respective macrophage state. We applied an additional filter step to increase the stringency of the analysis. Therefore, we retained only gene sets in which the sum of the genes contained in the set were expressed in more than 30% of a macrophage state. The GSVA results per donor were combined for the respective macrophage state using a Borda rank and the top 250 ranked gene sets per subtype were visualized in an UpSet plot using the R package 'UpSetR' [version 1.3.3 (88)].

AUCell for gene set enrichment analysis
Enrichment of gene sets was performed using the 'AUCell' method (89) implemented in the package (version 1.4.1) in R. We set the threshold for the calculation of the area under the curve (AUC) to the top 3% of the ranked genes and normalized the maximum possible AUC to 1. The resulting AUC values were subsequently visualized in a violin plot. For statistical testing, a Dunn's Post-Hoc test using the "dunn.test" R package (version 1.3.5) was performed. Resulting p-values were corrected for multiple testing using the Benjamini-Hochberg method. This approach was used, for example, in Figure 6A to assess the enrichment of monocyte-derived macrophage signature genes provided by Wohnhaas (unpublished results). This signature was obtained from scRNA-seq data of monocyte-derived macrophages that were identified in BALF of a murine 12-week smoke model. Human orthologues (obtained from BioMart [version 2.42.0 (90)] of the murine marker genes were used for the enrichment analysis. In a similar way, we also performed the enrichment of monocyte-

Distribution-free DE analysis across patient groups
To analyze the differences between the patient and control cohort, we employed a distribution-free test that preserves patient and cell information and thus considered possible individual donor effects. In contrast to available methods, it avoids the use of mini-bulk, the pooling of cells from different patients, and distribution assumptions. As input, we use the afore-computed macrophage state information and the normalized (non-integrated) scRNA-seq data.
For each macrophage state, a DE analysis between patient and control cohort was performed. Therefore, donors not possessing cells in a clusterwhich happened in a few casesand genes expressed in less than 10% of cells were disregarded for the analysis of this cluster. For each gene, the differences between all possible pairs of patients and controls were assessed using the non-parametric Wilcoxon rank sum test. To assess the differences between patient groups, the median Wilcoxon score of the pairwise tests was considered as a test statistic.
The Wilcoxon rank sum test was chosen because it does not rely on a specific distribution assumption. This is beneficial as the distribution of single-cell expressions is often skewed or shows multiple modes. Furthermore, benchmarking studies revealed that the Wilcoxon rank sum test performs well for the comparisons between two single cell data sets (92,93).
To assess if the observed value of the test statistic was significant, the probability of observing an equal or more extreme value of the test statistic under the null hypothesis was evaluated. The null hypothesis was that there is no difference between the two groups. The null distribution was evaluated with the permutation test, taking all possible permutations into account. For all permutations the afore-described test statisticthe median Wilcoxon scorewas evaluated. The distribution of the test statistic over all permutations provided the null distribution, since reshuffling of patients should not be significant under the null hypothesis. The p-value for the observed group assignment was then the fraction of permutations that led to an equal or more extreme value of the test statistic than the value of the test statistic of the observed patient arrangement.

Testing/Simulation study of the DE method
The DE analysis method was evaluated using simulation data. A first evaluationdenoted as (I)showed a good detection of differences in distributions across groups (patients and controls) with a similar mean. A second evaluationdenoted as (II)indicated that there is no tendency to false positive discoveries if the distributions across groups are similar.
The simulation study was performed on the basis of the here examined COPD dataset. The number of individuals per group and sample sizes per individual were adopted from the original dataset. Sample sizes (number of cells per patient) were taken from the macrophage clusters 0, 1, 3 (Table S3).
The mean of the read count data per gene per individual is sampled from the same log-normal distribution, to ensure variability between the individuals, For simulating differences in the distributions between the two groups (I), distinct success probability parameters (s 1 , s 2 ) were used. Various combinations of s 1 and s 2 were considered to explore the properties of the methods. For simulating similar distributions (II) the success probability parameters were set to the same value (s 1 =s 2 ). s values were set within the interval [0.1, 0.9].
For each combination of s 1 and s 2 , read count data was simulated for 50 m's per individual, which is in the following called a 'set'. In total, for each s-combination, three sets of read count data were simulated. DE analysis was performed with the proposed DE method and for a comparison with the widely used method edgeR [version 3.28.1 (94)], for each set of simulated read count data.
For the case of different distributions between the groups (I), the false negative rate (FNR) was calculated for each set of simulated read count data while evaluating the percentage of genes with p>0.05. In general, the proposed method identified the differences between the distributions for the two groups (low FNR), whereas edgeR failed to identify these differences for all cases with a FNR of over 90%. If the discrepancies between s 1 and s 2 were sufficiently large, which led to a clear difference between distributions (e.g. combinations 1-5 in Table S5), the proposed method performed well and achieved low FNRs. For similar values of s 1 and s 2 (e.g. combinations 6-8 in Table S5), the differences areas expectedmore difficult to identify and the FNR increases. The combination 9 was an exception, since for very small s-values (here s 1 = 0.086), sampling from the negative binomial distribution results in many zero counts, which then also leads to a clear characteristic of the distribution. Comparing the results between the distinct sample size combinations, the DE analysis on the simulated data sets with the sample size combinations 1 and 2 performed comparably well, whereas with sample size combination 3, which contained the lowest numbers of cells, performed slightly worse. For the case of similar distributions between the groups (II), the false positive rate (FPR) was evaluated with the percentage of genes with p<0.05, thus falsely detecting a difference, whereby no difference exists. For all implemented combinations both methods performed comparably well, whereas the proposed method showed slightly higher FPR, on average 1.47% higher (mean of edgeR: 3.77%, mean of proposed method: 5.24%).
To confirm the validity of this method, we performed a simulation study (see detailed methodological description below) in which we simulated gene expressions with similar mean expression values between COPD and control donors, but categorized gene expressions into two groups: 1) with different distributions between COPD and controls; or 2) with equal distributions between COPD and controls ( Figure S4B). In comparison to the widely used DE method edgeR 23 , the rate of false-positively identified DE genes was comparable to our DE method ( Figure S4C), however, our method showed a significantly lower rate of false-negative results ( Figure S4D). This means that our proposed DE method is able to detect differences in distributions across patient groups even if the differences in mean values are small and mostly the shape of the distributions changed.
Application of the novel DE analysis approach and GSEA DE analysis was performed for all macrophage states and the results are provided in Table S5. For the classification of genes being significantly DE, a test statistic cutoff of 0.75 was chosen. Additionally, for each macrophage state, the DE genes were sorted ascendingly according to their p-values and the 300 top ranked genes were chosen. The visualization of which DE genes are found and shared in which macrophage state was performed using the UpSetR package in R.
Gene set enrichment analysis (GSEA) was performed to identify shared common biological functions by groups of DE genes. The web-tool 'g:Profiler' (version e98_eg45_p14_ce5b097 (95)) was used to perform the functional profiling of the DE genes of interest (genes fulfilling the cutoff criteria for DE genes in >2 macrophage states). As multiple-testing correction method, g:Profiler's in-house g:SCS algorithm was chosen, which corrects for multiple tests that are dependent on each other, which holds true for the hierarchically arranged GO terms. The analysis was done using the Gene Ontology (77,78) database, as well as biological pathway databases, like KEGG (72), Reactome (74) and WikiPathways (79).

Use of publicly available bulk data for validation of results
To investigate whether human leukocyte antigen (HLA) genes, which we found downregulated in macrophages in our single cell data, showed the same trend in a second cohort, we used a bulk transcriptome dataset of human macrophages (GSE13896 (17)), comprising samples from 39 non-smokers, 49 smokers, and 12 COPD patients. We filtered the normalized genes for HLA genes and visualized them as a box plot comparing non-smokers vs. smokers, non-smokers vs. COPD, and COPD vs. smokers using a Wilcoxon rank sum test, as provided in the R package stats. Additionally, to show whether the HLA genes shows statistically significant differences in enrichment between non-smokers and smokers, non-smokers and COPD, and COPD and smokers, we performed GSEA for the respective comparisons using the function GSEA with 10,000 permutations and Benjamini and Hochberg to control the false discovery rate from the package clusterProfiler (version 3.16.1 (96)). The normalized enrichment score (NES) was plotted on the x-axis where a negative NES shows an enrichment on the left hand side and a positive NES shows an enrichment on the right hand side of the plot. The significance of the enrichment was color coded using a negative log10 scale where values above 1.3 were considered as significant.

Cell-Cell Communication
Potential cell-cell-interactions were inferred using 'CellPhoneDB' [version 2. 1.1 (19, 97)]. As input, we used the normalized gene expression matrix of control and COPD patients that was filtered separately for cell types, which were defined by the four-step cell type annotation approach and identified in at least three patients of any group (COPD or control) and contained ≥ 10 cells per patient. Genes were filtered for being expressed in ≥ 5% of a respective cell type. To run CellPhoneDB, the following parameters were set: -iterations=1,000 -pvalue=0.1 -result-precision=10.
In order to visualize the cell-cell communication, we filtered for significant interactions (adjusted (Holm) p-value < 0.05) and summarized the interactions per cell type pair. Network visualizations were done with the 'ggraph' package (version 1.0.2) setting the layout to "fr". To visualize single receptorligand pairs, we filtered for group-specific interactions (-log10 (p-value) > 1) and visualized the resulting interactions for control and COPD.
To evaluate the downstream transcriptomic changes caused by cell-cell-interactions, we applied 'NicheNet' [version 0.1.0 (20, 98)]. As the CellPhoneDB analysis revealed a central role of the C1Q and monocyte-like macrophages in the cellular communication in BALF, we focused on these cells for the subsequent analysis. As the model in NicheNet is based on a different collection of databases than CellPhoneDB, we defined potential sender cell-receiver cell interactions independently of CellPhoneDB. As potential ligands, we accepted all genes that were expressed in >5% of any cell type within the COPD group and which matched at least one receptor from the genes expressed in > 5% of the C1Q macrophages or monocyte-like macrophages in the COPD group, respectively. As input genes to infer the ligand activity score from, we defined all DE genes with a median Wilcoxon score < (-0.75) and p-value of the median Wilcox score <0.05 for each state separately. As background genes, we defined all genes that are not DE in monocyte-like macrophages (or C1Q macrophages) and expressed in > 5% of monocyte-like macrophages (or C1Q macrophages). For ligand prioritization, we selected the top 3 genes with the highest AUPR from each of the comparisons resulting in 6 top ligands.
The expression of these ligands for each cell type was visualized in a heat map scaled by each gene. The target genes of all top ligands were visualized in a heat map with their regulatory potential score for each ligand and their mean expression in C1Q macrophages or monocyte-like macrophages for either COPD or control patients (scaled by gene). To further decipher the exact connection between the ligand and the target genes, we visualized the transcriptional network based on which NicheNet associated the target genes with TGFB1 in a network with free topology. This network was subdivided into receptors for TGFB1, transcriptional regulators between TGFB1 and the target genes. The connections were subdivided into signaling (which does not induce a direct transcriptional change) and transcriptional regulation.

Monocyte-to-macrophage trajectory analysis
To generate a joint embedding of BAL and blood samples, the data were jointly pre-processed using 'Scanpy' [version 1.4.3 commit 0075c62 (68)] on AnnData (version 0.6.22.post2 commit 72c2bde). In concordance with previous analysis, cells from BALF were filtered out if the fraction of mitochondrial reads exceeded 0.1, and a threshold of 0.05 was used for blood samples. Genes that were expressed in fewer than 200 cells were also filtered out. Following previously published best-practices (99) we used scran normalization via the computeSumFactors function on the joint object. Spliced and unspliced counts were mapped to this object using scVelo [version 0.1.24 commit e45a65a (26)]. Quality control for spliced and unspliced counts was performed by removing cells with fewer than 20 spliced and/or 10 unspliced counts. Subsequent normalization by total counts and log-transformation was performed via the filter_and_normalize function from scVelo. Subsetting only relevant monocyte and macrophage populations from blood and BAL datasets (according to the coarse mapping shown in Figure S6A) resulted in a dataset of 57,280 cells and 11,530 genes.
The joint embedding of BAL and blood cells was generated by taking the top 4000 highly variable genes (HVGs) that were shared by most batches. This was done using the hvg_batch function from the single-cell data integration benchmarking package scIB [https://www.github.com/theislab/scib (100)]. This function computes the top 4000 HVGs per batch (here: donor) using Scanpy's highly_variable_genes function with method cell_ranger. These genes are ranked by the number of batches in which each gene is highly variable, and by their mean index of dispersion across all batches. Using this ranked list, we selected the top 4000 genes as a representation of HVGs that are shared across batches. This constitutes a weak integration across batches without direct alteration of the transcriptome data.
Due to an observed batch effect when performing RNA velocity analysis across patients, we ran scVelo per patient and aggregated the individual patient velocities to create a joint velocity embedding. For each donor spliced and unspliced counts were smoothed using the moments function, velocity genes were selected by a stringent log likelihood threshold of 0.1 (between 45 and 172 genes per donor), and the dynamical scVelo model was fit. The resulting inferred single-cell velocities were projected onto the joint UMAP computed from all donors by running velocity_graph on the concatenated object.
Furthermore, partition-based graph abstraction [PAGA (27)] was used to assess the connectivity of cell identity clusters that were suggested to show transitions by RNA velocity. To robustly assess the connectivity of cell identity clusters across donors, we performed PAGA analysis per donor. We computed a kNN graph with Scanpy's neighbors function (k=15) per donor using the joint PCA embedding across donors and ran the paga function on this graph. We used the resulting PAGA connectivities as a statistical test of kNN-graph connectivity between clusters. The median of PAGA connectivities over all donors with both blood and BAL samples was used as a PAGA distance metric.

Data visualization
In general, Seurat and the ggplot2 package [version 3.1.0 (101)] was used to generate figures. For the monocyte-tomacrophage analysis Scanpy, UMAP and scVelo packages were used to generate figures. The graphical summary was created with BioRender.com.

Quantification and statistical analysis
If not otherwise stated, the statistical evaluation was carried out in relation to the total sample size n. A t-test (two-sided) was used for n ≤ 10, otherwise a Wilcoxon rank-sum test was used.

Code availability
We deposited the code for the novel DE analysis approach used in this study on Zenodo (https://doi.org/10.5281/zenodo. 3717776). The analysis code used to generate the majority of the figures are available via FASTGenomics (https://beta. fastgenomics.org/p/bassler_scCOPD).

Data availability statement
The original contributions presented in the study are publicly available. This data can be found here: https://egaarchive.org/studies/EGAS00001004369.

Ethics statement
The studies involving human participants were reviewed and approved by Ethics committees of the University of Bonn and University hospital Bonn. The patients/participants provided their written informed consent to participate in this study.
All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

SUPPLEMENTARY FIGURE 1
Characterization of BALF immune cells using appropriate scRNA-seq technology and four-step cell-type annotation strategy (related to ). (A) UMAP representation of integrated blood and BALF data from different patients and two scRNA-seq technologies (10x Chromium and Seq-Well). The UMAP is split by technology and colored according to identified clusters. Clusters were additionally assigned to cell-types based on marker gene expression (according to Table S6). (B) Stacked bar plots of the relative cell-type proportions for MCFC, which served as ground truth, and cell-type proportions based on the assigned clusters of the two scRNA-seq technologies. (C) Overview of the cell types contained in the reference files used for cell-type annotation. The orange color indicates that the respective cell type is included in the reference file. (D) Confusion plots showing the concordance between the respective cell-type annotations across different annotation methods (SingleR with the Blueprint + ENCODE or HPCA as reference and GenSigPro with signatures from the LM22 dataset as reference). Only cell types that can be found in all reference files as shown in Figure S1C are displayed. (E) Scheme of the gradient boosted decision tree-based machine learningapproach for cell-type annotation. (F) UMAP representation of a benchmarking blood immune cell dataset (according to Table S1). The cells in the UMAPs are colored according to the respective cell annotation methods. The ground truth is derived based on the unique cell-type marker gene expression of each cell. Accumulation of cells that are annotated by the respective annotation methods, but show a deviation in the annotation with respect to the ground truth, are marked with an arrow. (G) UMAP representation of integrated BALF data obtained from COPD patients and control donors via the Seq-Well technology. Coloring according to machine learning-based cell-type annotation. (H) UMAP representation of the integrated dataset with the coloring of DCs, as predicted by the machine learning-based cell-type annotation. Non-DCs are colored light gray. MCFC = multicolor flow cytometry; BALF = bronchoalveolar lavage fluid; DC = dendritic cell; MDM = monocytederived macrophage; SR1 = SingleR (HPCA); SR2 = SingleR (Blueprint + ENCODE); macro = macrophage; mono = monocyte.

SUPPLEMENTARY FIGURE 2
Characterization of identified macrophage states (related to ).  Figure  S2C). Terms of cellular functions found in the same clusters are grouped into bins and the size of the bins is represented as a bar plot on the right, with bins containing more than 25 terms (dashed line) colored red. On the left side, dots indicate which clusters contain and share the binned terms. Frequently occurring terms of cellular functions within the bins containing more than 25 terms are shown. (E) Violin plots (with marked median enrichment values) displaying enrichment of different gene sets across clusters based on the Area Under the Curve (AUC). BALF = bronchoalveolar lavage fluid; sign. = signature; MФ = macrophage; GSVA = gene set variation analysis; mito. = mitochondrial; degrad. = degradation; mod. = modification; present. = presentation.

SUPPLEMENTARY FIGURE 3
Characterization of altered lipid metabolism in macrophages of COPD patients (related to ). (A) Schematic workflow of the GO-shuffling approach. (B) Heat map of lipid metabolism-associated genes predicted by the GO-shuffling approach. The mean gene expression per donor is represented as a z-transformed value (across all donors). Columns and rows of the heat map are sorted by hierarchical clustering. Genes that have been described as causing cholesteryl accumulation through the dysfunctionality of their protein products are marked in yellow. (C) Schema of the key steps in cholesterol metabolism and storage. Metabolites predicted by Compass are highlighted with a gray background. Enzymes involved in metabolism are abbreviated with a number that identifies them also in Figure S3B. (D) Heat map of NOTCH-signaling associated genes predicted by the GO-shuffling approach. The mean gene expression per donor is represented as a ztransformed value (across all donors). Columns and rows of the heat map are sorted by hierarchical clustering. (E) Schema of the time-dependent course of the oxygen consumption rate (OCR) and the inferred mitochondrial parameters based on the injection of different compounds (shown at the top of the plot). dist = distance; visual = visualization.

SUPPLEMENTARY FIGURE 4
Benchmarking of the novel DE-analysis approach (related to ). (A) Schematic workflow of the permutation test-based DE analysis approach. (B) Workflow for the simulation study used for the evaluation of the performance for the proposed DE-method. Two cases are considered, single cell count data was simulated for multiple patients within two groups, with I) different distributions between the groups, II) similar distributions between the groups. For each patient, the mean of the read counts is sampled from the same log-normal distribution with m = 1 and s = 0.15. The read counts are sampled from a negative binomial distribution with success probability (s) parameters dependent on the considered case. When simulating different distributions between the groups (I), distinct s-values are assigned to each group. When simulating the same distributions for the groups (II), the same s-values are assigned to the groups. The number of cells per patient were chosen according to the Seq-Well dataset, where three sample size combinations were considered, originating from three macrophage state datasets. (C) False positive rate (FPR) in simulation study (II) (percentage of genes with a p-value < 0.05) for the proposed method and edgeR for 5 s combinations and for each sample size realization. Boxplots (with marked median values) comprise results of three repeated sets of simulated data (one set of simulated data: 50* m per patient, per s combination). Low FPR denotes a higher number of correct equally expressed gene classifications. (D) False negative rate (FNR) in simulation study (I) (percentage of genes with a p-value > 0.05) for the proposed method and edgeR for 9 scombinations and for each sample size realization. Boxplots (with marked median values) comprise results of three repeated sets of simulated data (one set of simulated data: 50* m per patient, per s combination). Low FNR denotes a higher number of correct differential expressed gene classifications. Investigation of cell-to-cell interactions to infer important signaling pathways in macrophages (related to ). (A) Dot plot representation of monocyte-like macrophage-dependent ligand-receptor interactions predicted by CellPhoneDB that show significant enrichment (represented by the p-value) of the interacting pair in the interacting cell types either in COPD or in the control. Depicted are only selected interactions. (B, C) Illustration of the selection of potential upstream ligands of monocyte-like macrophages or C1Q+ macrophages based on the NicheNet analysis. The histograms show distributions based on ligand activity derived from the area under the precision recall curve (AUPR, upper histogram) and the Pearson correlation coefficient (PCC, lower histogram). The ligand activity of the highest ranked ligands is displayed in a color code together with the names of the 20 highest ranked ligands. The ligands predicted by the CellPhoneDB analysis (according to Figure S5A) are highlighted in red and the top 3 ligands based on AUPR for either monocyte-like macrophages or C1Q+ macrophages (as presented in ) are underlined. (D) Expression of ligand targets from in macrophage subsets comparing COPD and control patients (z-transformed by gene). AUPR = area under the precision recall curve; PCC = Pearson correlation coefficient; MФ = macrophage; DC = dendritic cell.

SUPPLEMENTARY FIGURE 6
Modeling the association of blood monocytes and BALF macrophages (related to ). (A) UMAP of embedded macrophages/monocytes from BALF and blood monocytes with coloring according to the cell types derived from the combined labels. The dendrogram on the right side illustrates the transcriptional relationship between the macrophage subtypes and shows how several subtypes were summarized in the combined labels. (B) Projection of computed RNA velocity vectors onto the UMAP of the embedded data. BALF = bronchoalveolar lavage fluid; mono = monocyte; MФ = macrophage; proliferat. = proliferating.

SUPPLEMENTARY FIGURE 7
Instructions for accessing the COPD Seq-Well dataset and scripts via the FASTGenomics platform.