Macrophages represent an important component of the tumor microenvironment and play a complex role in cancer progression. These cells are characterized by a high degree of plasticity, and they alter their phenotype in response to local environmental cues. Whereas the M1/M2 classification of macrophages has been widely used, the complexity of macrophage phenotypes has not been well studied, particularly in lung cancer. In this study we employed an orthotopic immunocompetent model of lung adenocarcinoma in which murine lung cancer cells are directly implanted into the left lobe of syngeneic mice. Using multimarker flow cytometry, we defined and recovered several distinct populations of monocytes/macrophages from tumors at different stages of progression. We used RNA-seq transcriptional profiling to define distinct features of each population and determine how they change during tumor progression. We defined an alveolar resident macrophage population that does not change in number and expresses multiple genes related to lipid metabolism and lipid signaling. We also defined a population of tumor-associated macrophages that increase dramatically with tumor and selectively expresses a panel of chemokine genes. A third population, which resembles tumor-associated monocytes, expresses a large number of genes involved in matrix remodeling. By correlating transcriptional profiles with clinically prognostic genes, we show that specific monocyte/macrophage populations are enriched in genes that predict outcomes in lung adenocarcinoma, implicating these subpopulations as critical determinants of patient survival. Our data underscore the complexity of monocytes/macrophages in the tumor microenvironment, and they suggest that distinct populations play specific roles in tumor progression.
Lung cancer remains the leading cause of cancer death in men and women, and most patients die of complications of metastatic disease (1). Whereas tumor initiation is driven by somatic mutations in oncogenic drivers and tumor suppressor genes, progression and metastasis involve critical crosstalk between cancer cells and the tumor microenvironment (TME) (2). Among the diverse cell types of the TME, macrophages have been implicated as important participants and are the most abundant non–tumor cell type in most cancers (3). Studies in several types of cancer, including lung cancer, have demonstrated that macrophage depletion either through pharmacological or genetic approaches results in slower tumor growth, implicating these cells as mediators of tumor progression (4, 5). In lung cancer, several studies have described correlations in macrophage number or phenotype with clinical outcomes (6). Many of these studies have employed the model of M1 and M2 macrophages. This model, which is based on activation of macrophages in vitro, designates a macrophage phenotype as being either proinflammatory M1, which is proposed to inhibit cancer progression, or alternatively activated M2, which is proangiogenic, immunosuppressive, and promotes cancer progression (7, 8). Although the M1/M2 classification of macrophages serves as a useful starting point, it does not take into account the complexity and plasticity of these cells. Recent studies have demonstrated that distinct populations of macrophages exist in tissues such as the lung, and that these different populations not only express distinct markers but also have different developmental origins (9). For example, the resident alveolar macrophages in the lung are derived from the embryonic yolk sac, whereas the recruited monocytes and macrophages are derived from the bone marrow.
These different populations are likely to play distinct and potentially opposing roles in cancer progression. In fact, tumors are infiltrated by different dynamically changing populations of monocytes/macrophages (8), and identifying key features of each population may help to target them therapeutically.
To study the role of the TME in lung cancer, we have developed an orthotopic immunocompetent mouse model, in which Lewis lung carcinoma (LLC) cells, which were derived from a spontaneous lung adenocarcinoma in C57BL/6 mice, are directly implanted into the lungs of fully immunocompetent syngeneic C57BL/6 mice. Within 2–3 wk these cells form a primary tumor that metastasizes to the other lung lobes, mediastinal lymph nodes, liver, and brain. We have previously shown that tumor progression in this model is associated with increases in several monocyte/macrophage populations (10). However, functional differences between these populations have not been identified, and the role of specific monocyte/macrophage subtypes in tumor progression is unknown. In this study, we have used multimarker flow cytometry to identify and isolate resident populations as well as dynamically expanding populations of monocytes and macrophages. Using RNA sequencing (RNA-seq) we have been able to profile differences in gene expression between these macrophage populations and determine how these profiles change during tumor progression. Furthermore, using transcriptional profiles, we were able to correlate specific monocyte/macrophage populations to clinical outcomes in lung adenocarcinoma.
Materials and Methods
LLC cells stably expressing luciferase were obtained from Caliper Life Sciences and maintained in DMEM (Corning CellGro, no. 10-017-CV) containing 10% FBS, penicillin/streptomycin, and G418 (500 ng/ml).
Orthotopic mouse model
C57BL/6 mice were obtained from The Jackson Laboratory and bred and maintained in the Center for Comparative Medicine at the University of Colorado Denver. Experiments were performed in 10- to 16-wk-old males. LLC cells stably expressing luciferase (1 × 105 in 25 μl per injection) were suspended in PBS containing 15% growth factor reduced Matrigel (BD Biosciences, catalog no. 354230) and injected into the parenchyma of the left lung lobe through the rib cage using a 30-gauge needle, as previously described (10). During injection, the lung was directly visualized through a 4- to 5-mm incision in the skin. After the procedure, the incision was closed using veterinary adhesive. All procedures were performed under protocols approved by the Institutional Animal Care and Use Committee at the University of Colorado Denver.
Preparation of single-cell suspensions and flow cytometry
Tumor-bearing mice were sacrificed either 2 or 3 wk after cancer cell injection, along with uninjected control mice. The lung circulation was perfused with PBS/heparin (20 U/ml, Sigma-Aldrich). Left lung lobes containing tumors were excised and weighed. For control mice, all the lung lobes were harvested for analysis. Tissues from three to five mice were pooled, mechanically dissociated with scissors, and incubated for 30 min at 37°C with 3.2 mg/ml collagenase type 2 (Worthington Biochemical, 43C14117B), 0.75 mg/ml elastase (Worthington Biochemical, 33S14652), 0.2 mg/ml soybean trypsin inhibitor (Worthington Biochemical, S9B11099N), and 40 μg/ml DNAse I (Sigma-Aldrich). During incubation, samples were placed in a shaking water bath and dispersed by pipetting every 10 min. Resulting single-cell suspensions were filtered through 70-μm cell strainers (BD Biosciences) and washed with staining buffer (PBS containing 1% FBS, 2 mM EDTA, 10 mM HEPES). Samples were subjected to RBC lysis, washed with staining buffer, and filtered through 40-μm cell strainers (BD Biosciences). Prior to staining, FcγR was blocked with anti-CD16/CD32 Ab (BD Biosciences) for 10 min. Cells were stained for 30–45 min at 4°C with the following Abs: CD11b-FITC (clone M1/70, BD Biosciences), Siglec-F–PE or Alexa Fluor 647 (clone E50-2440, BD Biosciences), Ly6G-PE-Cy7 (clone 1A8, BD Biosciences), CD64-PE or Alexa Fluor 647 (clone X54-5/7.1, BD Biosciences), CD11c-allophycocyanin-Cy7 (clone HL3, BD Biosciences), and CD11c-PerCP-Cy5.5 (clone N418, BioLegend). Flow cytometry analysis and cell sorting were conducted at the University of Colorado Cancer Center Flow Cytometry Core Facility using a Gallios flow cytometer (Beckman Coulter) for analysis and XDP-100 or Astrios (Beckman Coulter) cell sorters for sorting. The sorting strategy involved exclusion of debris and cell doublets by light scatter and dead cells by DAPI (1 μg/ml). Absolute cell counts were obtained using Sphero AccuCount ultra rainbow fluorescent particles 3.8 (Spherotech) per the manufacturer’s instructions. Data were analyzed using Kaluza software (Beckman Coulter).
RNA extraction and RNA-seq
Total RNA was isolated from flow cytometry–sorted cells using QIAshredders and an RNeasy Plus micro kit (Qiagen). RNA quality and quantity were analyzed using a NanoDrop and Bioanalyzer. RNA-seq library preparation and sequencing were conducted at the Genomics and Microarray Core at the University of Colorado Denver–Anschutz Medical Campus. Libraries were constructed using a NuGen Ovation human FFPE RNA-seq multiplex system kit customized with mouse specific oligonucleotides for rRNA removal. Directional mRNA-seq was conducted using the Illumina HiSeq 2000 system, using the single-read 100 cycles option.
RNA-seq reads were obtained using Illumina HiSeq analysis pipeline. Reads quality was checked using FastQC (http://www.bioinformatics.bbsrc.ac.uk/projects/fastqc). The median number of reads per condition was 24 million. Reads were then processed and aligned to the University of California Santa Cruz Mus musculus reference genome (build mm10) using TopHat v2 (11). TopHat incorporates the Bowtie v2 algorithm to perform the alignment (12). TopHat initially removes a portion of reads based on quality information accompanying each read and then maps reads to the reference genome. The prebuilt M. musculus University of California Santa Cruz mm10 index was downloaded from the TopHat homepage and used as the reference genome. The default parameters for TopHat were used. The aligned read files were processed by Cufflinks v2.0.2 (13). Reads were assembled into transcripts and their abundance was estimated. Cufflinks uses the normalized RNA-seq fragment counts to measure the relative abundances of transcripts. The unit of measurement is fragments per kilobase of exon per million fragments mapped (FPKM). Confidence intervals for FPKM estimates were calculated using a Bayesian inference method (13). EdgeR is a Bioconductor software package for examining differential expression of replicated count data (14). Briefly, it used an overdispersed Poisson model to account for both biological and technical variability. Additionally, it used empirical Bayes methods to moderate the degree of overdispersion across transcripts, thereby improving the reliability of inference. Default parameters for EdgeR were used.
Hierarchical clustering was performed using Cluster 3.0 (15). Gene expression levels were centered on the mean and normalized. Clustering was done using complete linkage and Euclidean distance as the metric. Heat maps were generated using Java TreeView (16). Pathway overrepresentation analysis was conducted using ConsensusPathDB (http://cpdb.molgen.mpg.de/MCPDB) (17, 18). Searches were conducted in KEGG and Reactome databases with cut-off values of p < 0.01 and a minimum of three overlapping genes. The set of initially selected 2458 differentially expressed genes was used as background. The analyses were additionally confirmed by DAVID (http://david.abcc.ncifcrf.gov/) (19, 20) and Gene Ontology (21, 22).
For comparison with profiles of immune cells published by the Immunological Genome Project (ImmGen), the following populations from naive mice were used: lung alveolar macrophages (MF_Lu), lung CD11b+ macrophages (MF_11b+_Lu), blood monocytes classical (Mo_6C+II+_Bl, Mo_6C+II-_Bl), and blood monocytes, nonclassical (Mo_6C-II-_Bl, Mo_6C-II+_Bl, Mo_6C-IIint_Bl). The top 100 genes differentially expressed in any pairwise comparison between these populations were downloaded from the Web site (http://www.Immgen.org) and combined into a list of 645 genes. The expression of these 645 genes among the populations in this study was examined and subjected to hierarchical clustering. Clusters of genes highly expressed selectively in one cell type were identified. Expression of genes in these clusters was mapped back to the reference ImmGen populations using the ImmGen Web site tool My Gene Set.
For the analysis of clinically prognostic genes, the list of ∼23,000 genes and their associated survival meta-Z scores for lung adenocarcinoma were downloaded from https://precog.stanford.edu/index.php. The survival scores were used as the metric to rank the genes, and the gene set enrichment analysis (GSEA; preranked tool) (23, 24) was used to examine enrichment of good or bad prognostic genes in the gene clusters highly expressed in MacA or MacB cells as defined in Fig. 3.
The RNA-seq data reported in this study have been deposited in the National Center for Biotechnology Information Gene Expression Omnibus repository (http://www.ncbi.nlm.nih.gov/geo/), with the accession no. GSE76033.
Quantitative real-time PCR
Total RNA was isolated from flow cytometry–sorted cells using QIAshredders and an RNeasy micro kit (Qiagen), and converted into cDNA using an iScript cDNA synthesis kit (Bio-Rad). Quantitative real-time PCR was conducted as previously described (10).
Mice were sacrificed 2.5 wk after cancer cell injection, along with uninjected control mice, and macrophage populations were recovered by flow cytometry. Immediately after sorting, cells were stimulated with a calcium ionophore, and eicosanoids were extracted and measured by liquid chromatography/tandem mass spectrometry previously described (10). Results were normalized to total protein.
Identification of myeloid cell subsets in the lung during tumor development
LLC cells were injected directly into the left lung lobe of fully immunocompetent syngeneic wild-type C57BL/6 mice. We harvested the whole tumor-bearing left lung lobes 2 and 3 wk after injection, along with lungs from uninjected control mice. Analysis of weights revealed rapid tumor growth within the lung parenchyma (Fig. 1A). Tissues were processed into single-cell suspensions and analyzed by flow cytometry using myeloid-specific cell-surface markers.
In addition to the CD11b and CD11c markers, which are traditionally used for identifying macrophages in the lung (25, 26), our strategy employed Siglec-F and CD64 markers (9). Siglec-F is a marker specific for murine lung-resident alveolar macrophages that is not expressed by interstitial or inflammatory macrophages. Although Siglec-F also stains eosinophils, their lack of CD11c expression allows us to exclude eosinophils from the monocyte/macrophage pool (9, 27–29). CD64 is a macrophage-specific marker, and in contrast to F4/80, its expression is low or negative on monocytes, dendritic cells, and granulocytes (9).
Using a modification of our previously published gating strategy (10) shown in Fig. 1B, we identified multiple subsets of myeloid cells. In both naive and tumor-bearing lungs, we identified MacA cells, which bear markers of alveolar macrophages (Siglec-F+/CD11c+). Furthermore, we excluded eosinophils (Siglec-F+/CD11c−, Eos gate) and neutrophils (Ly6G+/CD11b+, Neu gate). The remaining CD11b+ cells (MacB gate) were separated into three distinct populations based on CD11c and CD64 expression: MacB1 cells (CD11b+/CD64lo/CD11c+), MacB2 cells (CD11b+/CD64int/CD11c−), and MacB3 cells (CD11b+/CD64hi/CD11c+). As shown in the flow cytometry plots, the MacB compartment in naive lungs was dominated by MacB1 and MacB2 cells, with MacB3 present at very low numbers. Experiments quantifying absolute cell numbers of these populations in the lung revealed that the total numbers of MacA and MacB1 cells in the left lobe did not change with tumor growth, whereas the MacB2 and MacB3 populations rapidly expanded with increasing tumor size, with the most rapid expansion seen in MacB3 cells (Fig. 1C).
We further characterized MacA and MacB cells using Ly6C and MHC class II (MHC II), which are markers frequently used to phenotype murine myeloid populations by flow cytometry (Fig.1D). MacA cells, in the presence or absence of tumor, were negative for Ly6C and stained weakly for MHC II, as expected for alveolar macrophages. MacB1 cells expressed low levels of Ly6C in the presence or absence of tumor. In naive lung, this gate contained both MHC II− and MHC II+ populations, indicating that the MacB1 gate was composed of at least two distinct cell types, with the MHC II+ cells likely representing CD11b+ dendritic cells (9, 30). In the presence of tumor, the MHC II expression on the cells in the MacB1 gate was heterogeneous, ranging from negative to high. MacB2 cells, both with and without tumor, expressed high levels of Ly6C but had little to no expression of MHC II, which is typical for monocytes (31). MacB3 cells, which in naive mice were present in very low numbers, were negative for Ly6C and expressed high levels of MHC II, which is consistent with interstitial/infiltrating lung macrophages (31). In the presence of tumor, the MacB3 gate was composed of cells with low expression of Ly6C and intermediate expression of MHC II.
Recovery of cells and transcriptional profiling
To gain further insight into the identity and functional properties of myeloid cells in the lung TME, we performed transcriptome profiling using RNA-seq. We used flow cytometry–based sorting to recover cells from tumor-bearing lungs 2 and 3 wk after cancer cell injection, as well as from uninjected control mice. RNA was extracted from three pools of mice that were independently injected and harvested, with three to five mice contributing to each pool. We obtained sufficient cell numbers and RNA quantity to sequence the following populations isolated from three independent experiments: from naive mice, MacA cells (MacA-N), MacB1 cells (MacB1-N), and MacB2 cells (MacB2-N); from 2 wk tumor-bearing mice, MacA cells (MacA-2wk), MacB2 cells (MacB2-2wk), and MacB3 cells (MacB3-3wk); and from 3 wk tumor-bearing mice, MacB2 cells (MacB2-3wk) and MacB3 cells (MacB3-3wk). We performed RNA-seq profiling on these samples, compared their transcriptional profile to published databases, and analyzed differential gene expression between these populations.
Comparison to ImmGen profiles
To verify the identity of MacA and MacB populations, we compared their transcriptional profiles to immune cells from naive mice analyzed by the ImmGen (32). As a reference we focused on nine populations of myeloid cells identified by the ImmGen: lung macrophages (alveolar macrophages and CD11b+ lung macrophages), lung dendritic cells (CD11b+ and CD11b−), and blood monocytes (Ly6C+/MHC II−, Ly6C+/MHC II+, Ly6C−/MHC II−, Ly6C−/MHC II+, Ly6C−/MHC IIint). We selected 645 genes that were differentially expressed between these nine reference populations and analyzed their expression in MacA and MacB cells. Using hierarchical clustering, we delineated five gene clusters that were highly expressed in a specific MacA or MacB population (Fig. 2). To examine whether these clusters mapped back to specific ImmGen cell types, we plotted their expression in the nine ImmGen populations using the My Gene Set tool (http://www.ImmGen.org). This analysis indicated that genes specific to MacA-N and MacA-2wk (cluster I) are highly expressed in ImmGen alveolar lung macrophages. For MacB1 cells, we defined two separate clusters (IV and V). Genes in cluster V map to CD11b+ lung dendritic cells, whereas genes in cluster IV map to Ly6C−/MHC II− and Ly6C−/MHC IIint monocytes. This indicates that MacB1-N is a mixed population composed of lung CD11b+ dendritic cells and Ly6C− monocytes. MacB2 cells (MacB2-N, MacB2-2wk, and MacB2-3wk) express genes in cluster III, and these genes correlate with Ly6C+/MHC II+ and Ly6C+/MHC II− blood monocytes. Lastly, MacB3-2wk and MacB3-3wk express high levels of genes in cluster II, and this cluster maps to CD11b+ lung macrophages. Thus, the genes that are highly expressed in MacA or MacB populations allowed us to map them to specific immune cells based on published transcriptional profiles in naive mice. Furthermore, this analysis is in concordance with our immunophenotyping and Ly6C and MHC II expression by flow cytometry.
Global analysis of differentially expressed genes
The analysis of RNA-seq data allowed us to identify differentially expressed genes between the cell populations and between time points for respective cell populations. In brief, in a pairwise analysis a differentially expressed gene was expected to have expression level FPKM of >5 in at least one condition, and have >2.5-fold change in log2 expression with a false discovery rate–adjusted p value <0.05 between the two conditions. In all pairwise comparisons between cell types and/or time points, we identified 2458 genes that were differentially expressed (Supplemental Table I). The large number of differentially expressed genes suggests that although all the analyzed cell populations belong to the mononuclear phagocyte lineage, significant biological differences exist between them.
To identify relationships between genes and among MacA and MacB populations, we performed unbiased hierarchical clustering of all 2458 differentially expressed genes and visualized the data using a dendrogram and a heat map (Fig. 3). The top dendrogram indicates that the individual replicates of the cell populations cluster together. Interestingly, the replicates of MacB2-N and MacB2-2wk were closely clustered. This may reflect the fact that in 2 wk samples, the harvested tissue was composed of smaller tumors with a larger contribution from the uninvolved lung. Otherwise, genes could be grouped based on their expression into five broad clusters: upregulation (or downregulation) of gene expression in MacA-N and MacA-2wk (cluster A); MacB1-N (cluster B1); MacB2-N and MacB2-2wk (cluster B2); MacB2-3wk (cluster B2-3wk); and MacB3-2wk and MacB3-3wk (cluster B3).
We conducted pathway overrepresentation analysis of the gene clusters using Web-based bioinformatic tools at ConsensusPathDB (http://consensuspathdb.org/) (17, 18), searching in the KEGG and Reactome databases. The analyses were additionally confirmed by DAVID (http://david.abcc.ncifcrf.gov/) (19, 20) and Gene Ontology (21, 22) analyses independently. Below, we discuss characteristic pathway signatures of individual clusters in greater detail.
We first focused on genes highly expressed by MacA-N and MacA-2wk (see Supplemental Table I). Although MacA cells do not increase in number during tumor growth in our model, they represent pre-existing resident alveolar macrophages that are unique to the lung and may play a critical role in early stages of lung cancer or in the setting of slowly developing tumors. In cluster A, pathway overrepresentation analysis highly ranked multiple gene sets related to lipid metabolism (Fig. 4A). To follow up on the pathway overrepresentation results, we selected consensus genes from the top-ranking five lipid metabolism pathways, and we examined their expression across the populations in the study. The heat map in Fig. 4B shows 28 lipid metabolism genes that are differentially expressed. Eighteen of these genes formed a cluster that was highly expressed in MacA cells but was expressed at low levels in the other populations. Thus, both the overrepresentation analysis and the analysis of expression of lipid pathway genes across MacA and MacB populations revealed that lipid metabolism is a specific gene expression signature of MacA cells.
Another pathway that ranked highly in cluster A was peroxisome proliferator receptor (PPAR) signaling (KEGG). A heat map visualizing all differentially expressed genes belonging to this pathway shows that out of 16 PPAR signaling genes, 12 were highly expressed in MacA cells (Fig. 4C).
We have previously shown that MacA cells in tumor-bearing lungs produce leukotrienes, which are a subfamily of lipid inflammatory eicosanoids derived from arachidonic acid (10). Thus, in Fig. 4D we examined the expression of known genes in the eicosanoid pathway (KEGG arachidonic acid metabolism) (33, 34), of which 16 were differentially expressed across the populations in this study. The key leukotriene metabolism genes 5-lipoxygenase (Alox5) and leukotriene C4 (LTC4) synthase (Ltc4s) were highly expressed in MacA cells from both tumor- and non–tumor-bearing lungs. In contrast, the key enzyme in PG production pathway, cyclooxygenase (COX)-2 (Ptgs2), was specifically downregulated in MacA cells compared with other populations. COX-1 (Ptgs1) was expressed at low levels in MacA cells from lungs without tumor; however, its expression was upregulated in MacA-2wk. For these key enzymes, we confirmed the mRNA-seq data in a separate group of mice using quantitative RT-PCR (Fig. 4E). To confirm the gene expression data functionally, we recovered MacA and Mac B cells from control and tumor-bearing mice by flow cytometry, stimulated these cells in vitro with a calcium ionophore, and measured eicosanoid production by liquid chromatography/tandem mass spectrometry (Fig. 4F). The production of LTC4 (a product of 5-lipoxygenase and LTC4 synthase) was specific to MacA cells, both in tumor-bearing and non–tumor-bearing lungs. Furthermore, the production of PGE2 (a product of Ptgs1/COX-1 or Ptgs2/COX-2) was increased in MacA cells from tumor-bearing mice versus control, which is in agreement with increased Ptgs1 expression by RNA-seq. Surprisingly, MacB cells produced low levels of PGE2, but this may be due to the nonphysiological stimulation. Among the receptors, MacA expressed high levels of Ptger2, the receptor for PGE2, and Fpr2, the receptor for lipoxins. In contrast, the receptor for LTB4 (Ltb4r1, BLT1) was specifically downregulated in MacA cells compared with other populations, as was the receptor for cysteinyl leukotrienes (Cysltr1).
Such a combination of enzymes and receptors suggests that among the studied populations, MacA cells would be the major producers of leukotrienes, but would not respond to these products. However, they would be responsive to lipoxins, which are the anti-inflammatory products of the lipoxygenase pathway, and to PGE2. The data also suggest that MacA cells in naive lung do not produce PGs; however, they gain that capability in the setting of tumor through the upregulation of COX-1. These data are consistent with increased levels of PGs seen in tumor-bearing mice, which are abrogated when cytosolic phospholipase A2, the rate limiting enzyme in eicosanoid production, is deleted (10, 35).
Comparing MacA-2wk to MacA-N, we identified 89 genes that were upregulated in tumor-bearing compared with control lungs, and 1 gene that was downregulated (Supplemental Table II). However, even though the differential expression was statistically significant, the absolute expression levels of the upregulated genes remained significantly lower in MacA-2wk than in MacB populations (Supplemental Fig. 1). It is possible that the increased levels of these genes reflect overall dysregulation of gene expression in MacA cells from tumor-bearing lungs, rather than biologically significant induction. Out of the 89 upregulated genes, only 5 (Ncam1, Slc40a1, Il6, Tsku, Ptgs1) had above average absolute expression levels in MacA-2wk cells compared with other populations (Fig. 4G). As discussed above, the increased expression of Ptgs1 (COX-1), which was confirmed by quantitative RT-PCR and by increased PGE2 production in vitro (Fig. 4E, 4F), indicates that MacA cells from tumor-bearing lungs may gain the potential to produce PGs. IL-6 has been implicated in many types of cancer, and increased levels are associated with poorer overall survival in lung cancer (36). We have previously demonstrated increased levels of IL-6 in our model and shown that cancer cells can induce expression in bone marrow–derived macrophages in vitro (35). Interestingly, IL-6 appears to be selectively expressed by MacA in the setting of tumor.
We next focused on MacB3 cells, which increase rapidly with tumor growth and constitute the major component of the lung TME, particularly at the 3 wk time point. As indicated by pathway overrepresentation analysis, cluster B3 (genes highly expressed both in MacB3-2wk and MacB3-3wk) was enriched in pathways related to chemokine and cytokine signaling (Fig. 5A, Supplemental Table I). To confirm this result independently of KEGG and Reactome databases, we chose a published panel of chemokine genes (37) and examined these in MacA and MacB populations. As shown in the heat map in Fig. 5B, of 21 differentially expressed chemokine genes, 16 clustered as highly expressed in MacB3s. Of those 16, 6 were also highly expressed in MacB2-3wk (Fig. 5B). Thus, the analysis of gene cluster B3 suggests that MacB3 cells play a critical role in communication between the diverse cells of the TME through the secretion of multiple chemokines. The other highly ranking pathways (Meiotic recombination, Systemic lupus erythematosus, Alcoholism, Chromosome maintenance) all contained 14 histone genes that overlapped with cluster B3, without any other genes pointing to a specific function.
In naive lung, the MacB3 cells are present in very low numbers, and thus we were unable to recover sufficient numbers of cells for RNA-seq analysis. However, these cells increased rapidly with tumor growth, between 2 and 3 wk. We identified 35 genes that were upregulated and 4 that were downregulated in MacB3 cells at 3 wk compared with 2 wk (Supplemental Table II). Overrepresentation analysis indicated enrichment of pathways related to extracellular matrix (ECM) (Fig. 5C). Interestingly, whereas upregulation of the ECM-related genes in MacB3 at 3 wk was statistically significant, these genes had the highest expression in MacB2 cells at 3 wk (Fig. 5D). This suggests that ECM production is upregulated in the TME at 3 wk, with MacB2 cells being the major producers, but MacB3 cells also contributing.
MacB2 and MacB1 cells
MacB2 cells also increase substantially during tumor growth. In our initial clustering of 2458 differentially expressed genes (Fig. 3), genes that were highly expressed in MacB2 at 3 wk formed a cluster separate from genes highly expressed in MacB2 from uninjected lung or from 2 wk tumor. In this cluster (B2-3wk), pathways that were highly ranked by overexpression analysis were related to ECM, particularly collagen (Fig. 6A, Supplemental Table I).
To follow up on this result, we examined a published collection of ECM proteins (38, 39) of which 61 were differentially expressed in this study (Fig. 6B). Interestingly, this analysis revealed two distinct gene clusters, one highly expressed in MacB2-3wk (34 genes) and one highly expressed in MacB1-N cells (21 genes). The MacB2-3wk cluster contained several genes previously associated with tumor progression, such as collagens, fibronectin, tenascin C, and versican (40, 41). As discussed above, these genes were also upregulated in MacB3 cells at 3 wk (Fig. 5D).
We examined changes in gene expression in MacB2 as a function of tumor progression. In comparison with MacB2-N, 87 genes were upregulated and 3 were downregulated in MacB2-2wk cells. However, most of these genes were not selectively expressed in MacB2 cells. Comparison of MacB2-N to MacB2-3wk found that 472 genes were upregulated and 97 were downregulated in MacB2-3wk cells (Supplemental Table II). Overrepresentation analysis indicated that the set of upregulated genes was enriched in pathways related to ECM and pathways related to cell cycle regulation (Fig. 6C, 6D). Because expression of ECM proteins was a distinctive feature of MacB2-3wk, this result is consistent with upregulation of these genes during tumor progression. The downregulated genes contained genes related to biological oxidation, and enriched pathways included Phase I—Functionalization of compounds, Drug metabolism—cytochrome P450, and Biological oxidations.
Analysis of genes highly expressed in MacB1 (cluster B1) or in MacB2 from uninjected lung and from 2 wk tumor (cluster B2) did not provide robust and repetitive results, such as those observed for MacA, MacB2-3wk, and MacB3 cells.
Relationship to M1/M2 macrophages
Because studies focusing on clinical outcomes frequently use the M1/M2 designation of macrophages, we examined the relationship of our populations with commonly used M1 and M2 markers. This analysis has shown that ArgI (M2 marker) was expressed at highest levels in MacB3 (at both 2 and 3 wk), whereas Nos2 (M1 marker) was expressed in both MacB2 (2 and 3 wk) as well as MacB3 (2 and 3 wk), with the highest levels in MacB2-3wk. Analysis of a broader set of M1/M2 markers (42) indicated that none of our populations preferentially expresses either M1 or M2 markers (Supplemental Fig. 2). This analysis demonstrates the complexity of macrophage phenotypes in vivo and poor correlation with the M1/M2 designation, which was based on in vitro activation of bone marrow–derived macrophages (42).
Correlation with genes prognostic in human lung adenocarcinoma
Finally, a recent study has used a bioinformatics approach to establish a system to examine associations between gene expression in human tumors and patient survival across multiple types of cancer (PRECOG; https://precog.stanford.edu/) (43). We used the data from the PRECOG database to test whether MacA or MacB populations are enriched in genes correlated with good or bad clinical outcomes in lung adenocarcinoma. We ranked ∼23,000 human genes according to their survival scores, and we used GSEA (preranked tool) (23, 24) to examine whether the clusters associated with MacA or MacB populations (defined in Fig. 3) are enriched in genes correlated with good or bad prognosis. The GSEA plots and normalized enrichment scores (Fig. 7) indicate that MacB2-3wk and MacB3 (-2wk and -3wk) are enriched in genes correlated with bad clinical outcomes. In contrast, MacB1-N, MacA, and the early stage MacB2 (-N and -2wk) were enriched in genes predicting good clinical outcomes, with MacB1-N having the highest enrichment score. These data suggest that prognostic genes identified in the analysis of whole tumors are highly expressed in specific macrophage populations, suggesting that enrichment of these populations in a tumor may predict clinical outcome.
Macrophages have been implicated as critical components of the TME, controlling tumor initiation, progression, and metastasis. In this study we have taken an unbiased approach to begin to define the complexity of macrophage phenotypes in the setting of lung cancer progression. We have employed an immunocompetent orthotopic model of cancer progression in which interactions between macrophages, cancer cells, and adaptive immune cells are present. Our study has defined at least four separate monocyte/macrophage populations and used transcriptional profiling to begin to define their function. Our results reveal that each population has a distinct gene expression signature, suggesting that it plays a specific role in cancer progression.
Our data underscore that the TME is shaped both by increases in specific cell populations, as well as by gene expression changes in each of these populations (Fig. 8). The population designated MacA represents resident alveolar macrophages. This is confirmed by correlating these cells with the ImmGen signature. This population does not change significantly in number during tumor progression. However, there are modest but potentially important changes in gene expression associated with the presence of the tumor. The gene expression signature implies that this MacA population may be the major mediator of leukotriene production, and the induction of COX-1 suggests that these cells also contribute to the increase in PG production observed in lung cancer progression (10). Although we confirmed the gene expression data by measuring LTC4 and PGE2 production in MacA cells, one limitation of this experiment is that the cells were analyzed outside of their physiological environment. Another gene selectively expressed in the MacA from tumor-bearing lungs is IL-6, which plays a complex role in cancer progression.
Cells that we have designated as MacB, especially MacB2 and MacB3, represent populations that increase during tumor progression in our model. Our analysis suggests that these cells are involved in specific functions regulating tumor progression. In the absence of tumor, these cells are present in very low numbers in the lung, suggesting that they are recruited from the circulation in response to factors produced by cancer cells. MacB3 cells, which bear markers of infiltrating macrophages, expressed high levels of multiple cytokines at both time points examined. This would suggest that these cells may orchestrate the recruitment of other stromal cells during tumor progression. This may induce a feed-forward loop, in which early recruitment of MacB3 cells leads to additional waves of recruitment of other stromal cells. Interestingly, at the later time point, we observed increases in genes encoding ECM and enzymes that remodel the ECM. Similar ECM-related pathways were also detected in MacB2 cells, which represent tumor-associated monocytes, and which show similar increases during tumor progression. Thus, these distinct populations of macrophages can alter the TME both by recruiting additional cells and by remodeling the physical properties of the TME.
Correlating our data with analysis of human tumors, it appears that genes associated with better clinical outcomes are expressed at earlier time points, whereas genes associated with poor outcomes are expressed at later time points. Although studies have suggested that these changes are a result of a phenotypic switch in macrophage populations interacting with cancer cells, our data indicate a more complex situation. Distinct populations of macrophages expand during cancer progression. Thus, a change in a specific marker will reflect alterations in existing populations due to phenotypic modulation, as well as increases in populations that express these markers constitutively. For example, studies have demonstrated increased production of TGF-β1 by macrophages during cancer progression. From our analysis, the level of TGF-β1 production is not altered in any of our populations (not shown), but marked increases in the numbers of MacB2 and MacB3 cells during progression will result in increased numbers of TGF-β1–producing cells. However, increases in several ECM proteins will be a result of increased numbers of ECM-producing macrophages, as well as upregulation of production in these cells during tumor progression.
The populations we have defined do not easily conform to the M1/M2 paradigm. Although some M1/M2 markers are enriched in certain populations, using a more extensive list of markers shows that these genes are expressed across multiple populations. Additionally, a model in which macrophages recruited to the site of the tumor as M1/antitumorigenic macrophages undergo phenotypic changes to an M2/protumorigenic phenotype is not supported by our data. Based on the PRECOG analysis, most of the populations we have defined retain an expression signature associated with either good or bad clinical outcomes throughout tumor progression. An exception is the MacB2 population, which appears to be associated with good outcomes at early time points, but undergoes changes in gene expression leading to an association with a bad outcome at later time points.
To our knowledge, this is one of the first studies to molecularly define the diversity of macrophages in lung cancer in an immunocompetent setting. The major conclusion of our analysis is that distinct macrophage populations express specific gene signatures, which predicts that they play specific and distinct roles in cancer progression. However, the actual function of these cells will need to be validated in experiments that selectively deplete a specific macrophage subpopulation in the TME. Our studies indicate that recruited populations of macrophages generally aid tumor progression, suggesting that selectively blocking their recruitment may represent a therapeutic strategy to inhibit lung cancer progression. A limitation of this study is the use of only one lung cancer model. However, we would argue that our orthotopic injection model is an appropriate model to answer questions concerning progression of advanced lung cancer and the role of specific cells in the lung TME. Other models would less accurately reflect the TME of the lung (e.g., s.c. flank injections) or would not result in a well-defined primary tumor (e.g., tail vein injections). Another limitation of these studies is the use of a single lung cancer cell line. Targeted therapies in lung adenocarcinoma are based on identification of specific oncogenic drivers that drive proliferation of cancer cells in a cell-autonomous fashion. However, there is little known as to how specific drivers affect the TME. It will be of interest to determine whether the changes in number and gene expression profiles of the macrophage populations we have defined are dependent on the oncogenic drivers of the cancer cells. For our orthotopic model, this will require a panel of murine lung cancer cells that express the drivers that have been identified in human lung cancer. Thus, our report provides a framework for studies of macrophage populations in other models of lung cancer, including spontaneous models in which specific oncogenic drivers are activated.
Finally, it is important to extend these studies to human lung tumors. In this regard, it is encouraging that we were able to correlate transcriptional profiles of specific macrophage populations with either bad or good prognostic genes derived from bioinformatic studies in human lung adenocarcinoma. Thus, our studies form the basis to develop a panel of markers to characterize macrophage populations in human tumors. Earlier studies have shown that the location of macrophages in lung tumors may determine their role in cancer progression. Increased numbers of macrophages within the tumor itself was associated with a favorable prognosis, whereas increased numbers of macrophages in the stroma surrounding the tumor was associated with an unfavorable prognosis (44). The development of a panel of markers specific for MacB3 and MacB2 cells will be required to define their localization, and may be of value in predicting clinical outcomes in human disease.
We thank the members of the University of Colorado Cancer Center Flow Cytometry Core for assistance with flow cytometry and the Genomics and Microarray Core for assistance with RNA-seq.
This work was supported by National Institutes of Health/National Cancer Institute Grants R01 CA162226, R01 CA164780, and R01 CA108610 (to R.A.N.); the Colorado Lung Specialized Program of Research Excellence (National Cancer Institute Grant P50 CA058187); Ruth L. Kirschstein National Research Service Award T32CA17468 (to T.R.S.); and by Department of Veterans Affairs Grant CDA 1IK2BX001282-01A1 (to H.L.). The Flow Cytometry and the Genomics and Microarray Shared Resources receive support from the National Institutes of Health/National Cancer Institute (University of Colorado Cancer Center Support Grant P30 CA046934). S.D. acknowledges support from the Boettcher Foundation and the United Against Lung Cancer Foundation.
The RNA-seq data presented in this article have been submitted to the National Center for Biotechnology Information Gene Expression Omnibus repository (http://www.ncbi.nlm.nih.gov/geo/) under accession number GSE76033.
The online version of this article contains supplemental material.
Abbreviations used in this article:
fragments per kilobase of exon per million fragments mapped
gene set enrichment analysis
Immunological Genome Project
Lewis lung carcinoma
- MHC II
MHC class II
from naive mice
peroxisome proliferator–activated receptor
from 2 wk tumor-bearing mice
from 3 wk tumor-bearing mice.
The authors have no financial conflicts of interest.