Macrophages with reduced expressions of classical M1 and M2 surface markers in human bronchoalveolar lavage fluid exhibit pro-inflammatory gene signatures

The classical M1/M2 polarity of macrophages may not be applicable to inflammatory lung diseases including chronic obstructive pulmonary disease (COPD) due to the complex microenvironment in lungs and the plasticity of macrophages. We examined macrophage sub-phenotypes in bronchoalveolar lavage (BAL) fluid in 25 participants with CD40 (a M1 marker) and CD163 (a M2 marker). Of these, we performed RNA-sequencing on each subtype in 10 patients using the Illumina NextSeq 500. Approximately 25% of the macrophages did not harbor classical M1 or M2 surface markers (double negative, DN), and these cells were significantly enriched in COPD patients compared with non-COPD patients (46.7% vs. 14.5%, p < 0.001). 1886 genes were differentially expressed in the DN subtype compared with all other subtypes at a 10% false discovery rate. The 602 up-regulated genes included 15 mitochondrial genes and were enriched in 86 gene ontology (GO) biological processes including inflammatory responses. Modules associated with cellular functions including oxidative phosphorylation were significantly down-regulated in the DN subtype. Macrophages in the human BAL fluid, which were negative for both M1/M2 surface markers, harbored a gene signature that was pro-inflammatory and suggested dysfunction in cellular homeostasis. These macrophages may contribute to the pathogenesis and manifestations of inflammatory lung diseases such as COPD.

