Phagocytosis is a key macrophage function, but how phagocytosis shapes tumor-associated macrophage (TAM) phenotypes and heterogeneity in solid tumors remains unclear. Here, we utilized both syngeneic and novel autochthonous lung tumor models in which neoplastic cells express the fluorophore tdTomato (tdTom) to identify TAMs that have phagocytosed neoplastic cells in vivo. Phagocytic tdTompos TAMs upregulated antigen presentation and anti-inflammatory proteins, but downregulated classic proinflammatory effectors compared to tdTomneg TAMs. Single-cell transcriptomic profiling identified TAM subset-specific and common gene expression changes associated with phagocytosis. We uncover a phagocytic signature that is predominated by oxidative phosphorylation (OXPHOS), ribosomal, and metabolic genes, and this signature correlates with worse clinical outcome in human lung cancer. Expression of OXPHOS proteins, mitochondrial content, and functional utilization of OXPHOS were increased in tdTompos TAMs. tdTompos tumor dendritic cells also display similar metabolic changes. Our identification of phagocytic TAMs as a distinct myeloid cell state links phagocytosis of neoplastic cells in vivo with OXPHOS and tumor-promoting phenotypes.
Macrophages are innate immune cells that play crucial homeostatic, defensive, and reparative roles in all tissues of the body (Gessain et al., 2020). These cells respond to tissue environmental cues and possess a spectrum of activation states that contribute to their considerable heterogeneity (Mosser and Edwards, 2008). A key function of macrophages is to phagocytose extracellular material including pathogens, dead/dying cells, and neoplastic cells (Metchnikoff, 1883; Watanabe et al., 2019). The phagocytosis of apoptotic cells, termed efferocytosis, is critical for homeostasis, but can also suppress immune surveillance and promote tumor growth (A-Gonzalez et al., 2009; Hanayama et al., 2004; Mukundan et al., 2009; Roszer et al., 2011; Scott et al., 2001; Uribe-Querol and Rosales, 2020). Phagocytic macrophages at steady state have been associated with increased cytoskeletal rearrangements, phago-lysosomal activity, antigen presentation, and lipid metabolism (A-Gonzalez et al., 2017). Additionally, pathway analyses relating to phagocytosis have helped identify putative phagocytic macrophages across healthy and diseased tissue states (Jaitin et al., 2019; Lantz et al., 2020; Sanin et al., 2022). Despite these efforts, the functional identification and validation of phagocytic immune cells in the tumor microenvironment remains elusive, and thus it is unclear whether the engulfment of neoplastic cells by tumor-associated macrophages (TAMs) contributes to immune-permissive or immune-suppressive macrophage cell states.
TAMs are one of the most abundant immune cell types in tumors, and their numbers often inversely correlate with clinical outcome (Cassetta and Pollard, 2018). In lung cancer, TAMs can be derived from both tissue resident macrophages, such as alveolar macrophages (AMs) of primarily embryonic origins, and monocyte-derived macrophages, such as interstitial macrophages (IMs; Casanova-Acebes et al., 2021; Zhu et al., 2017). Understanding TAM heterogeneity has been an area of intense interest; however, how macrophage functions such as phagocytosis alter their phenotypes, the tumor microenvironment, and ultimately clinical outcomes remain poorly understood.
Macrophage metabolism greatly influences their phenotype and function (Koo and Garg, 2019; Wculek et al., 2022; Yan and Horng, 2020). In vitro treatment of macrophages with defined factors can either induce a proinflammatory/anti-tumor (M1) phenotype, which is dependent on glycolytic metabolism, or an anti-inflammatory/pro-tumor (M2) phenotype that is characterized by a more flexible metabolism, utilizing both glycolysis and mitochondrial oxidative phosphorylation (OXPHOS; Viola et al., 2019). The OXPHOS pathway is responsible for synthesizing and/or catabolizing metabolites or macromolecules like amino acids or lipids. Lipid accumulation in the cell, derived for example from phagocytosed material, needs to be either stored, effluxed, or catabolized (Yurdagul, 2021). During wound healing, fatty acids derived from phagocytosed cells can lead to increased OXPHOS to fuel reparative macrophage phenotypes (Zhang et al., 2019). Additionally, efferocytosis can induce macrophages to secrete suppressive cytokines (Yin and Heit, 2021). Therefore, efferocytosis and the consequent metabolic reprogramming of macrophages may play an important role in driving the phenotypes of TAMs in vivo.
TAMs are often described as having a wound-healing, anti-inflammatory phenotype, and the identification of functional and metabolic cell states of macrophages has improved our understanding of TAMs beyond the M1/M2 paradigm (Katzenelenbogen et al., 2020; Loke and Lin, 2022; Mulder et al., 2021; Sanin et al., 2022). Connecting the phagocytic and metabolic states of TAMs in the context of a growing tumor may uncover novel biological mechanisms in TAMs that could inform therapeutic design. Historically, studies of macrophage metabolism have employed short-term, in vitro assays, often using bone marrow–derived macrophages (BMDMs), and it is unclear whether these systems are relevant to TAMs in vivo (A-Gonzalez et al., 2017; Schulz et al., 2019; Zhang et al., 2019). Studies of macrophage identity and heterogeneity have shown that microenvironmental cues imprint phenotypic changes after chronic tissue residency (Guilliams et al., 2020; Gundra et al., 2017; Yona et al., 2013). Autochthonous mouse models, which rely upon engineered genetic mutations to drive tumor development in tissues, reproduce many features of human cancer (Singh et al., 2012). Importantly, these models provide unique systems to study how signals, like phagocytosis of neoplastic cells, influence TAM heterogeneity.
Here, we employed syngeneic and autochthonous mouse models of lung tumorigenesis to study the impact of phagocytosis of neoplastic cells on TAMs. We identify phagocytic TAMs as a distinct subset of macrophages in vivo, show that phagocytosis of neoplastic cells reprograms macrophage metabolism and sustains immune-suppressive TAM phenotypes, and uncover a unique phagocytic tumor (PHAT) signature enriched in OXPHOS genes. This signature is shared among monocytes, macrophages, and dendritic cells in murine lung tumors, and is enriched within phagocytic macrophages in previously described models of inflammation and cancer. The PHAT signature is enriched in lung tumor myeloid cells and correlates with poorer clinical outcome in several human cancer types. These expression signatures were validated by protein-level phenotyping and functional metabolic assays. Interestingly, non-phagocytic TAMs were phenotypically similar to macrophages in normal tissues, suggesting that engulfment of neoplastic cells is a major driver of TAM phenotypes. This work connects a key cellular function of TAMs—the phagocytosis of neoplastic cells—with changes in intracellular metabolism and TAM heterogeneity in lung tumors.
Phagocytic TAMs have distinct phenotypes from non-phagocytic TAMs in a syngeneic transplant lung tumor model
To initially investigate the phagocytosis of lung cancer cells by macrophages in vivo, we intravenously injected an oncogenic KrasG12DTrp53-deficient murine lung cancer cell line that expresses the red fluorescent protein tdTomato (tdTom) to seed orthotopic tumors in the lung (KPT cell line; Fig. 1 A). We identified CD45pos tdTompos immune cells as those that had recently phagocytosed cancer cells (Fig. 1 B). The frequency of CD45pos tdTompos cells significantly increased over time as did tdTom mean fluorescence intensity (MFI) in tumor-infiltrating CD45pos cells (Fig. 1, B and C). By 3 wk after transfer, alveolar macrophages (AMs; defined as CD64pos CD24neg CD11bneg SIGLEC-Fpos cells) and CD11bpos interstitial macrophages (IMs; defined as CD64pos CD24neg CD11bpos SIGLEC-Fneg cells) were the most abundant phagocytic cell type. IMs were the most numerous among tdTompos CD45pos myeloid cells, which is consistent with previously described models (Maier et al., 2020; Fig. 1 D). CD11bneg dendritic cells (DCs) and monocytes were also tdTompos, albeit at low frequencies. As expected, non-phagocytic B and T lymphocytes were tdTomneg within tumors (Fig. S1 A). Thus, in this syngeneic lung tumor model, IMs are the most prevalent phagocytes of neoplastic cells.
The expression of several key immune modulatory proteins (PD-L1 and CD80) and phagocytic receptors (MERTK and ITGAV) was quantified on TAMs by flow cytometry (Fig. 1, E and F). CD80 expression on tdTompos macrophages significantly increased over time on both AMs and IMs. tdTompos IMs had significantly higher expression of the immune checkpoint ligand PD-L1 than tdTomneg IMs at 2 and 3 wk after cell line transplantation. In contrast, PD-L1 was downregulated in tdTomneg AMs 3 wk after transfer, relative to naive AMs. Interestingly, MERTK and ITGAV expression on AMs followed a similar trend to PD-L1, in which tdTomneg AMs down-regulated cell surface expression. In contrast, the expression of MERTK and ITGAV on IMs steadily increased over time, thus highlighting key differences between tissue resident AMs and monocyte-derived IMs. Monocyte-derived macrophages that have phagocytosed KPT cells in vitro also upregulate CD80, MERTK, ITGAV, and PD-L1, consistent with the in vivo IM phenotypes (Fig. S1 B).
We next profiled tdTompos and tdTomneg AMs and IMs using a panel of mRNA probes for genes related to innate myeloid immunity. Phagocytic tdTompos IMs and AMs downregulated genes related to TLR signaling, cytokine signaling and T cell activation (Fig. 1, G and H). Relative to tdTomneg IMs, tdTompos IMs had higher expression of genes related to extracellular matrix (ECM) remodeling and lower expression of those related to angiogenesis. Conversely, tdTompos AMs had positive enrichment of both angiogenesis and ECM remodeling (Fig. 1, I and J). Together, these results suggest that phagocytosis of neoplastic cells alters TAM cell surface and transcriptional phenotypes in vivo in a manner distinct from those not actively engaged in phagocytosis, and that these phenotypes differ between IMs and AMs.
Tissue-resident AMs are the most prevalent phagocyte in autochthonous lung tumors
Genetically engineered mouse models recapitulate the genetics and histology of human lung tumors. Importantly, the slower growth of these tumors allows the microenvironment to develop over time, which more closely resembles the long-term tumor residency of macrophages in human tumors. Previous viral-Cre–induced lung cancer models have documented transduction of lung-resident AMs and consequently modification of conditional allele by Cre recombinase in AMs (Tippimanchai et al., 2018). Thus, we developed a novel model that specifically induces tdTomato protein expression and oncogenic alterations in neoplastic cells to enable tracking of phagocytic activity (Fig. 2, A and B). This autochthonous lung tumor model specifically induces genetic alterations in type II alveolar cells which are a cell type of origin of lung adenocarcinoma (Sutherland et al., 2014). Transduction of KrasFSF-G12D;Rosa26FSF-LSL-tdTomato(ai65);SftpcCreER (KTai65;SpcCreER) mice with a lentiviral vector containing an inverted FLPo flanked by pairs of heterotypic loxP sites enables tamoxifen-induced expression of oncogenic KRAS specifically in type II alveolar cells (Fig. 2 A; unpublished data). Nine months after tumor initiation, KTai65;SpcCreER mice developed numerous tdTomato-positive lung tumors. Phagocytic tdTompos myeloid cells could be visualized by flow cytometry similar to the cell transfer model (Fig. 2 B). Importantly, unlike in the Lenti-Cre initiated autochthonous lung tumor models (Fig. S1, C–H), tdTombright myeloid cells were not observed, consistent with recombination events only occurring in epithelial cells (Fig. 2 C).
Within the KTai65;SpcCreER model, the majority (∼70%) of tdTompos myeloid cells were AMs, whereas only 3% were IMs (Fig. 2 D). This suggests contrasting properties in myeloid surveillance in syngeneic transplant models, in which IMs accounted for 25–50% of the phagocytic myeloid population across the 3 wk of lung tumor growth (Fig. 1 D). Compared to the syngeneic model, similar trends were seen in the cell surface expression of PD-L1, CD80, MERTK, and ITGAV on tdTompos AMs and IMs; however, the fold-change differences were higher in IMs from syngeneic tumor models compared to the KTai65;SpcCreER model (Fig. 2, E–H; and Fig. 1, E and F). Moreover, tdTompos AMs from KTai65;SpcCreER tumors most closely resembled the tdTompos AMs from syngeneic tumors observed at 2 wk after transfer.
When lung tumors were initiated with Adeno-Spc-Cre in KrasLSL-G12D;Trp53fl/fl;Rosa26LSL-tdTomato (KPT) mice, there were also very few tdTombright AMs compared to when mice are transduced with non-specific viral vectors (Fig. S1, I and J). In these tumor-bearing KPT mice, there was only a slightly higher percentage of phagocytic AMs compared to IMs (Fig. S2, A–C). Thus, it appears that IMs play a larger phagocytic role in these tumors compared to tumors in the KTai65;SpcCreER model. However, the expression of key cell surface proteins on AMs and IMs closely resembled the phenotypes in the KTai65;SpcCreER model (Fig. S2, D and E). Collectively, the data from these models highlight the heterogeneity of phagocytic myeloid cells in the lung tumor microenvironment and provide insight into the phenotypes that correlate with the phagocytosis of neoplastic cells in vivo.
Phagocytosis of neoplastic cells imprints distinct transcriptional patterns in monocytes and macrophages
To investigate cell state changes that occur in phagocytic myeloid cells in lung tumors, we microdissected tumor and uninvolved tissue from KTai65;SpcCreER tumor-bearing mice, as well as lung tissue from tumor-naive mice, and performed single-cell RNA-sequencing (scRNA-seq) on sorted tdTompos and tdTomneg myeloid cells (enriched for myeloid cells; CD45pos CD11bpos and/or CD11cpos) and CD45neg EPCAMpos non-immune cells (Fig. 3 A and Fig. S3, A and B). Cell hashtagging using antibodies with conjugated oligonucleotide barcodes was used to multiplex cells from different animals and tissue regions, filter singlets, and enable direct comparisons between tdTompos and tdTomneg cells from tumor and uninvolved tissues (Fig. 3 B). After implementing a multi-step procedure for filtering transcriptome quality and removing multiplets (Fig. S3, C–G), we recovered 44,016 high-quality single-cell transcriptomes that could be assigned tissue and animal identities based on hashtag labeling (Fig. 3, C–E).
Graph-based clustering after optimizing cluster resolution identified 33 distinct cell states with unique expression profiles, including two actively proliferating cell states (Fig. S3, H–J). These grouped into 16 major cell types from myeloid, lymphoid, endothelial, or epithelial origins based on known cell markers (Fig. 3 F). We identified seven cell types of myeloid origin (Fig. 3 G). Although the identified myeloid cell types were present in all tissue types, AM, IM, monocyte, and DC cell types were enriched in sorted tdTompos myeloid cells from tumors (Fig. 3 H). This corroborated our flow cytometry results and underscores the role of these cell types in neoplastic cell engulfment.
To elucidate how cells that have engulfed neoplastic cells are transcriptionally distinct from other non-phagocytic myeloid cell types, we performed differential gene expression analyses comparing tumor-associated tdTompos monocytes, IMs, and AMs to their tumor-associated tdTomneg counterparts (Fig. 3 I). Genes involved in phagocytosis, metabolism, and co-stimulation were generally upregulated in tdTompos cells of all three types, although some genes associated with these functions were only significantly upregulated in IMs (e.g., Axl, Cpt1a, and Cd40; Fig. 3 J). Antigen presentation genes were upregulated in tdTompos AMs and IMs but not in tdTompos monocytes. Finally, certain mediators of inflammation (e.g., Il1b, Tnf) and tissue repair (e.g., Mgl2, Chil3) had contrasting expression patterns across cell types. These data suggest that phagocytosis of neoplastic cells induces conserved and distinct transcriptional programs in IMs, AMs, and monocytes in autochthonous lung tumors.
PHAT gene signature is associated with OXPHOS metabolism and translation
Next, we sought to identify a gene signature representing common transcriptional changes associated with phagocytosis of neoplastic cells. Because macrophages are the most prevalent tdTompos myeloid cell in KTai65;SpcCreER tumors and monocytes are the precursor to IMs, we analyzed the transcriptional changes occurring in AMs, IMs, and monocytes. tdTompos AMs, IMs, and monocytes shared a common gene signature of 123 upregulated genes, which we term the PHAT signature (Fig. 4 A and Tables S1, S2, S3, S4, and S5). A comparison of pathways from the Kyoto Encyclopedia of Genes and Genomes (KEGG) database with the PHAT signature identified significant enrichment for pathways relating to ribosomal translation, OXPHOS, phagosome, and metabolic pathways (Fig. 4 B). Additionally, Ingenuity Pathway Analysis revealed shared activation of the OXPHOS, phospholipases, glycolysis, and gluconeogenesis metabolic pathways in tdTompos macrophages and monocytes (Fig. 4 C). We also defined a 233-gene signature to uncover pathways shared by both macrophage populations, independent of monocytes, termed the PHAT-Macrophage signature. This analysis detected antigen processing/presentation and lysosomal pathways in addition to pathways identified from the overall PHAT signature (Fig. 4 D). Macrophage metabolism is known to be a key driver of the phenotypes and functions of TAMs and a source of TAM heterogeneity within the tumor microenvironment (Li et al., 2022; Wang et al., 2021; Xia et al., 2020). Our data identify a common metabolic gene signature of phagocytic monocytes and macrophages that have recently phagocytosed neoplastic cells in vivo.
To discriminate whether these genes originate from the phagocytes themselves or internalized neoplastic cells, we performed bulk RNA-seq on tdTompos AMs and tdTompos neoplastic cells from KTai65;SpcCreER tumors. Known marker genes for alveolar macrophages and type II alveolar cells/neoplastic cells showed expected cell-specific expression patterns. Importantly, tdTompos AMs had extremely low (>400-fold lower) relative expression of type II/neoplastic cell markers, suggesting there is minimal contribution of mRNA from engulfed neoplastic cells (Fig. 4 E). Therefore, the PHAT signature is comprised mainly of genes expressed by phagocytes.
Comparison of PHAT signatures to other in vivo macrophage datasets shows both conserved and context-dependent phagocytic macrophage states
Since our model enables characterization of phenotypic changes of immune cells following the in vivo engulfment of tdTompos neoplastic cells, we compared our functionally defined PHAT signatures to several studies of macrophage heterogeneity related to phagocytosis. One study directly follows efferocytosis-driven gene expression changes in peritoneal macrophages after isolation and co-culture with apoptotic cells; phagocytosis-dependent changes were defined by using a mouse model that expresses a kinase-deficient form of Mertk, termed MerTK−/−, thus preventing apoptotic cell uptake (Lantz et al., 2020). In this dataset, our PHAT signature score was higher in both the monocyte-derived F4/80neg MHCIIhi and tissue-resident F4/80pos MHClo macrophages co-cultured with apoptotic cells. However, MerTK−/− macrophages had decreased the PHAT scores, demonstrating that our signature is significantly upregulated in efferocytic macrophages (Fig. 4 F). Another study characterized monocyte-derived macrophage heterogeneity across a broad set of inflammatory states and defined certain cell clusters as “phagocytic” based on in vitro–defined key pathway genes (Sanin et al., 2022). We focused our signature to include only genes upregulated in the monocyte-derived IMs comparing tdTompos IMs relative to tdTomneg IMs, termed the PHAT-IM signature. In two different analyses, these putative phagocytosis-related clusters show enrichment for our PHAT-IM signature, but other differently labeled clusters also show enrichment (Fig. 4 G). This result highlights the importance of defining phenotypes functionally and suggests that although some pathways are conserved, the tumor context may drive phagocytic macrophages to resemble activation states in other inflammatory settings. Interestingly, a study of adipose tissue macrophages defined populations of lipid-rich Trem2pos macrophages in mouse (Mac3) and human (LAM; Jaitin et al., 2019). The PHAT-IM signature was highly enriched in LAMs compared to macrophages with low lipid content (Fig. 4 H). Additionally, flow cytometry revealed increased TREM2 expression on tdTompos AMs and IMs from autochthonous lung tumors (Fig. S2, D and E). Together, these comparisons show that gene expression changes in TAMs are shared in other non-tumor contexts, although certain changes are uniquely imprinted by the mechanism of macrophage phagocytosis, tissue milieu, and inflammatory status.
PHAT gene signature translates to human TAMs and clinical outcome in lung cancer
To investigate the clinical relevance of these phagocytic myeloid cells, we queried The Cancer Genome Atlas (TCGA) bulk gene expression datasets. Given the immune suppressive phenotypes of phagocytic TAMs, we anticipated that our PHAT signature could negatively correlate with clinical outcomes. We found that our PHAT signature could stratify patients with lung adenocarcinoma, with patients whose tumors were enriched for the PHAT signature having significantly decreased probability of survival (Fig. 4 I). Breast cancer, hepatocellular carcinoma, and colon adenocarcinoma patients with high PHAT signature had only slightly lower survival; however, the PHAT signature significantly stratified head and neck squamous cell carcinoma patients (Fig. 4 J and data not shown). These data suggest that certain molecular features are shared, although idiosyncrasies exist in different tumor settings.
To query our PHAT signature specifically within TAMs from human cancer, we investigated the enrichment of our PHAT signature in bulk RNA-seq dataset on sorted myeloid cells from tumors and normal lung tissue from non-small cell lung cancer patients (Combes et al., 2022). Importantly, the PHAT signature was significantly enriched in tumor-associated myeloid cells compared to myeloid cells from normal lung (Fig. 4 K). This contrasted with classical M1 and M2 myeloid gene signatures, which were not significantly enriched in tumor-associated myeloid cells compared to those from normal lung. Together, these findings demonstrate that our murine-derived PHAT signature associates with the gene expression of human tumor myeloid populations and connects the transcriptional changes that occur in phagocytic monocytes and macrophages in autochthonous mouse tumors with clinical outcome.
Phagocytic AMs mature through a distinct lineage associated with tissue remodeling
We next evaluated how phagocytosis influences AM maturation at single-cell resolution using trajectory inference and cell state composition analyses across tissue and phagocytic (tdTompos) status. Cell embedding using diffusion maps contextualized local maturation relationships between the six identified AM cell states within the lung microenvironment (Fig. 5 A and Fig. S4 A; Haghverdi et al., 2015). This embedding revealed that naive and tdTomneg AMs from uninvolved tissue and tumors grouped separately from tdTompos tumor-associated AMs, underscoring the unique transcriptional cell states driven by phagocytosis of neoplastic cells (Fig. 5 B). Furthermore, naive, uninvolved, and tumor-associated tdTomneg had a similar distribution of AM cell states, while the tumor-associated tdTompos AMs had a high proportion of AM-2 cells (Fig. S4 B). Trajectory inference using Slingshot predicted three trajectories of AMs (termed AM-Trajectory A, AM-Trajectory B, and AM-Trajectory C; Fig. 5 C; Street et al., 2018), and dynamic RNA velocity modeling confirmed the relative progression across AM cell states (Fig. S4 D). The majority of tdTompos tumor-associated AMs were positioned relatively late along the maturation continuum and were found in AM-Trajectory A (Fig. 5, D and E). Interestingly, the presence of tumor-associated tdTomneg AMs (but not AMs from normal lung) at later pseudotime values suggests that these cells may be phagocytic AMs that have degraded tdTomato.
To define the transcriptional phenotypes of tdTompos AMs between trajectories, we analyzed the position of AM cell states across pseudotime (Fig. 5 F). All three AM-Trajectories initiate from the AM-1 state (Fig. 5 G). While AM-Trajectory A and AM-Trajectory B then differentiate toward the AM-2 state (Plet1high and MHC class IIhigh) before diverging toward AM-3/AM-4 (tissue remodeling) or AM-5/AM-6 (inflammasome activation/cytoskeletal organization) states, AM-Trajectory C, which has the highest fraction of AMs from naive and uninvolved lung, directly transitions from AM-1 to AM-4. AM-4 appears in AM-Trajectory A and C, and differential expression testing identified higher immunosuppressive genes in AM-4 cells from Trajectory A, suggesting that convergence of AM maturation may arise from different environmental stimuli (Fig. S4 E).
Cells from AM-Trajectory A had the highest expression of PHAT signature genes, revealing shared molecular programs between different types of phagocytic macrophages (Fig. 5 H). Also, cells from AM-Trajectory A expressed higher levels of tissue remodeling regulators (e.g., Fabp5, Lgmn, Mmp12, Trem2, Calr, Il13ra1) and M2-like genes (e.g., Tgfbi, Chil3, Cd274) compared to AM-Trajectory B and AM-Trajectory C (Fig. 5 I and Fig. S4 F). This downregulation of inflammasome genes (e.g., Tnf, Neat1, Nlrp3, Jun, Il18) and upregulation of antigen tissue remodeling genes is indicative of an M2-like skewing of TAMs correlating with phagocytosis. However, TAMs within AM-Trajectory A also upregulate antigen presentation genes commonly associated with M1 macrophage phenotypes. This is in line with mixed M1/M2 TAM phenotypes commonly observed in human TAMs (Larionova et al., 2020). Metabolic genes, such as Atp5a1, a subunit of Complex V of the electron transport chain, and Slc2a1, which encodes the glucose transporter GLUT1, have high, sustained expression throughout AM maturation in Trajectories A and C but not in Trajectory B, suggesting utilization of both OXPHOS and glycolytic metabolism by cells from AM-Trajectory A and AM-Trajectory C (Fig. S4 F). In summary, trajectory analysis of AMs, combined with fluorescent tracking of neoplastic engulfment, highlights the diversity of AM responses and suggests that neoplastic engulfment promotes AM cell states with increased metabolic gene expression, tissue repair, and anti-inflammation properties as well as an inability to sustain proinflammatory gene expression.
Phagocytic IMs diverge from monocytes and exhibit plasticity after engulfment
To understand how phagocytosis influences monocyte and IM functions, we embedded these populations in a lower-dimensional space using diffusion maps and then performed trajectory inference (Fig. 6, A and C). tdTompos monocytes/IMs accumulated separately from cells in naive and uninvolved tissues along the differentiation continuum, suggesting that differential transcriptional regulation occurs in cells that have phagocytosed neoplastic cells. Tumor-associated IMs, which could be further classified into two cell states based on the expression of complement system genes (e.g., C1qa and C3ar1) or oxidative stress genes (e.g., Ftl1 and Mt1; Fig. S4 G), accounted for a higher proportion (∼80%) of tumor tdTompos monocytes/IMs compared to monocytes (∼20%; Fig. S4, H and I). Trajectory inference using Slingshot predicted two trajectories (Mono-IM-Trajectory and Mono-Trajectory) and ordered cells by pseudotime. The directionality of the trajectory was confirmed using RNA velocity (Fig. S4 J). The Mono-IM-Trajectory contained a significantly higher percentage of tdTompos cells compared to the Mono-Trajectory (Fig. 6 D). Like tumor-associated tdTompos AMs, the majority of tdTompos monocytes/IMs were associated with later pseudotime values (Fig. 6, E and F).
Both trajectories initiate from the Cd300ehigh monocyte state (Monocyte-1) and an intermediate Monocyte-2 cell state, followed by a bifurcation into a terminal Monocyte-3 state (Mono-Trajectory) or into two IM states (Mono-IM-Trajectory; Fig. 6, F and G). Trajectory inference predicted that IM-1, which accounts for the majority (∼65%) of tumor-associated phagocytic monocytes/IMs, precedes IM-2 in maturation and is less widely distributed across non-tumor tissues (Fig. 6 G and Fig. S4 I). The PHAT signature gene expression was highest in cells along the actively phagocytic Mono-IM trajectory (Fig. 6 H). Genes associated with phagocytosis were more highly expressed in the Mono-IM-Trajectory compared to the Mono-Trajectory, while genes associated with migration were more highly expressed in the Mono-Trajectory (Fig. 6 I). Markers of phagocytosis (e.g., C1qa, Mrc1, Trem2) and glucose metabolism (e.g., Atp5a1, Slc2a1, Cox8a) also correlated with the highest enrichment of tumor-associated tdTompos monocytes/IMs (Fig. S4 K), but expression of these genes was downregulated at later pseudotime values. In contrast, phagocytic AMs have persistent upregulation of the same genes. We observed a similar transient upregulation of pro-inflammatory (M1-like) markers (e.g., Tnf and Cd74) across pseudotime, though mixed expression patterns were observed in anti-inflammatory (M2-like) genes, as Chil3 was upregulated in both cell types, but Tgfbi and Cd274 were upregulated and eventually downregulated in more mature states (Fig. S4 K). This reveals potential differences in the response of AMs and IMs to neoplastic cell engulfment. Together, the transcriptional patterns associated with phagocytosis in IMs suggest that tumor-associated IMs exhibit greater transcriptional plasticity throughout maturation than AMs and do not sustain expression of immune, phagocytic, and metabolic markers that peak during phagocytosis of neoplastic cells.
Phagocytic dendritic cells acquire an mregDC-like phenotype
Although conventional dendritic cells (cDCs) represented a minor proportion of tdTompos cells, understanding phagocytic phenotypes in cDCs is relevant because they are potent antigen presenting cells that play a central role in innate immunity and priming adaptive immune responses (Murphy et al., 2016). A population of tumor-associated cDCs classified as “mregDCs” are transcriptionally distinct from previously described cDC1 and cDC2s and have been characterized by their expression of maturation and immunoregulatory molecules (Maier et al., 2020). From our scRNA-seq dataset, we identified Xcr1-cDC and Ccl22-cDC cell states that correspond to putative cDC1 and mregDCs, respectively (Fig. S3 I). Although we identified cDC2 cells in the lung by flow cytometry (CD11bpos DCs in Fig. 2 D), we did not observe cDC2 cells in our scRNA-seq data likely due to their relatively low frequency. mregDCs had almost double the percentage of tumor associated tdTompos cells compared to cDC1s (Fig. 7 A). mregDCs also contained the highest percentage of tdTompos cells from uninvolved tissue of any myeloid cell type, suggesting a potential migratory capacity of these cells from tumor tissue once they engulf neoplastic cells. Both findings are consistent with previous studies that report increased antigen uptake and high migratory potential in tumor-associated mregDCs and reveal that mregDCs increase these functions after phagocytosing neoplastic cells (Kryczanowsky et al., 2016; Maier et al., 2020; Ness et al., 2021). The PHAT gene signature, defined in monocytes and macrophages, was most highly expressed in tdTompos cDCs, particularly cDC1s, compared to tdTomneg cDCs demonstrating that this signature could be applied beyond monocytes and macrophages (Fig. 7, B and C). Overall, half of the cDC DEGs overlapped with the PHAT-Macrophage signature (Fig. 7 D).
Trajectory inference of cDCs predicted maturation from cDC1s toward mregDCs (Fig. 7, E–G). Analysis of differentially expressed genes across pseudotime not only revealed a concomitant expression shift from DC1 markers (e.g., Ppt1 and Irf8) toward mregDC markers (e.g., Ccl22, Fscn1, and Ccr7), but also an increase in proinflammatory signaling mediators (e.g., Stat4 and Rel) despite higher Cd274 and lower class II MHC (Fig. 7 H). Therefore, tumor-associated tdTompos cDCs in lung tumors from KTai65;SpcCreER mice appear to acquire an mregDC-like phenotype upon engulfment of neoplastic cells in vivo. By tracking neoplastic engulfment in cDCs, our findings reveal an association between phagocytosis and the coordination of migratory, inflammatory, and immunoregulatory functions in cDCs.
Phagocytic TAMs have increased mitochondrial content, oxygen consumption, and expression of OXPHOS-related genes
To understand the functional consequences of the upregulated metabolic gene expression signatures in tumor-associated myeloid cells, we investigated the metabolic features of phagocytic tumor-associated myeloid cells from both the syngeneic and autochthonous KPT:Adeno-SpcCre lung tumor models. Mitochondria are the location of OXPHOS metabolism and their abundance has been shown to affect macrophage phenotypes and functions (Tur et al., 2017). Intracellular staining of mitochondria in myeloid cells (CD11bpos and/or CD11cpos) revealed significantly increased mitochondrial content in tdTompos AM, IM, and DCs compared to tdTomneg cells in both tumor models, suggesting a potential functional role of mitochondrial metabolism (Fig. S5 A and Fig. 8, A and B). Therefore, we next performed extracellular flux analysis of sorted macrophages from autochthonous tumors to measure metabolic activity (Fig. 8, C and D). tdTompos AMs and IMs had significantly increased oxygen consumption rate (OCR) at baseline and at maximal respiration (induced by FCCP treatment), indicative of ATP generation via the electron transport chain in mitochondria (Fig. 8, E and F). However, no significant differences in the levels of spare respiratory capacity were detected despite increased OCR (Fig. S5, B and C). Non-mitochondrial oxygen consumption was significantly higher in tdTompos AMs and IMs compared to tdTomneg AMs and IMs (Fig. 8, E and F). Naive AMs and IMs had modest OCR and tdTomneg AM respiration levels resembled that of the naive AMs. Interestingly, tdTomneg IMs have nearly undetectable levels of cellular respiration. This suggests that IMs that are not activated by phagocytosis of neoplastic cells are in a metabolically quiescent state, unlike AMs which have relatively higher baseline levels of oxygen consumption.
The extracellular acidification rate (ECAR) is indicative of the secretion of lactic acid, a by-product of the glycolytic pathway. ECAR was dramatically elevated in tdTompos AMs and IMs compared to their respective naive or tdTomneg counterparts. A concomitant increase in both OCR and ECAR at baseline suggests that phagocytic tdTompos macrophages have increased glycolysis, in addition to increased OXPHOS (Fig. 8, G and H). This combinatorial increase in OCR and ECAR is indicative of an elevated energetic state in phagocytic tdTompos macrophages. Interestingly, naive IMs have a skewing toward aerobic respiration and lack of glycolytic respiration indicating an energetically quiescent state.
The final step of cellular respiration involves the production of ATP by the mitochondrial ATP synthase, which consists of several subunits, including the protein, ATP5A (Jonckheere et al., 2012). We used ATP5A as a marker of mitochondria and cellular respiration in histological samples from autochthonous lung tumors. To evaluate the spatial localization of macrophages with increased metabolic activity, we performed immunofluorescence staining for ATP5A in autochthonous lung tumors. AMs localized near and within the tumors had increased expression of ATP5A compared to AMs in normal adjacent tissue (Fig. 9 A). To further validate these metabolic features at the single-cell protein level, we quantified the levels of key metabolic enzymes relating to OXPHOS and glycolysis utilizing metabolic flow cytometry (METFLOW) on sorted tdTompos and tdTomneg myeloid cells from autochthonous lung tumors (Ahl et al., 2020). tdTompos IMs AMs and cDCs had significantly increased expression of proteins associated with key metabolic pathways, including GLUT1 (glycolysis) ATP5A (OXPHOS) and CPT1A (fatty acid oxidation [FAO]; Fig. 9, B–E; and Fig. S5 D). Higher expression of CPT1A and ATP5A suggest that mitochondrial lipid transport and respiration is active in phagocytic TAMs and cDCs. Higher expression of GLUT1 is indicative of the sensing and transport of glucose into the cell for glycolysis.
In parallel to METFLOW, we measured the glycolytic activity of phagocytic tumor-associated myeloid cells by treating these cells ex vivo with a fluorescently labeled, non-catalyzable glucose analog (2-NBDG). Glycolytic activity was significantly higher in tdTompos AMs, IMs, and cDCs compared to tdTomneg and naive lung macrophages (Fig. 9, F–I; and Fig. S5 E). These data demonstrate that the transcriptional signatures relating to OXPHOS and glycolytic metabolism produce altered metabolic states in phagocytic TAMs and DCs. Taken together, these data identify tdTompos macrophages as a distinct TAM state in lung tumors and highlight OXPHOS metabolism as a key feature that correlates with immune suppressive phenotypes and functions in lung tumors (Fig. 9 J).
In this study, we used a novel autochthonous genetic lung tumor model that minimizes undesirable, direct labeling of immune cells to identify, quantify, and profile phagocytic myeloid cells at the transcript, protein, and functional levels. Prior studies characterizing TAMs in vivo have not specifically isolated phagocytic cells and most have used either cell transfer models, in which the dominant phagocytic cell is monocyte derived, or autochthonous models that include macrophages that are directly transformed (Casanova-Acebes et al., 2021; Gordon et al., 2017). Our models provide a unique path to define how phagocytosis alters tumor-associated myeloid cells. We find that IMs engulf more neoplastic cells in syngeneic transplant models, while AMs are the dominant phagocytic cell type in autochthonous lung tumors. This could be due to localization of AMs in the alveolar space versus IM localization in the lung parenchyma and the differential source of tumors, either seeding from the blood or initiated from transformed alveolar type II cells.
Within autochthonous lung tumors, we also found key differences between AMs and IMs in the type of response to the engulfment of neoplastic cells in vivo. The inducible response observed in IMs supports prior findings of increased plasticity of monocytes and monocyte-derived macrophages compared to tissue resident macrophages such as AMs, which harbor poised enhancer landscapes shaped by early lung development and/or local microenvironment cues (Lavin et al., 2014). Our scRNA-seq analysis of phagocytic IMs partially correlated with previous studies identifying phagocytic monocyte-derived macrophages in inflamed and cancer tissues (Sanin et al., 2022). The discrepancy of overlap between their signature and ours could be a result of our defining actively phagocytosing cells versus enrichment of genes associated with phagocytosis. Furthermore, we were able to corroborate gene expression of phagocytic TAMs at the transcriptional level with protein detection within AMs and IMs in autochthonous tumors. For example, PD-L1 surface expression was upregulated in tdTompos IMs but was also detected at high levels in AMs irrespective of phagocytic status, consistent with a described tolerogenic role for most tissue resident macrophages (Fig. 2, E and F; and Fig. S2, D and E; Allard et al., 2018). One possible explanation for this is that AMs actively phagocytose lipid-rich lung surfactant at steady-state and in doing so, have high expression of genes involved in tolerogenic pathways at baseline (de Aguiar Vallim et al., 2017).
While our model identifies tdTompos macrophages that have recently internalized neoplastic cells, it cannot explicitly distinguish between macrophages that have degraded the tdTomato fluorophore or have never phagocytosed a neoplastic cell. Trajectory inference of AMs and IMs, allowed us to model the diverse activation paths of macrophages during tissue homeostasis and in the presence of tumors. Notably, the upregulation of Cd274 (which encodes PD-L1) in AM-Trajectory A revealed an association between trajectory-specific neoplastic cell engulfment and immunosuppression (Fig. S4 E). Of the three AM trajectories, the highest percentage of tdTompos AMs were in AM-Trajectory A, suggesting that neoplastic engulfment drives AM maturation toward the cell states enriched in tissue remodeling genes, culminating in a Chil3high AM state (AM-4; Fig. 5 I and Fig. S4 A). Our results suggest that the AM-4 cell state, despite having a lower fraction of tdTompos cells than the AM-2 state, derives from actively phagocytic AMs and represents an imprinted consequence in AMs following phagocytosis of neoplastic cells. AM-4 cells were also detected in AM-Trajectory C, although AM-4 cells in AM-Trajectory A had higher immune-suppressive gene expression. The high fraction of AMs in naive and uninvolved tissues in AM-Trajectory C suggests that certain homeostatic processes may polarize AM maturation toward phenotypes similar to those that arise after neoplastic cell engulfment.
In addition to defining cell type–specific phagocytosis-driven TAM phenotypes, we identified a common 123-gene phagocytic signature. Importantly, this PHAT signature was enriched in human lung cancer tumor-associated myeloid cells and correlates with a lower probability of survival in patients with lung adenocarcinoma (Fig. 4 I). In human tumors, direct identification of myeloid cells that have recently phagocytosed neoplastic cells is not feasible. However, one study used morphological features to characterize putative phagocytic large TAMs (L-TAM) and compared them to non-phagocytic small TAMs (s-TAM; Donadon et al., 2020). The enrichment of lipid metabolic pathways in L-TAMs is consistent with our findings in tdTompos TAMs. Interestingly, our PHAT signature better represents human TAM gene expression than traditional M1 or M2 phenotypes suggesting that gene signatures associated with macrophage functions, such as the phagocytosis of neoplastic cells, more closely identifies TAMs in clinical samples.
Phagocytic tdTompos cDCs were identified both in tumors and adjacent tissues (Fig. 7), and Ccl22 mregDCs were enriched in phagocytic cells and represented a more differentiated state compared to Xcr1 cDC1s, consistent with other studies (Maier et al., 2020). These cells represented the majority of tdTompos cells from non-tumor tissue, potentially due to their enhanced migratory capacity. Additionally, tdTompos cDCs have significant enrichment of the PHAT Signature, in which OXPHOS metabolism is a key feature. Functional assays are consistent with phagocytic DCs being more metabolically active, with increased OXPHOS. These DCs have higher mitochondrial content and ATP5A and CPT1A expression, as well as increased glycolysis, with higher glucose uptake and GLUT1 expression. A link has been proposed between lipid-fueled FAO and cancer-associated DCs with tolerogenic functions in tumors (Herber et al., 2010; Zhao et al., 2018). Thus, we provide further evidence that engulfment of neoplastic cell may promote a tolerogenic mregDC program which is sustained by metabolic reprogramming.
OXPHOS metabolism is associated with immune suppressive phenotypes of macrophages in vitro (Jin et al., 2014; Viola et al., 2019). Macrophages can acquire lipids from the engulfment of dying cells (Tabas and Bornfeldt, 2016; A-Gonzalez et al., 2017); these lipids are transferred to the mitochondria to replenish the TCA cycle with acetyl-CoA for the generation of NADH and FADH2 to support the electron transport chain (ETC; Mehta et al., 2017). Our functional metabolic assays demonstrated that phagocytic tumor-associated AMs and IMs undergo metabolic reprogramming toward increased OXPHOS metabolism. Earlier reports demonstrated that type-2 cytokines in vitro or in vivo can maintain an intact TCA cycle coupled with mitochondrial OXPHOS (Jha et al., 2015). Recent analyses of macrophage metabolism describe how lipid accumulation occurs within macrophages in chronic inflammatory environments due to defective ETC when tissue oxygen levels are low, ultimately leading to an inefficient ETC and subsequent inhibition of OXPHOS metabolism (Eming et al., 2021). However, we find that the OXPHOS metabolism is intact in phagocytic TAMs (Fig. 8, C–H). Thus, phagocytosed material, including lipids derived from neoplastic cells, can be catabolized via FAO and the TCA cycle in TAMs from autochthonous tumors. Although we connect OXPHOS metabolism with phagocytosis in TAMs, more work will be needed to elucidate the requirement of this energetic shift to sustain TAM phenotypes and functions in vivo and, most importantly, how these phenotypes change with therapeutic intervention.
The connection of TAM and DC phagocytosis with metabolic changes toward increased OXPHOS and an immune suppressive/tumor-promoting phenotype links previous studies defining metabolic phenotypes of phagocytic macrophages in vitro with characterization of TAM metabolism in vivo (Lantz et al., 2020; Li et al., 2022). Similarly, disease associated microglia (DAM) and lipid-associated macrophages (LAM) have been shown to utilize lipid metabolic pathways and OXPHOS and share significant gene expression overlap with each other (Jaitin et al., 2019; Keren-Shaul et al., 2017). Interestingly, we found significant transcriptional similarity between PHAT macrophages and lipid-rich Trem2pos LAMs (Fig. 4 I). This finding underscores a recurrent theme among macrophages in inflammatory microenvironments, in which the magnitude of lipid metabolism and oxidation may direct macrophage fate.
Lipid metabolism has been reported to be essential for macrophage alternative activation of macrophages and efficient efferocytosis (Huang et al., 2014). In line with the downstream catabolism of lipids potentially derived from recently phagocytosed cells, key components of the ETC are significantly upregulated in tdTompos TAMs including ATP5A, a subunit of the complex V of the ETC, and CPT1A, a critical transporter of fatty acids from the cytosol into the mitochondria (Fig. 9, B–D). Using immunofluorescence, we show that ATP5A is detected at higher levels in TAMs that are in close proximity to neoplastic cells in tumor tissue, showing in situ the increased metabolic state of macrophages positioned to phagocytose cancer cells. These data suggest that myeloid cell phagocytosis of neoplastic cells results in dynamic metabolic changes that could inform the design of novel therapies targeting TAM functions and induce macrophage repolarization.
Materials and methods
Autochthonous mouse models were performed at Stanford animal facilities and used KrasLSL-G12D;Trp53flox/flox;Rosa26LSL-tdTomato (KPT) mice and KrasFSF-G12D;Rosa26FSF-LSL-tdTomato(ai65);SftpcCreER (KTai65;SpcCreER) mice. KrasLSL-G12D (Jackson et al., 2005; Jackson et al., 2001; Tuveson et al., 2004), Trp53flox/flox (Jax# 008462; Marino et al., 2000; Jonkers et al., 2001), Rosa26LSL-tdTomato (Jax# 007905; Madisen et al., 2010), KrasFSF-G12D (Jax# 023590; Young et al., 2011), Rosa26FSF-LSL-tdTomato(ai65) (Jax# 021875; Madisen et al., 2015), and SftpcCreER (Jax# 028054; Rock et al., 2011) have been previously characterized. Intravenous cell transfer experiments were conducted at Amgen animal facilities and used 8-12-wk-old F1 (first filial generation) hybrid C57BL/6J and 129S1/Svlmj mice (JAX Strain#: 101043; B6129SF1/J) using a primary lung tumor cell line derived from a KPT autochthonous lung tumor. All animal experimental protocols performed at Amgen were approved by the IACUC of Amgen and were conducted in accordance with the guidelines set by the Association for Assessment and Accreditation of Laboratory Animal Care. The Stanford Institute of Medicine Animal Care and Use Committee approved all animal studies and procedures performed at Stanford University.
Lenti-PGK-Cre and Ad5-SPCcre generation
Lentiviral plasmids were co-transfected with packaging vectors (delta8.2 and VSV-G) into 293T cells using PEI. The supernatant was collected at 36 and 48 h after transfection, ultracentrifuged at 25,000 rpm for 90 min, and resuspended in PBS. Tumors were initiated by intratracheal infection of mice with lentiviral or adenoviral vectors expressing either inverted-FLPo or Cre recombinase (as indicated for KTai65;SpcCreER or KPT-AdenoSPCcre tumor models). Ad5-SPCcre adenoviral vectors were purchased from University of Iowa Carver College of Medicine Viral Vector Core, deposited by Dr. Anton Berns from Netherlands Cancer Institute (VVC-Berns-1168). For Ad5-SPCcre transductions, mice were given an intratracheal dose of 1 × 109 TFU in 60 μL of PBS.
Lentiviral backbone with inverted FLPo recombinase (Lenti-InFlpo), was generated through modification of pLL3.3 lentiviral vector by addition of PGK promoter and inverted FLPo recombinase surrounded by pairs of heterotypic lox sites (LoxP and Lox 2272). These additional fragments were generated by Integrated DNA Technologies and assembled through NEB HiFi DNA assembly system (E2621). In this design, there are two different lox sites (LoxP and Lox 2272) on each side of FLPo Recombinase. Identical lox sites (one on each side of FLPo recombinase) are facing each other and they invert the FLPo recombinase as a result of Cre mediated recombination. To prevent expression of FLPo expression in the absence of Cre recombination even when the Lentiviral vector integrates next to an endogenous promoter, we included an HSV TK poly A between the 3′ LoxP and Lox 2272 sites (unpublished data).
Tumor dissociation and cell sorting
Primary tumors were dissected and dissociated as previously published (Yu et al., 2016). Briefly, (n = 3 tumor-bearing mice and n = 3 tumor-naive mice) using an enzymatic digestion mixture of 1.5 mg/ml of Collagenase A (Sigma-Aldrich: COLLA-RO) and 200 units/ml DNaseI in DPBS (with Ca2+/Mg2+) plus 5% heat-inactivated fetal bovine serum (HIFBS) and 10 mM HEPES. Lungs were inflated with digestion mixture and incubated at 37°C for 30 min shaking at 120 rpm. The lungs were mashed and sequentially strained, washed with media (RPMI, no Glut containing 10% HIFBS), and centrifuged at 450 ×g through 100-, 70-, and 50-μm cell strainers. The resulting cell suspensions were treated with ACK RBC lysis solution (A1049201; Gibco) for 5 min at room temperature and washed/centrifuged with media. For cell sorting, antibody labeled cells were sorted based on detection of indicated antibody fluorescence into 1.7-ml Eppendorf tubes containing media and placed on ice for downstream applications.
Flow cytometry and staining
After cell counting, at least 1 × 106 cells per sample were stained with Near-IR Fixable LIVEDEAD dye (L10119; Life Technologies) according to the manufacturer’s instructions. Cells were then incubated in FACS buffer (554657; BD) containing 1% TruStain FcX anti-mouse CD16/32 antibody (101320; Biolegend) for 15 min on ice. Subsequently, cells were washed with FACS buffer, centrifuged, and stained with immunophenotyping antibodies (Anti-CD11b, clone: M1/70; Anti-CD45, clone: 30-F11; Anti-CD11c, clone: N418; Anti-Ly6G, clone: 1A8; Anti-I-A/I-E, clone: M5/114.15.2; Anti-CD64, clone: X54-5/7.1; Anti-CD24a, clone: M1/69; Anti-PD-L1, clone: 10F.9G2; Anti-CD80, clone: 16-10A1; Anti-MERTK, clone: 2B10C42; Anti-ITGAV, clone: RMV-7; Anti-IL-4Ra, clone: mIL4R-M1; Anti-TREM2, clone: FAB17291A; Anti-AXL, clone: MAXL8DS; Anti-EpCAM, clone: G8.8) for 25 min on ice. After staining, cells were washed with FACS buffer and either analyzed on the flow cytometer or cell sorter directly or fixed with DPBS containing 2% paraformaldehyde for 10 min at room temperature, protected from light. Cells were then washed with FACS buffer and stored at 4°C until analyzed on the flow cytometer. For Metabolic Flow cytometry (METFLOW), staining procedures were described elsewhere (Ahl et al., 2020). Briefly, cell suspensions were stained with cell surface antibodies as detailed above then cells were fixed and permeabilized using eBioscience Foxp3/Transcription Factor Staining Buffer Set (Cat# 00-5523-00; Invitrogen) according to manufacturer’s instructions. Cells were then washed with cold DPBS and centrifuged at 450 ×g, 5 min, three times. Subsequently, cells were stained with intracellular antibodies (ATP5A, clone: EPR13030(B); CPT1A, clone: 8F6AE9; GLUT1, clone: EPR3915) in permeabilization buffer for 1 h at room temperature. Then, cells were washed once with permeabilization buffer and then once with cold DPBS. Samples were acquired on a X-30 FACSymphony (BD) with FACS Diva Version (BD, Version 8.0.1) software. Analysis was performed using FlowJo (Treestar, version 10.5.2). TdTomato-positive and -negative population percentages only displayed cell populations above a 5% threshold of CD11bpos and/or CD11cpos cells.
Nanostring gene expression analyses
RNA was purified from ∼200,000 FACS sorted tdTompos or tdTomneg lung macrophages (AMs: CD45pos Ly6Gneg MHCIIpos SSC-Ahigh CD64pos CD24neg CD11blow; IMs: CD45pos Ly6Gneg MHCIIpos SSC-Ahigh CD64pos CD24neg CD11bhigh) using the Qiagen Rneasy kit, per the manufacturer's recommendations. RNA quantification of selected genes was performed using the nCounter myeloid innate immunity panel for mouse consisting of 770 target genes from 19 pathways with housekeeping and control genes (NanoString Technologies). Pathway scores and P values were calculated from quality-controlled, normalized raw expression data using the nSolver software (v4.0) with the nCounter Advanced Analysis add-on (v2.0.134).
Phagocytosis in vitro co-culture assay
BMDMs were derived per previous well-established protocol (Weischenfeldt and Porse, 2008). Briefly, whole bone marrow was flushed from hindlimb bones of C57BL/6 mice (female, between 6–8 wk of age) using PBS. Red blood cells were lysed (R7757-100ML; Sigma-Aldrich) and remaining bone marrow cells were cultured in BMDM media (DMEM plus 10% HIFBS and 1% penicillin/streptomycin media that contains 100 ng/ml recombinant murine M-CSF; Peprotech: 315-02) for 6 d with a media change on day 3. Fully differentiated BMDM were detached, counted, and plated into 12-well plates at 2 million cells per well (n = 5) and allowed to settle for 2 h. KPT cancer cells (cell clone: SP110P) were counted and resuspended in BMDM media to overlay onto BMDMs at 1:1 ratio. The co-culture was incubated at 37°C overnight (16–18 h), at which point cells were detached and subject to flow cytometry–based analyses to evaluate the status of their surface marker expression and intracellular metabolic activity (see Flow cytometry and staining; and Metabolic characterization using extracellular flux analysis and cellular dyes).
Metabolic characterization using extracellular flux analysis and cellular dyes
ECAR (mpH/min) and OCR (pmol/min) were measured using the XFe96 extracellular flux analyzer (Agilient Technologies). Glycolytic function correlated with ECAR and mitochondrial respiration correlated with OCR. 200,000 cells were plated per well in complete media (RPMI 10% HIFBS 1× pen/strep/glut) in a 96-well plate after cell sorting and PBS wash. The assay was performed using XF Assay Modified Media with L-glutamine (2 mM; 103579-100; Agilient Technologies), sodium pyruvate (1 mM; 103578-100; Agilient Technologies) with 11 mM glucose (103577-100; Agilient Technologies). Mitochondrial respiration was measured using the mitochondrial stress test kit (103015-100; Agilient Technologies), by sequentially adding oligomycin (2 µM), carbonyl cyanide-4 (trifluoromethoxy) phenylhydrazone (1.5 µM), and rotenone and antimycin A (1 µM). For measurement of glycolysis using a dye as described elsewhere (Bala et al., 2021), 200,000 isolated tumor cells were resuspended in complete media containing 2-(N-(7-Nitrobenz-2-oxa-1,3-diazol-4-yl)Amino)-2-Deoxyglucose (2-NBDG) at 60 µM for 1 h at 37°C. Cells were then analyzed for via flow cytometry for fluorescence intensity on the FITC detector. For mitochondrial staining, MitoTracker Deep Red FM (M22426; Invitrogen) was used per the manufacturer’s recommendations and analyzed for MFI via flow cytometry.
Histology and immunofluorescence
Tumor tissue samples were fixed in 10% neutral buffered formalin for 12–24 h at room temperature. Fixed samples were embedded in paraffin blocks, and 5-μm sections were cut. H&E and staining procedures were performed by Pathology core at Stanford University. Fluorescent staining was performed with anti-RFP (600-401-379; 1:200, Rockland Immunochemicals) with secondary anti-rabbit Alexa Flour 647 Ab, ATP5A (ab196198, 1:200; Abcam) with secondary anti-rabbit Cy3 Ab. Pictures were taken on Keyence BZ-X710 using ×4 PlanFluor NA 0.13 PhL, ×20 PlanApo NA 0.75, and ×40 PlanFlour NA 0.60 Ph2 objective lens.
Sample preparation and FACS for scRNA-seq
Single cell suspensions of lung tumor cells were prepared as described above. Sample hashing was performed using TotalSeq-A (BioLegend) reagents. Specifically, single cell suspensions were stained with oligo-tagged antibodies specific for B2M and CD45 at a concentration of 1:100 for 25 min on ice and subsequently washed and stained for surface markers as described above. Prior to sorting on the FACS machine, adjacent normal and tumor samples of the same biological replicate were combined equally and myeloid cells (Singlet Live CD45+CD11b+CD11c+) and epithelial/neoplastic cells (Singlet Live CD45-EpCAM+) cells were sorted based on their tdTomato fluorescence (tdTompos or tdTomneg).
scRNA-seq library preparation and sequencing
Sorted single-cell suspensions were pooled into scRNA-seq sample groups according to the sample hashtag labeling strategy to identify cells by animal and tissue type of origin. Samples were centrifuged for 400 ×g × 7 min and resuspended in RPMI + 10% FBS at a concentration of 400–1,000 cells/μl for samples containing over 50,000 cells and at the resuspension concentration (200–300 cells/μl) for samples with <50,000 cells. Cell encapsulation and transcriptomic library preparation were performed using the 10X Genomics Chromium 3′ V3 chemistry according to manufacturer’s instructions (protocol CG000185 RevD). Targeted cell recovery was 10,000 cells for samples with over 50,000 cells sorted, and 46.6 μl of cell resuspension was loaded for samples with <50,000 cells sorted.
After encapsulation of cells and hydrogel beads into Gel Bead-In Emulsions and reverse transcription, 11–12 cDNA amplification cycles were performed with standard 10X Genomics gene expression cDNA primers with 2 nM HTO additive primer to amplify native transcripts and hashtag oligo (HTO) barcodes, respectively. SPRI clean-up of cDNA amplification products was performed to separate HTO products (in the supernatant) from gene expression (GEX) products (in the pellet). GEX products were then subject to fragmentation, end-repair, A-tailing, adapter ligation, and an additional round of 11-15 PCR cycles to attach library-specific 8-bp i7 indexes. After double-sided SPRI bead purification, this resulted in the generation of ∼500 base pair libraries for next-generation sequencing. Libraries were quantified using Agilent TapeStation 4200 D5000 high-sensitivity screen tapes.
HTO libraries were constructed using the BioLegend TotalSeq-A HTO protocol (https://www.biolegend.com/en-us/protocols/totalseq-a-antibodies-and-cell-hashing-with-10x-single-cell-3-reagent-kit-v3-3-1-protocol). SPRI-purified HTO fractions were amplified using the SI-PCR primer and TruSeq D70×_LONG 8-bp library-indexing primers for 13–15 cycles. Products were purified again using SPRI purification to generate final libraries and quantified using Agilent TapeStation 4200 D1000 high-sensitivity screen tapes.
To sequence the libraries, each GEX and HTO library was diluted to 5 nM and pooled at a volume ratio of 4:1 for GEX:HTO, resulting in four times more GEX reads than HTO reads. Libraries were then sequenced on two lanes of an S4 version v1.5 flow cell on an Illumina NovaSeq6000 instrument to achieve an average sequencing depth of 50,000 reads per cell. The sequencing parameters used were 28:8:0:91 cycles (Read1:i7 index:i5 index:Read2).
scRNA-seq alignment, preprocessing, and normalization
Raw sequencing files were demultiplexed by GEX and HTO library i7 barcodes to generate FASTQ files for each GEX and HTO library. Reads from each library were aligned using cellranger count (CellRanger version 4.0.0 from 10X Genomics). More than 90% GEX reads aligned to the Ensembl 93 mm10–3.0.0 mouse reference genome containing the tdTomato gene (https://www.ensembl.org), and ∼65% of the total reads aligned to the transcriptome. HTO reads were aligned to a feature reference (using the argument --feature-ref) containing the TotalSeq-A barcodes used, and ∼60% of reads contained a valid HTO barcode. The unfiltered UMI count matrices generated by CellRanger were used for pre-processing of data.
Quality control was performed using a multi-step procedure. For preliminary removal of cell barcodes containing only ambient mRNA or poor-quality transcriptomes, outlier barcodes that were more than three median average deviations below the median read count were removed. Cell barcodes with <250 mRNA features, <500 UMI, sequencing saturation <50%, and percent of mtRNA reads >20% were also considered to be poor-quality transcriptomes or dying cells and removed (Fig. S3 C).
HTO counts were normalized across cells using the centered log-ratio transformation and demultiplexed using the HTODemux() function with default settings in the Seurat v4.0.0 R package (Hao et al., 2021; Stoeckius et al., 2018; Fig. S3 D). A single HTO barcode could be assigned to 72.3% of cells, and HTO labels were later used to filter cell multiplets. Alignment of sequencing reads generated 2.7 billion total genome-mapped reads across and 951 million hashtag barcode-mapped reads (89.7% associated with gene expression cell barcodes) from 10 gene expression and antibody hashtag libraries.
For preliminary removal of multiplets by read count and number of features, barcodes with a read cutoff of 10 median average deviations above the median read number and >6,000 mRNA features were removed. Next, cell barcodes with two or more HTO barcodes identified by the HTO classifier were removed. Finally, the R package DoubletFinder (McGinnis et al., 2019) was applied to each 10X sample library to remove doublets not removed by either method. Briefly, artificial doublets containing averaged expression values from two cell barcodes were simulated and integrated into the dataset. Normalization, scaling, and dimension reduction using principal-components analysis were used to embed cell barcodes and artificial doublets. The proportion of artificial nearest neighbors (pANN) were then computed using the k nearest neighbors for every cell barcode in principle-component space, and the barcodes with top pANN values were classified as doublets. Overall, we observed a 77.7% concordance in predicted singlets between HTO labeling and DoubletFinder, with 15.6% of doublets detected by either HTODemux() or DoubletFinder (Fig. S3, E and F).
Cell cycle stage was calculated for each cell using the machine learning approach previously described (Scialdone et al., 2015). This approach trains a classifier on pairs of genes that change expression directionality across different cell cycle phases. When applied to our data, each cell’s phase could be predicted by observing the calculated gene expression difference for either G1, G2M, or S phases (Fig. S3 J).
All 10 libraries were normalized using the procedure outlined previoulsy (Ellwanger et al., 2021). In short, scaling factors for each library were first calculated. Then, similar cells were binned into pools containing at least 10% of the total number of cells per cluster and a minimum average gene expression of 1 using the quickCluster function in the R package scran (Lun et al., 2016). Pool-based size factors were calculated and deconvolved to yield cell-based size factors. The cell-based size factors were normalized between libraries based on the ratio of average UMI counts between libraries using the R package batchelor. The resulting UMI count matrix was log2-transformed and used for downstream analyses.
Macrophages were distinguished from monocytes and dendritic cells based on established combinatorial markers (e.g., ItgaxlowNr4a1highLy6c2high for monocytes and Flt3highLyz2low for dendritic cells). Alveolar macrophages could be separated from interstitial macrophages based on marker genes such as Ear2, Siglecf, and Lpl.
Variance modeling, spectral embedding, and cell projection
To quantitatively find the subset of genes that have significantly high variance between cells, we used the approach implemented in (Ellwanger et al., 2021). Briefly, the variance of each gene was fitted as a function of mean expression. Genes were tested for significance using the R package scran by modeling the residuals with an F-distribution, and genes with a P value <0.05 were considered significant. To retain tissue and FACS-specific variance, the union of all significant variable genes in all sample libraries was used.
Data dimensionality was reduced using principal component analysis (PCA) on the highly variable genes (R package irlba). Library- and animal-specific variations in cellular expression were detected and determined to be the main confounding variables to mask tumor and tdTomato-based variation using the plotExplanatoryVariables function in the R package scater (version 1.0.4). Therefore, Harmony was used to integrate PCA embeddings across library and animal (Korsunsky et al., 2019). After integration, we aimed to increase the signal-to-noise ratio of the spectral embedding. A factor analysis using a scree plot determined the optimal number of components to use for computing a two-dimensional uniform manifold approximation and projection of single cells.
Clustering and cell classification
To cluster cells for cell type and subset identification, a shared nearest neighbor graph was first constructed using the top 40 integrated principal components. The Louvain algorithm in the R package Seurat was used to cluster cells. The optimal resolution parameter for unsupervised clustering was empirically determined using a systematic assessment of clustering robustness provided in the R package chooseR (Patterson-Cross et al., 2021). Briefly, cells were subsampled and iteratively clustered across 20 different resolution parameters, ranging from 0.2 to 12. For each resolution, 100 bootstrap replicates were performed. Next, co-clustering frequencies for each cell were summed at each resolution and used to calculate mean per-cluster silhouette scores. From the distribution of Silhouette scores, a resolution of 1.6 was selected. Clustering resulted in the detection of 33 cell clusters, which could be combined into 16 major cell types.
Differential expression using the Wilcoxon rank sum test with a cutoff of FDR-adjusted P value <0.05 and log2FC >0.25 (FindAllMarkers function in Seurat) was utilized to identify marker genes for each cell. Cells were assigned cell types by matching cluster marker genes to previously reported cell marker genes and by expert annotation.
To manually confirm that minor populations were not missed at this clustering resolution, cells were clustered at resolutions of 1.8 and 2.0 and FindAllMarkers() was applied to corroborate that additional clusters identified did not express unique marker genes.
Differential expression and analysis of functional scores
To calculate differential gene expression between any two populations of cells, we applied MAST (Finak et al., 2015) to evaluate the significance of each gene, using the Bonferroni correction for multiple hypothesis correction. Genes with an adjusted P value < 0.05 and log2-fold change (log2FC) > |0.25| were considered to be differentially expressed.
To score cells using gene sets pertaining to shared cell function, we used the approach reported in Tirosh et al. (2016). We used the log-normalized expression of signature genes to calculate functional scores using the AddModuleScore function in Seurat. Differential expression of functional scores was tested using the Wilcoxon rank sum test in R.
For determination of alveolar macrophage, monocyte-to-interstitial macrophage, and classical dendritic cell development and differentiation trajectories, relevant classified cell types were first subsetted. Variable features were selected using the F-distribution (see Variance modeling). Dimensionality reduction and cell embedding were calculated using diffusion maps (Coifman and Lafon, 2006), which resolves both linear and non-linear substructures across cell states in a deterministic manner. Diffusion components were calculated from a cell-by-cell expression matrix of the distance function representing each cell’s relative dissimilarity was first generated. Then diffusion components were calculated using multiple local sigmas in the diffusion kernel as performed in Haghverdi et al. (2015). Correction of batch effects across libraries and animals was performed using Harmony. The number of components for lower-dimensional embedding was selected using a scree plot analysis, and the resulting manifold was visualized to approximate the relationship between cells from the same lineage. Pseudotime values and relative ordering of cell subsets were then assigned to cells using the getLineages function in the R package slingshot (Street et al., 2018).
RNA velocity analysis
RNA velocity was determined to understand the directionality and differentiation potential for cells within lineages. Spliced and unspliced mRNA UMI counts were tallied for the 10X scRNA-seq CellRanger output from each sample library using the python package velocyto (version 0.6). Spliced and unspliced count matrices were then merged across libraries and normalized by library size factors (see Normalization section under scRNA-seq pre-processing). Genes were filtered to select genes expressed in only 3% or more cells for each lineage in both count matrices. The python package scvelo (version 0.2.3) was used to derive the velocity kinetics. First, principal components analysis was performed on the filtered data, and a k-nearest neighbor graph (k = 30) was built using the top 30 principal components to connect similar cells. For each cell, first- and second-order means and uncentered gene variances (moments) were calculated using scvelo.pp.moments, and these moments were passed to scvelo.tl.recover_dynamics to derive splicing kinetics of selected genes. Dynamical modeling of velocity was calculated using scvelo.tl.velocity (mode = “dynamical”), which was then used in scvelo.tl.velocity_graph to construct a velocity graph based the cosine similarities between potential state transitions and velocities. Last of all, velocity streams were visualized from the velocity graph using scvelo.pl.velocity_embedding_stream and embedded onto DiffusionMap cell projections.
TCGA analysis of clinical outcomes
Clinical phenotypes and expression data from the TCGA Clinical Data Resource (Liu et al., 2018) were downloaded using the Diseaseland Explorer (from Qiagen OmicSoft Studio V11.5). To test for the correlation between the PHAT-related gene expression and overall survival (OS) in lung adenocarcinoma (LUAD) and head and neck squamous cell carcinoma (HNSC), a signature score was calculated using the mean mRNA expression values for PHAT genes shared between monocytes, macrophages, and dendritic cells. The mean expression was then divided by the expression of PTPRC to adjust for CD45+ cell levels for each sample. Samples were then grouped into PHAThigh or PHATlow categories based on expression above the 60th (PHAThigh) or below the 40th (PHATlow) percentiles. The Cox proportional hazards model was used in GraphPad Prism 9.3.1 to test for significance between the selected genes and OS after correcting for the effect of gender, race, age, and pathologic stage.
Meta-analyses with published scRNA-seq datasets
For the Sanin et al. (2022) datasets, cluster-defining genes were extracted from supplemental data files 1 and 2. Ribosomal proteins were excluded from the author’s gene lists and were subsequently removed from the interstitial macrophage signature gene list. For the Jaitin et al. (2019) datasets, differentially expressed gene lists provided in dataset S5 and dataset S6 were filtered to remove genes with P value≥0.05. Gene enrichment in the Mac3/LAMs or Mac1/Macs clusters was determined by further filtering by log2 fold change ≥0.25 or log2 fold change ≤0.25, respectively. For both Sanin et al. (2022) and Jaitin et al. (2019) datasets, overlap of the interstitial macrophage signature genes and the cluster-defining genes was determined from the Jaccard index using GeneOverlap (version 1.34.0).
For the Lantz et al. (2020) dataset, normalized and pre-processed single cell data were imported from BioTuring as a Seurat object (Hao et al., 2021). Cells were scored for expression of the PHAT signature using UCell (version 2.2; Andreatta and Carmona, 2021). Statistical significance of signatures scores was determined using the unpaired Wilcoxon test.
PHAT gene signature analysis of human TAMs
Tumor samples from patients with non-small cell lung cancer, with a primary diagnosis of lung adenocarcinoma or lung squamous carcinoma were collected and analyzed as described in Combes et al. (2022). Briefly, myeloid cells were sorted from single-cell suspension of freshly harvested tumors or normal adjacent tissue using the expression of the following protein markers: negative for viability dye; CD45+ CD3− CD19− CD56− HLADR+ CD11b+/−; CD90+ CD44+. Cells were sorted into lysis buffer and mRNA was isolated. Libraries were created and sequencing was performed as described. Median expression of genes comprising PHAT signature; M2 or M1 CIBERSORT signature (Newman et al., 2015) was calculated for each myeloid sorted sample. Statistical analysis used was a one-tailed t test with two-sample unequal variance, tumor vs. normal.
Online supplemental material
Fig. S1 shows the frequency and phenotype of tdTompos macrophages in syngeneic, autochthonous, and in vitro tumor models. Fig. S2 shows the frequency and phenotype of tdTompos macrophages in Adeno-SPCcre autochthonous lung tumors. Fig. S3 shows a summary of quality control metrics for scRNA-seq, including tissue dissections, cell sorting gating strategy, library sample preparation, doublet discrimination, heat map of cell clusters and cell cycle scores. Fig. S4 shows the identification, frequency, and gene expression trajectory of AM and IM subsets in naive lung or dissected autochthonous tumors at the single cell level. Fig. S5 shows the mitochondrial levels in tdTompos or tdTomneg immune cells from AdenoSPcre or syngeneic tumors and the metabolic phenotypes of AMs, IMs, and DCs from AdenoSPCcre autochthonous tumors. Table S1 shows summary of differentially expressed genes between tumor tdTom-pos vs. tumor tdTom-neg cells (Fig. 4 A). Table S2 shows human orthologs of PHAT signature genes. Table S3 shows differential expression results for tumor tdTom-pos alveolar macrophages vs. tumor tdTom-neg alveolar macrophages. Table S4 shows differential expression results for tumor tdTom-pos interstitial macrophages vs. tumor tdTom-neg interstitial macrophages. Table S5 shows differential expression results for tumor tdTom-pos monocytes vs. tumor tdTom-neg monocytes.
Data and availability
scRNA-seq metadata and sequencing files can be found in the Gene Expression Omnibus (accession ID: GSE217368).
We thank Stanford Veterinary Animal Care Staff, the Amgen In Vivo Pharmacology department, and the Amgen Flow Cytometry Services department. We thank Tyler Jacks and Tushar Desai for providing mouse lines. We also thank Sarah Pierce for generating the KPT cell lines in the Winslow lab.
M. Yousefi was supported by a National Institutes of Health (NIH) Ruth L. Kirschstein National Research Service Award (NIH F32-CA236311). C.W. Murray was supported by the NSF Graduate Research Fellowship Program and an Anne T. and Robert M. Bass Stanford Graduate Fellowship. E.L. Ashkin was supported by Howard Hughes Medical Institute Gilliam Fellowship for Advanced Study (GT14928). This work was supported in part by NIH R01-CA231253 (to M.M. Winslow) and NIH R01-CA234349 (to M.M. Winslow).
Author contributions: M.A. Gonzalez and D.R. Lu contributed equally to the completion of the manuscript. M.A. Gonzalez and D.R. Lu performed data curation, formal analysis, investigation, methodology, visualization, and writing. M.A. Gonzalez performed all flow cytometry analyses. H. Zhou, V. Arias, and D.R. Lu performed sample preparation for Nanostring and scRNA-seq cell capture and downstream analysis. M Yousefi contributed to methodology and development of autochthonous lung tumor mouse models, data generation and writing. A Kroll contributed analysis to compare signatures across datasets. C.W. Murray, E.L. Ashkin, M.K. Tsai provided technical assistance and editing. C.G. Briseño performed scRNA-seq analysis of DC phenotypes and contributed to editing. J.E.V. Watson performed data curation and analysis for Immunoprofiler dataset and contributed to writing. S. Novitskiy performed methodology and visualization of immunofluorescence and contributed to writing. A.P. Stapper performed QC and data curation for the Immunoprofiler dataset. C.M. Li provided project administration and resources. K.V. Tarbell and M.M. Winslow contributed through study conceptualization, methodology, writing, project administration, and supervision.
M.A. Gonzalez and D.R. Lu equally contributed to this paper.
Disclosures: M.A. Gonzalez is a current employee of Gilead Sciences. D.R. Lu reported currently being an employee of Amgen, Inc. M. Yousefi reported other from Pionyr Immunotherapeutics outside the submitted work. C.G. Briseño reported other from Amgen, Inc. during the conduct of the study; other from Amgen, Inc. outside the submitted work; and being an employee of Amgen, Inc. J.E.V. Watson reported being an employee of Amgen, Inc. and holding Amgen, Inc. stock. V. Arias reported currently being employed at Amgen, Inc. C.M. Li reported non-financial support from Amgen, Inc. during the conduct of the study; and personal fees from Amgen, Inc. outside the submitted work. M.M. Winslow reported grants from NIH during the conduct of the study; and is a cofounder, advisor, and has equity in D2G Oncology, Inc. K.V. Tarbell reported personal fees from Amgen, Inc. during the conduct of the study; personal fees from Amgen, Inc. outside the submitted work; and being an employee of Amgen, Inc. In addition, Amgen, Inc. financially supported the research included in this manuscript. No other disclosures were reported.