RNA-sequencing. Next, we performed RNA-sequencing on BAL macrophages collected from 10 consecutive participants between March 2019 and June 2019. After exclusion of 3 samples (i.e. 1 DN and 2 DP's, which were obtained from three patients with asthma) because of poor RNA quality or low RNA yield, we subjected the remaining samples (n = 37) to bulk-RNA sequencing. The demographic data of the study participants are summarized in Table S1.
Differentially expressed genes (DEGs) between subtypes. We compared the transcriptomic expression pattern across all four subtypes of macrophages. The greatest number of differentially expressed genes was observed with the DN subtype (1886 differentially expressed genes versus all other subtypes at a 10% false discovery rate, FDR). There were 498 differentially expressed genes between the DP subtype and the others; 15 genes between the M1 subtype and the others; and 52 genes between the M2 subtype and the others (Fig. 1). All differentially expressed genes at 10% FDR are shown in Table S2.
Among the differentially expressed genes, 602 were up-regulated in the DN subtype, 292 were up-regulated in DP, 12 were up-regulated in M1 and 11 were up-regulated genes in the M2 subtype. The top 20 up-regulated genes for each subtype are shown in Fig. 1. Importantly, the top 20 up-regulated genes for the DN macrophages included 15 mitochondrial genes and 4 mitochondrial pseudogenes; for other macrophage subtypes, these mitochondrial genes were not differentially expressed.
Enrichment analysis was performed using only the up-regulated genes for each macrophage subtype at 10% FDR. Up-regulated genes in DN were enriched in 86 GO biological processes. The results of the enrichment analysis including the top 5 strongest p-values are shown in Fig. 2. These included genes involved in inflammatory responses (FDR = 3.89E−04).
Likewise, we compared the transcriptomic expression pattern on the basis of CD163 positivity. Among 128 differentially expressed genes at 10% FDR, 39 genes were up-regulated and 89 genes were down-regulated in macrophages positive for CD163 (Fig. S1). All differentially expressed genes at 10% FDR are shown in Table S4. However, none of the up-or down-regulated genes were enriched in GO biological processes at 10% FDR.
Weighted gene co-expression network analysis (WGCNA). Based on detectable expression levels of all genes on the RNAseq platform (18,314 genes in total), 13 gene expression modules were constructed using a weighted gene co-expression network analysis (WGCNA). The WCGNA-derived modules ranged in size from 19 genes in module 13 to 1,831 genes in module 1. The 13 modules as well as the "garbage module" (i.e. module 00) derived from WGCNA are shown in Table S5.
Among these module eigengenes, nine of them were differentially expressed in at least one macrophage subtype at a 10% FDR threshold (Fig. 3). Gene signatures in the DN subtype were significantly distinct from those of other subtypes. In contrast, M1 macrophages were not associated with a distinct module. Importantly, module 12, which was up-regulated in DNs, was down-regulated in the DP and M2 subtypes. This module was composed of 15 mitochondrial genes, 13 protein subunit genes, 2 ribosomal RNA and 4 mitochondrial pseudogenes. The top genes with the highest membership in this module were MT-ATP6, MT-CYB and MT-ND4. www.nature.com/scientificreports/ Enrichment analysis was performed at 10% FDR to identify GO biological processes that were associated with these modules. The top three GO processes are shown in Fig. 4. Module 2 was associated with the function of ATP production in mitochondria including oxidative phosphorylation (FDR = 6.32E−44) and respiratory electron transport chain (FDR = 2.62E−31). Module 3 was associated with fatty oxidation and metabolism (FDR = 1.49E−03). Module 4 was associated with RNA splicing (FDR = 4.88E−08) and mitotic nuclear division (FDR = 4.88E−08). Module 6 was associated with regulation of protein localization to telomeres (FDR = 2.81E−04). Mitochondrial genes in module 12 were not included in the enrichment analysis because of the absence of a reference gene set for these genes. All enriched GO biological processes at 10% FDR and the top 20 genes with the highest membership for each of the modules are shown in Tables S6 and S7, respectively. To further explore potential changes in mitochondrial respiration in the DN subtype, we examined the relative gene expression in the GO biological process of oxidative phosphorylation, which overlapped with those in www.nature.com/scientificreports/ module 2, and the mitochondrial genes in module 12. The oxidative phosphorylation genes were down-regulated in the DN subtype ( Fig. 5a), whereas the mitochondrial genes were up-regulated in the DN subtype (Fig. 5b).
Macrophage sub-phenotype distributions across diseases. Among 25 participants including 8 with COPD and 13 with asthma, we investigated the distribution of macrophage subtypes based on disease. Because the presence of asthma did not significantly alter the macrophage subtype distribution (Fig. S2), we compared the subtypes between participants with and without COPD. Patients with COPD were more likely to be males (p = 0.003) and current smokers (p = 0.001) and had lower FEV 1 /FVC (p = 0.006). None of the patients with COPD had been previously diagnosed with asthma. A majority of patients with COPD had moderate to severe airflow limitation. Baseline characteristics are shown in Table 1. Macrophages in patients with COPD were less likely to express classical surface markers including CD 163 compared to those without COPD (39.3% vs. 71.6%, p = 0.004). For the four macrophage sub-phenotypes, the DN subtype was enriched in patients with COPD (46.7% vs. 14.5%, p < 0.001). In contrast, the DP subtype was significantly reduced in COPD (16.5% vs. 44.3%, p = 0.001). There were no significant differences in terms of M1 and M2 subtypes (Table 1).

Discussion
Here, we showed that many macrophages in human BAL did not conform to the M1/M2 paradigm, and that these cells, especially those that did not harbor M1 or M2 cell surface markers, demonstrated distinct transcriptomic signatures. Interestingly, the double negative subtype contained differential expression of mitochondrial genes, which were significantly enriched in patients with COPD. Although previous studies have attempted to characterize the gene signatures in the framework of dichotomized classification in vitro, or in vivo, to the best The circle size represents the number of overlapping genes between each GO process and up-regulated genes according to the subtype. The colour scale represents the extent to which the up-regulated genes are significantly enriched in each GO process. The GO processes associated with DN, DP and M1 subtype included inflammatory response, complement activation and response to virus, respectively, whereas none of the up-regulated genes for M2 macrophages were enriched in the GO processes. Figure created with the R package ggplot2. https:// cran.r-proje ct. org/ web/ packa ges/ ggplo t2/ index. html.
To date, the prevalence of M1 and M2 macrophages in vivo has not been well characterized. A priori we decided to use CD40 and CD163 based on previous studies, which showed that they were distinct markers for human M1 and M2 macrophages [9][10][11] . A number of studies have shown that CD40 and CD163 are expressed on 4-20% and 45-60%, respectively, of human lung macrophages in non-smoking individuals 9,15,16 . We extend these findings by showing that approximately 25% of all macrophages in human BAL fluid could not be phenotyped with CD40 and CD163 (and thus were deemed "double negative" macrophages); and that the percentage of these double negative cells increased significantly in BAL fluid of patients with COPD. This observation is consistent with previous studies which reported reduced expression of CD163 and CD40 on alveolar macrophages in COPD lungs 9,17 .
In WGCNA, interestingly, we found that module 2, which was associated with oxidative phosphorylation in mitochondria, was down-regulated in double negative macrophages. Mitochondria is the main producer of cellular energy by means of oxidative phosphorylation, which involves electron-transferring respiratory chain (complexes I-IV) and adenosine triphosphate (ATP) synthase (complex V). Although mitochondria has its own DNA, which contains 37 genes, over 98% of the mitochondrial proteins are encoded by the nuclear genome 25,26 . We The origins of double negative macrophages are unknown. Previous studies have shown that aged cells have reduced expression of surface markers including MHC class II molecules and co-stimulatory receptors such as CD40, which raises the possibility that the double negative macrophages may be senescent cells 28,29 . This notion is further supported by our data showing that these cells harbor differentially expressed genes involved in inflammation and mitochondrial (dys)function 30 . In the present study, we could not distinguish tissue resident macrophages from monocyte-derived macrophages 31 . It is noteworthy, however, MARCO was down-regulated in double negative macrophages (logFC = − 0.35, FDR = 0.021) compared to the other subtypes. MARCO is a marker of embryonically-derived resident macrophages, which raises the possibility that many of the double negatives originated from monocyte-derived macrophages 32 . Specific studies focused on cell lineages will be needed to pinpoint the exact source of these and other macrophage subtypes.
There were several limitations in this study. First, some of the up-regulated genes in M1 or M2 macrophages in our study did not fully align with classical patterns associated with M1-or M2-related genes. For instance, up-regulated genes in M1 macrophages included those encoding for proteins such as CCL22 and MMP12, which have been related to M2 macrophages 33,34 . Also, genes encoding CD40 and CD163 did not co-express with genes encoding other surface markers commonly used to characterize M1 or M2 macrophages such as CD80 or CD206, respectively 5 . However, these markers were identified using in vitro generated macrophages where M1-related markers were induced by LPS/IFN-γ and M2-related markers were induced by IL-4/IL-13. A comparative study of in vivo and in vitro macrophages showed that the gene signature of M1 macrophages activated by LPS in vivo shared many features of alternatively activated M2 genes in vitro including up-regulation of CCL22 35 . While not all M1 and M2 markers in vitro can be directly translated to the in vivo situation 35 , further studies are needed to validate our findings with other M1 and M2 markers using technologies such as flow cytometry and single cell sequencing. Second, although the unique gene expression signature of double negative macrophages was interesting, we did not conduct functional studies to validate the potential mitochondrial dysfunction in these macrophages, which were suggested by the RNA sequencing data. Third, it is possible that some of our samples may have contained dendritic cells as they have similar morphology and share common surface markers such as HLA-DR and CD40 as macrophages 36 . Although dendritic cells generally constitute only 0.5% of total cells in BAL fluid, future studies should gate out these cells for more precise assessment of alveolar macrophages 37,38 .
In conclusion, approximately one out of four macrophages in human BAL fluid stain negatively for both M1/M2 cell surface markers. These double negative macrophages harbor gene expression signature that is Bronchoscopy procedure. All procedures were performed by an experienced pulmonologist. Conscious sedation was first provided to the patient with the use of intravenous midazolam and fentanyl. Through the working channel of the bronchoscope (Olympus Corporation, Tokyo, Japan), topical 2% lidocaine was instilled, if needed, to prevent bronchospasm and cough. BAL was generally taken from the right middle lobe or lingula unless these lobes had disease noted on a pre-operative computed tomography (CT) imaging (e.g. pulmonary www.nature.com/scientificreports/ nodule). For these cases, the right or left upper lobe was used. After the bronchoscope was fully wedged into the desired segment, 20 ml of 0.9% saline was instilled (with a dwell time of 10 s) and then the fluid was manually aspirated out using a vacuum syringe. The first aliquot of the recovered solution was discarded because of concerns over contamination by bronchial lining fluid or mucus secretions. Aliquots of 40 ml of saline were then sequentially instilled to a maximum volume of 200 mL or until 30-50 mL of BAL fluid was recovered, whichever came first.
Sample preparation. The recovered BAL fluid was put in ice immediately after aspiration and filtered through sterile 70 µm DNase/RNase-free cell strainers to remove large clumps and debris. Cells were recovered by centrifugation at 500 g for 10 min at 4 °C and were washed twice with PBS. The concentration of viable cells was determined by trypan blue staining on a hemocytometer. These samples generally contained 5 × 10 5 to 5 × 10 6   Statistical analysis and RNA sequencing data processing. All statistical analyses were performed in R. For patient characteristics, continuous data including the distribution of macrophage subtypes were represented as mean ± standard deviation (SD) and categorical data as numbers (%) of observations. Statistical significance between COPD and non-COPD patients was assessed using a student's t-test for continuous variables, and a Fisher's exact test for categorical variables. P values < 0.05 were considered significant.
In RNA-seq data processing, raw sequencing reads were quality controlled using FastQC 41 . STAR (Spliced Transcripts Alignment to a Reference) was used to align the reads to GENCODE GRCh37 (version 31) genome reference and RSEM (RNA-Seq by Expectation Maximization) was used for quantification to obtain the counts and the transcript per million (TPM) 42,43 . The principal component analysis was used to check for potential batch effect and confounding factors. No obvious batch effect was observed but smoking status showed a potential confounding effect on the gene expression data (Fig. S4). Limma voom was used to normalize the count to log 2 counts per million (CPM) 44 . Log 2 CPM was used for the differential expression analysis while TPM, which normalizes for gene-length was used for the weighted gene co-expression network analysis (WGCNA). Genes with low abundance (log 2 CPM < 1 or TPM < 5 in more than one-fourth of the samples) were filtered out.
For differential expression analysis, limma's mixed effect model with adjustment of smoking status was used to compare one macrophage subtype versus the rest of the subtypes to identify characteristic gene expression signatures for each macrophage subtype. The Benjamini-Hochberg procedure was used to correct for multiple hypothesis testing and to control the false discovery rate (FDR); 10% FDR was used in line with previous studies which also have used 10% FDR [45][46][47] . A Gene Ontology (GO) enrichment analysis was performed on the significantly up-regulated genes at 10% FDR for each macrophage subtype using the R package, clusterProfiler 48,49 .
For the weighted gene co-expression network analysis (WGCNA), the R package WGCNA was used to construct gene modules with genes that were co-expressed with each other. The R code for WGCNA is available at https:// github. com/ yyola nda/ macro phage_ rnaseq. Given that gene length may have an effect on gene clustering, the gene-length normalized TPM was used for WGCNA 48 . We chose a soft-thresholding power (β) of 29 to obtain a scale free topology model fit index (R 2 ) of > 0.8 (Fig. S5). The minimum module size was set to 15 genes and modules with a distance < 0.2 were merged together. A signed network with thirteen modules, excluding the "garbage" module which contained genes that did not co-express with any other genes, was constructed and each module was assigned a number. The "garbage" module was discarded and not included in the downstream analysis. The eigengene of each module was obtained by calculating the first principal component of the module and was considered as a representative of the expression profile of the module. Figure S6 shows a dendrogram of the clustering of the module eigengenes and a heatmap of the correlation between the module eigengenes. Limma's mixed effect model was used to identify the module eigengenes that were associated with each macrophage subtype after adjusting for smoking status. For modules that were significantly associated with any of the macrophage subtypes at 10% FDR, a Gene Ontology (GO) enrichment analysis was performed using the R package clusterProfiler to identify biological processes that were associated with the genes of each module. To further investigate these modules, we ranked the module genes by their module membership. The top genes with larger and positive module membership were highly connected to the other genes in the module and exhibited the same direction of effect as the module eigengene. As a sensitivity analysis, we performed the analyses using a 5% FDR; these results are shown in Tables S8, S9, S10 and S11.

Data availability
The datasets used for the current study are available from the senior author (DDS) upon a reasonable request.