Single cell transcriptomics of human and mouse lung cancers reveals conserved myeloid populations across individuals and species

scRNA-Seq profiling maps immune cell gene expression in human and mouse lung cancer

We set out to define TIM cell sub-populations in human NSCLC and in a mouse lung adenocarcinoma (LA) model, and to establish similarities and discrepancies between the two organisms (Fig. 1A). We first obtained fresh tumor tissue from seven human patients undergoing resection of one of two NSCLC types: LA (5 of 7) and squamous cell carcinoma (SCC; 2 of 7). The samples were from patients who had not yet received cancer treatment in all cases but one (Table S1). The patients ranged from 61 to 83 years old; the tumors were of diverse tumor stages and split between female (n=4) and male (n=3) subjects. The biopsies were dissociated, washed, and rapidly processed using the inDrop scRNA-Seq platform (Zilionis et al., 2017). After filtering scRNA-Seq data to exclude putative cell doublets (Wolock et al., 2018) and stressed or dead cells, a total of 40,362 cell transcriptomes were retained for analysis across the seven patients. These data were visualized for exploration using SPRING, a two-dimensional force-layout embedding (Weinreb et al., 2018) (Fig. 1B). Of all cells, 15% of them (n=5,912) were classified as fibroblasts, endothelial, erythroid and epithelial cell types. The population structure of these non-immune cells is analyzed in Fig. S1A–E and data on these cells are available for further exploration (see Key Resource Table). The remainder of the study focused on immune cells (Fig. 1C; n=34,450).

A matched dataset was obtained by scRNA-Seq profiling of immune cells from healthy and tumor-bearing mice (Fig. 1D). To induce LA tumors, mice were injected with the KP1.9 tumor cell line (Pfirschke et al., 2016). We isolated and sequenced a total of ~17,000 CD45+ single cells from the lungs of healthy (i.e. tumor-free) and tumor-bearing mice (two mice per condition). After excluding cells of high mitochondrial load (an indicator of cell death) and putative doublets, 15,939 single cell transcriptomes were retained for analysis, and visualized using SPRING (Fig. 1D, Figs. S1F, G).

Both human and mouse immune cells partitioned into several major cell clusters, some with complex internal structure, as well as minor clusters representing <1% of the total cell population. To define the identity of the cells, we applied a Bayesian Classifier that assigns each single cell transcriptome to prior annotated transcriptional states (Figs. 1C,D, Fig. S1H). For human cell annotation, we used whole-transcriptome profiles of FACS-sorted sub-populations (Newman et al., 2015). For mouse, we used a comparable data set from the IMMGEN consortium (Heng et al., 2008). Our classification revealed that the clusters corresponded to almost all major known immune cell lineages. In humans, we identified various myeloid cells (mast cells, neutrophils, classical dendritic cells (DCs), plasmacytoid DCs (pDCs), monocytes, macrophages) and lymphoid cells (T cells, NK cells, B cells, plasma cells). In mouse, we identified the same lineages, except for plasma and mast cells. We detected basophils, which were not seen in the human biopsies. In neither organism did the classifier identify a distinct cell cluster as eosinophils, suggesting that these either did not represent a distinct cell state, were extremely rare or failed to survive sample preparation. Expression of genes defining major immune subsets is shown for humans and mice in Figs. 1E, F.

We make two technical notes about collecting these data. First, neutrophils showed very low transcript counts and required setting low filtering thresholds to allow their detection. These cells could be inadvertently excluded using the data filters commonly used in scRNA-Seq studies. Second, plasma cells in patient tumors expressed very low transcript counts for PTPRC, the gene encoding the antigen CD45 (Table S2). This may explain why plasma cells were not detected in the CD45-sorted mouse immune repertoire.

To establish confidence in the ability of scRNA-Seq to report on known cell states, we compared the abundances of cell types identified by scRNA-Seq to conventional gates used in flow cytometric analysis of mouse tissue. This comparison showed strong quantitative agreement between the two methods (Fig. 1G, p<10⁻⁵, Pearson’s R=0.98). We observed that all eight major cell types were represented by 6 of 7 patients (Figs. 1H,I and Fig. S1I), albeit in different proportions that may reflect differences in disease etiology or progression. These findings suggest that the patient cohort represented major aspects of the lung tumor immune microenvironment.

Unbiased comparison of human and mouse immune cells reveals a broad conservation of major gene expression programs

We generated a resource on the gene expression state of immune cells in mouse and human tumors, by cell type. In total, over 4,896 genes showed significant differences in expression in the major cell lineages in human, and a comparable number (4,545) in mouse (Figs. 2A,B, Tables S2,3). Although each immune cell type defined here contains further sub-structure, the gene expression profiles defined at the level of classical cell types reflect many known markers (Figs. 2A,B).

To systematically compare human and mouse immune cells, we quantified the similarity between the average transcriptome of the immune cell types identified in human and mouse (orthology method and gene list in Table S4). This ‘low resolution’ analysis ignored sub-structure within each lineage, and served as a starting point for a detailed analysis within each lineage. We observed that cell type identity dictated the similarity between gene expression profiles, rather than the source organism. This could be seen by clustering on gene expression correlations (Fig. 2C) and by inspecting the expression of cell type-specific genes across cell types in both species (Fig. 2D). With the exception of T and NK cells, each major cell lineage clustered with its ortholog. T and NK cells were more similar to each other within organisms than across organisms, reflecting their proximity in gene expression observed in SPRING visualization (Figs. 1C,D) and expression heatmaps (Figs. 2D, Table S3). The similarities in gene expression also reflected ontogeny, with myeloid cells separating from lymphoid cells. An apparent exception was the pDCs, which are considered a myeloid cell type, but associated here with lymphoid lineages. This observation is consistent with recent observations that pDCs are distinct from conventional DCs (See et al., 2017; Villani et al., 2017).

This correspondence of cell type orthologs reaffirmed that there exist coherent immune programs conserved between human and mouse, and it justified an examination of similarities and differences within each immune lineage.

Mouse and human lung tumors contain common and distinct neutrophil subsets

We analyzed the substructure of each major cluster, focusing on neutrophils, DCs, monocytes and macrophages (see Figs. S2A–R for T, NK, mast, and plasma cell sub-states). We present first data for neutrophils, which showcase examples of conserved and divergent gene expression patterns between mouse and human.

In patient tumor samples, neutrophils formed a continuum of states, which resolved into five subsets by spectral clustering (human N_1–5 or hN_1–5, Figs. 3A–C). The representation of these subsets varied among patients (Figs. 3D,E, Fig. S3A, Table S1). Neutrophils in mouse KP1.9 tumors also formed a continuum of states, which resolved into six subsets (mouse N_1–6 or mN_1–6, Figs. 3F–H). The proportion of these subsets markedly differed in tumor-free and tumor-bearing lungs, with mN_1,2 enriched ~ten-fold in healthy lungs, mN_4,5 increasing 10 to 20-fold in tumor tissue, and mN_3,6 exclusively present in tumor tissue (Fig. 3I, Fig. S3B).

We previously identified two mouse lung neutrophil subsets, characterized by variable expression of the sialic acid binding Ig-like lectin F (Siglecf). SiglecF^low cells are found in tumor-free lungs, whereas SiglecF^high neutrophils accumulate in tumor tissues and exhibit several pro-tumor functions as (Engblom et al., 2017). Here, we identified that three neutrophil subsets (mN_4–6) expressed Siglecf, and indeed these were highly tumor-enriched, whereas the remaining three subsets (mN_1–3) were Siglecf-low (Fig. S3C).

To build confidence in the observed sub-population structure, we sorted SiglecF^high and SiglecF^low cells from lung tissue (Fig. S3D) and examined predicted gene expression markers by qPCR (Fig. S3E). We expected this sorting strategy to enrich for cells in clusters mN_4–6, compared to mN_1–3, and indeed confirmed the predicted overexpression of the Car4, Clec4n, Clec5a, Csf1, Ltc4s, Nt5e, Runx1, Siglecf, Spp1, Vegfa and Xbp1 transcripts in SiglecF^high cells (Figs. S3E–G). We further analyzed Nt5e and Clec5a protein expression by flow cytometry and found strong correspondence between mRNA and protein gene expression patterns (Fig. S3H).

To define how neutrophil subsets relate to each other across species, we compared human and mouse neutrophils systematically on the single cell and whole transcriptome levels (as for Fig. 2C,D). We found conservation of several aspects of neutrophil population structure between species, as can be appreciated through an unsupervised comparison (Figs. 3J–L). In both human and mouse, we observed a common ordering of neutrophil subsets, from h/mN₁ neutrophils that expressed high levels of canonical neutrophil markers (MMP8,9; S100A8,9; ADAM8; synonymous lower case mouse gene symbols omitted here and below), via h/mN_3,4, to h/mN₅ that were tumor-specific in mouse and express cytokines CCL3 and CSF1, as well as CSTB, CTSB, and IRAK2 (Figs. 3B,G, Table S2,3). mN₆ neutrophils were most closely related to their mN₄ counterparts (e.g. they both expressed Hexb and Ptma) although mN₆ uniquely expressed Fcnb and Ngp (Fig. 3G). In both organisms h/mN₂ neutrophils formed a sub-population apart from the continuum of states; these cells were rare but exhibited a strong transcriptional signature, expressing type I interferon response genes, including IFIT1, IRF7 and RSAD2 (top 50 marker genes showed enrichment for the GO term ‘cellular response to Type I interferon’ with FDR=4×10⁻²⁶). We found N₂ cells in both healthy and tumor lung tissue in mice and will later show that they also appear in patients’ peripheral blood.

In summary, we identified three conserved modules of neutrophil gene expression within mouse and human neutrophils (Fig. 3L): a) neutrophils that express canonical neutrophil markers (h/mN₁); these progress continuously to b) neutrophil states that are tumor-specific and promote tumor growth in mice (h/mN₅); and c) a small but distinct neutrophil subset with an expression signature of type I interferon response (h/mN₂).

DCs contain four distinct subsets that are conserved between mouse and human

Using the same Bayesian classifier and reference datasets, the remaining TIM populations were assigned as having gene expression profiles of monocytes, macrophages and DCs (Fig. S1H). We analyzed each subset in separation. Spectral clustering of patient tumor DCs identified four subsets: pDC (human pDC, or hpDC) and hDC_1–3, which showed distinct gene expression programs across tens of genes (Fig. 4A–C, Table S1–3). The DC subsets were all present in 6/7 patients, albeit in variable proportions; one patient lacked several subsets, possibly due to the low cell numbers in this sample (Figs. 4D,E, Fig. S4A).

To assess hDC_1–3, we compared their gene expression profiles to those of bulk-sorted classical DCs, which comprise cDC1 (defined as XCR1⁺ CADM1⁺ by flow cytometry) and cDC2 (CD1A⁺ CD172A⁺) (Eisenbarth, 2019; Guilliams et al., 2016). cDC1 are efficient antigen cross-presenters to CD8⁺ T cells, whereas cDC2 preferentially interact with CD4⁺ T cells. The hDC₁ and hDC₂ subsets expressed gene signatures that mapped well to these DC populations (Fig. 4F). Accordingly, hDC1 uniquely expressed the cDC1 markers XCR1, CLEC9A, TBHD (the gene encoding for cell surface protein CD141) and CADM1, whereas hDC2 expressed cDC2 markers, including CD1A, CD1C, CD1E, CD207 (encoding for langerin) and FCER1A (Fig. 4B, Table S2).

hDC₃ showed an “activated” DC phenotype based on comparisons to previously published microarray datasets of in vitro LPS-stimulated DCs (Fig. 4G) (Zanoni et al., 2009) but lacked expression of key cDC_1/2 and pDC genes. For example, hDC₃ expressed the cDC1-associated genes BATF3 and IRF8 but neither XCR1 nor CLEC9A (Fig. 4H). Individual hDC₃ cells also failed to consistently associate with a single DC signature (Fig. 4I). Many appeared similar to hDC₁, however, and hDC₃ were much closer in expression to DCs than to monocytes or macrophages.

Six human blood DC subsets have been reported (Villani et al., 2017), most of which showed a correspondence to our data with some exceptions (Fig S4B–D). One of the blood clusters (denoted DC5 in that paper) was absent from our tumor data, while the tumor hDC₃ subset was absent in their blood data. These differences likely reflect variation in DC states between healthy blood and tumor tissue.

In mouse (Figs. 4J–O), we again found four distinct DC subsets, all of which were significantly increased in tumor-bearing lungs (Fig. 4M). These four subsets could be related one-to-one to those found in human, first by inspecting marker genes for each cell cluster (Figs. 4K,L, Table S2), and then by unsupervised hierarchical clustering using orthologous variable genes (Figs. 4P–R, Table S2,3). mDC_1–3 mirrored hDC_1–3 while mpDCs mirrored hpDCs.

The gene expression profiles of mDC₁ and mDC₃ respectively resembled those identified for Ccr7–and Ccr7+ cDC1 subsets in the thymus and peripheral lymph nodes of mice (Ardouin et al., 2016) (Fig. S4E, F). These data suggest that Ccr7– mDC₁ can give rise to Ccr7+ mDC₃ cells. mDC₃ also resembled tumor-infiltrating Ccr7+ DCs that secrete interleukin-12 and have anti-tumor activity (Garris et al., 2018) (Fig. S4G, H). CCR7 expression may enable mDC₃ to migrate to draining lymph nodes (Roberts et al., 2016) but some of these cells can remain within tumors (Garris et al. 2018).

In summary, the conserved structure of both the mouse and human DC populations, and its recurrence in other studies, supports the view that DCs can be consistently classified and compared between species (Dutertre et al., 2014). Tumors in both organisms decomposed into plasmacytoid and cDC1–2 subsets, and into resting and activated DCs, and we could identify no others. In both organisms, pDCs appeared to represent a distinct cell state, suggestive of alternate ontogeny to other DCs. Both organisms also presented monocytes showing DC-like signatures, which are described below.

Monocyte subsets are well conserved between mouse and human, while macrophage subsets show species-specific patterns

In human, we defined 14 transcriptional states of cells classified as monocytes and macrophages. Three showed gene signatures of monocytes (hMono_1–3), nine of macrophages (hMø_1–9), one of monocytes and DCs (hMonoDC), and one characterized by cell cycle gene expression (hMøCycl) (Figs. 5A–D, Fig. S5A, Tables S2,3). We assigned these identities to the clusters by comparing their gene expression to that of annotated reference gene expression profiles (Gentles et al., 2015; Newman et al., 2015) (Fig. 5D). Even though spectral clustering identified a relatively large number of transcriptional states, all of them could be detected using any 4/7 patients (Fig. 5C, Figs. S5B,C, Table S1). These results indicate substantial monocyte and macrophage heterogeneity in lung tumors defined by multiple gene expression programs, which are reproducible across individuals.

We examined gene signatures to assess whether any of the states observed correspond to canonical M1 and M2 states, which have been used to define anti-inflammatory and pro-inflammatory macrophages, respectively. In humans, no macrophage cluster exhibited only an M1-like phenotype, whereas hMø_3,6,8, and hMonoDC exhibited an M2 gene signature (Fig. 5E). Enrichment for the M2 signature can be expected for macrophages found in growing tumors, but these data show that M2 is not a distinct state.

In mouse, we defined a total of eight transcriptional states corresponding to monocytes (mMono_1–3), macrophages (mMø_1–4), and both monocytes and DCs (mMonoDC) (Figs. 5F–J, Fig. S5D, Tables S2,3). The monocytic states were well represented in healthy tumor-free tissue, as was a single cluster (mMø₄) identified as a resident alveolar lung macrophage by its expression of Krt79, Krt19 and Car4 (Gautier et al., 2012). The remaining mouse macrophage states were strongly tumor-enriched (Fig. 5H). As with human, we assigned the cluster identities using annotated reference gene expression profiles (Heng et al., 2008) (Fig. 5I). Similarly to our human data, we found that previously annotated M0/½ gene expression signatures (Hou et al., 2018) did not clearly associate with murine macrophage clusters, and that all macrophages exhibited an M2-like phenotype (Fig. 5J).

The complexity of these populations makes it a challenge to identify similarities in population structure between organisms by inspection alone. Nevertheless, a systematic unbiased comparison (Figs. 5K,L, Table S3) identified that hMono_1–3 corresponded respectively to mMono_1–3, with each of the populations sharing distinct gene expression patterns. The MonoDC subset was likewise similar between patient and murine tumors, and was characterized by high gene expression of DC-related genes, including CLEC10A, MHC class II genes and associated genes (e.g. CD74) (Figs. 5B,G,K, L). Both the human and murine MonoDC subsets were distinct from the seemingly more mature DC populations DC_1–3 and pDCs, e.g. human MonoDCs lacked or expressed lower levels of DC marker genes, including XCR1, CD1A, CD1E, FSCN1, and TCF4.

The correspondences noted above reflected canonical distinctions between monocyte subsets, as well as a less appreciated monocyte state. Historically, two types of monocytes have been described in both humans and mice albeit with different markers. In humans, these are identified in blood as CD14⁺ (so called ‘classical’) and CD14^int CD16⁺ (‘non-classical’) (Ginhoux and Jung, 2014). Corresponding ontogenically and functionally related subsets in mice are defined as Ly6C^high CCR2⁺ CX₃CR1^int and Ly6C^low CCR2^– CX₃CR1^high. Ly6C^high monocytes can extravasate into tissues and give rise to macrophages and DCs, whereas Ly6C^low monocytes can remain in the vasculature and patrol vessel walls. In our murine data, mMono₁ expressed the Ly6C^high monocyte markers Ly6c1, Ly6c2 and Ccr2, whereas mMono₂ expressed the Ly6C^low monocyte markers Cx3cr1 and Fcgr4 (Fig. S5E, Table S2), in agreement with previous work (Ingersoll et al., 2010; Mildner et al., 2017). Similarly, hMono₁ expressed classical monocyte-associated genes including CD14 and FCN1, and hMono₂ expressed non-classical markers CDKN1C, LILRB2, and ITGAL (Fig. S5F, Table S2). Our genome-wide comparison of h/mMono subsets (Figs. 5K,L), confirmed the correspondence of the classical and non-classical monocyte types between organisms.

These canonical descriptions left the third cluster, h/mMono₃, unaccounted for. Among all monocytes and macrophages, h/mMono₃ uniquely expressed a set of neutrophil-associated genes (Fig. 5B, D, G, I), including S100A8, S100A9, and CSF3R (Fig. S5G), consistent with previous reports of both bulk (Ingersoll et al., 2010; Schmidl et al., 2014) and single cell analysis of sorted CD14⁺ blood monocytes (Villani et al., 2017). In humans, hMono₃ also expressed markers of CD14⁺ CD16^– monocytes suggesting that this population may be an unappreciated subtype of ‘classical’ monocytes (Fig. S5F). In mice, mMono₃ did not express a clear Ly6C^high or Ly6C^low gene signature; instead, this cluster contained both Ly6C^high and Ly6C^low-like cells. It is therefore plausible that, in mice, the neutrophil-associated genes define a program orthogonal to that of classical and non-classical identity. A fourth monocyte subset was previously identified in peripheral blood (Villani et al., 2017). Many genes expressed by this subset were not present in any of our monocytes and other populations and we suspect that this subset might include physical doublets with NK cells (Fig. S5H).

Macrophages showed both conserved and divergent phenotypes between mouse and human. In both organisms, key macrophage-associated genes, such as MRC1, APOE, complement genes (C1QA,B,C) and cathepsins (CTSB,D) were enriched in macrophages compared to other TIMs and broadly expressed across clusters in continuous gradients (Figs. S5I, J). To relate macrophage heterogeneity across species, we asked whether gene signatures of mouse macrophage subsets varied across human macrophage clusters (Fig. 5M). We found that mMø₁ and mMø₂ gene signatures mapped to a coherent gradient across human macrophage subsets, especially subsets hMø_8,9 and hMø_3,4,5, respectively. Similarly, the mMø₄ (alveolar macrophage) gene signature mapped mainly to hMø₇. In contrast, the mMø₃ gene signature did not map well to human macrophage subsets, suggestive of a subset unique to the murine tumors. Conversely, hMø_1–2 did not appear well represented by murine macrophages. These complex relationships reinforce the difficulty in relating macrophage subsets between species. However, our analyses provide evidence that transcriptional programs broadly distinguishing tumor-associated macrophages from monocytes are well-conserved between mouse and human and that features of distinct murine tumor-infiltrating macrophage subsets may correspond to several cell subsets in patient tumors. We summarize the key similarities and differences in Fig. 5N.

Among the human macrophage subsets, we identified distinct chemokine expression (Fig. S5K, Tables S2,3) including hMø₁ expressing the neutrophil chemoattractant CXCL5, hMø₈ expressing the CXCR4 ligand CXCL12, and hMø₉ expressing the T cell recruiting chemokines CXCL9,10,11. The respective receptors for these macrophage-expressed chemokines displayed heterogeneous expression patterns across other immune cell types. For example, CXCR2, encoding the receptor for CXCL5, is expressed exclusively by neutrophils, whereas CXCR3, encoding the receptor for CXCL9,10,11 is mainly expressed by T cells, NK cells and pDCs. In mice, we did not observe the same pattern; instead, murine mMø (mMø₁ and mMø₃ in particular) expressed mainly monocyte and neutrophil chemoattractants, while expressing negligible levels of the T cell chemokines CXCL9,10,11 (Fig. S5L, Tables S2,3). These findings suggest that macrophage diversity might be explained in part by their migration cues directed at cell populations with distinct chemokine receptor expression.

Defining distinct marker genes for TIM subsets and their association with patient survival

By analyzing seven patients, a question was whether the observed states were representative of a wider population. In examining non-immune proliferating epithelial cell states, which likely includes tumor cells within our samples (Fig. 6A), we found that each patient exhibited distinct cell states, indicating that our sample size is not adequate to comment on the diversity of patient epithelial phenotypes. By contrast, TIM states overlapped between patients (Fig. 6B), indicating that the myeloid microenvironment is much more stereotyped than the tumor tissue that it infiltrates. This observation gives confidence that we captured reproducible subsets within the myeloid lineages.

The overlap of TIM states between patients motivated us to assess their association with patient survival, by examining the expression of genes specific to each sub-population in survival studies. For this, we determined which genes were enriched within each state relative to all other cells in human tumors, taking all tumor-associated cells (non-immune and immune) into account (Figs. 6C, Tables S2,3). Some specific genes showed clear association with patient survival: for example, both the hN₂ marker ISG15 and the hN₅ marker PI3 associated negatively with patient survival, whereas the hDC₂ marker CD207 associated positively with survival (Fig. 6D). We systematically scored marker gene prognostic association over multiple studies (n=1,127 patients) by PRECOG (Gentles et al., 2015) (Fig. 6E) and summarized the resulting associations in Fig. 6F. This analysis confirmed that genes enriched in CD8⁺ T cell and B cell subsets associated consistently with better patient survival, as expected (Fig. 6E, F). For some sub-populations, no single gene was expressed with sufficient specificity to directly relate its expression to population abundance (Fig. 6E). Among neutrophils, DCs and monocyte/macrophage clusters, many subsets showed specific gene expression that provided a remarkably consistent view of the prognostic association of each sub-population. The neutrophil subsets hN₂ and hN₅ both showed an abundance of marker genes associated with poor patient survival, as did hMø₁ and hMø_9, while among myeloid cells, hDC₂ showed the most consistent positive association with survival, with the other DC subsets and a mast cell subset also showing positive associations. These results support findings of previous association studies and extend them to resolution of subsets within each major cell lineage.

Access to marker genes (Figs. 6C,E) allows developing strategies to detect immune cell subsets in histological sections, opening an opportunity for future large-scale cytometric analysis, e.g. using existing patient sample banks. We tested feasibility of such analyses by examining markers specific to hN₅ neutrophils in situ. This subset represented one end of a continuum of neutrophil states (Fig. 3A), was associated most strongly with negative patient survival and was conserved between mouse and human. Within the patient cohort, we identified two patients (Pt 3 and Pt 7) that showed extreme differences in the abundance of hN₅ cells, with Pt 3 hosting few of them and Pt 7 hosting many (Fig. 6G, Table S1). We stained tissue sections for transcripts of PI3, which marks hN₅. We additionally stained for a pan-neutrophil marker, the primary granule protein myeloperoxidase, or MPO (Eruslanov et al., 2014). As predicted, Pt 7 hosted an abundance of PI3⁺ cells while Pt 3 hosted almost none (Fig. 6H). As predicted, both patients had detectable numbers of MPO⁺ cells (Fig. 6I) and PI3+ cells were MPO+ (Fig. 6J). These results validate the trends seen by scRNA-Seq, and demonstrate the utility of the latter for designing histology-based studies.

Genes targeted by mouse models, immunotherapies and immunophenotyping across murine and patient tumor-infiltrating immune cells

Defining the TIM landscape offered an opportunity to ask how genetic reagents, drugs and antibodies might target different myeloid cell populations in mice and humans. We cataloged 34 genes targeted in genetically engineered mouse models, e.g. used for fluorescent cell tracing, cell ablation, or other transgene expression (Figs. S6A–F, Table S5). We also examined 26 genes that encode proteins targeted by immunotherapeutic agents (Fig. S6G, Table S5) and 36 genes whose products encode surface antigens used in immunophenotyping (Fig. S6H, Table S5). For each of these gene sets, we assessed how well gene expression agreed between mouse and human across immune cells, and how specific gene expression was to defined cell types. Though many of the genes targeted in genetically engineered mouse models and by drugs showed a good correspondence in expression between organisms, there were notable exceptions that highlight inter-species differences. For surface antigens, as well as identifying similarities and differences between organisms, it was notable that many genes canonically used to identify specific cell types showed expression over multiple cell types. For all of these cases, we refer to Fig. S6 for notable examples. Overall, these results highlight the need to benchmark specific mouse immunotherapy models against human as performed here. The findings also indicate the potential shortcomings of quantifying or purifying cell types by immunocytometry, which may group distinct cell states together even with some of the best-established markers.

Tumor infiltrates only partially overlap with states of the peripheral blood

We asked whether TIM cell states can be appreciated by sampling peripheral blood (PB), which is a more easily accessible tissue and could be useful for diagnostic purposes. We began to address myeloid cell states’ tumor specificity by comparing them to PB myeloid cells of the same patients (we profiled 6 patients). Our goal was to define which tumor myeloid states had a circulating counterpart.

Several scRNA-Seq studies have previously examined blood immune populations focusing on PB mononuclear cells (PBMCs), which are fractionated by their density [e.g. (Zheng et al., 2017)]. PBMC preparations are deliberately depleted for polymorphonuclear cells, which include all granulocytes, and thus provide an incomplete view of myeloid cells. Therefore, we prepared blood for analysis by red blood cell depletion to obtain all other immune cells types. After filtering we collected 14,411 cell transcriptomes by scRNA-Seq.

SPRING embedding of cells from PB and tumors across all patients revealed that a few tumor and PB immune cells overlapped in state, although most did not (Fig. 7A). This separation did not reflect library batch effects, because all patient samples overlapped within the blood and tumor clusters. Therefore, PB cells were transcriptionally separable from those in tumors, but could still reflect specific tumor cell states. Applying the cell type classifier (as in Fig. 1C), we identified the expected array of blood leukocytes, including monocytes, neutrophils, pDC, T, NK and B cells (Fig. 7B). Average gene expression profiles for these blood populations are provided as a resource in Table S2 and Fig. S7A. The proportions of each cell type differed between the two tissues and some populations were only detected in the tumor (Fig. S7B).

We asked whether myeloid cell sub-clusters were similar across tumor and blood tissue. Clustering of PB monocyte populations revealed three monocyte states (hbMono_1–3) and one monocyte-like state expressing concomitant dendritic-cell associated genes (hbMonoDC) (Fig. 7C, Figs. S7C–F). PB neutrophils formed a continuum, which we divided into six subsets (hbN_1-6) by spectral clustering (Fig. 7D, Figs. S7G–J). For monocytes, hierarchical clustering revealed a one-to-one match between tumor and blood populations (Figs. 7E, F, Table S3), although the tumor and blood Mono₁ subsets shared few marker genes. Part of this correspondence was anticipated, considering that monocytes are defined as circulating cells. However, blood and tumor monocytes differed (409 genes >2-fold differentially expressed at 5% FDR, Fig. 7G, Table S2). PB monocytes, for example, expressed significantly more transcripts for the chemokine receptor CCR2 while tumor monocytes expressed CXCL3, a cognate ligand for the CXCR2 receptor expressed by neutrophils (Fig. 7H). For reference, we also compared the patient blood monocyte subsets to those from healthy donors described in (Villani et al., 2017) [here denoted Mono (Vil.)] (Fig. S4B, C). We found a close correspondence between Mono1,2,3 (Vil.) and hbMono_1,2,3. Mono4 (Vil.) expressing NK-associated genes, as previously defined by Villani et al., was not detected in our dataset. Additionally, one subset thought to be DCs in that study (DC4) showed a correspondence to hbMono₂, consistent with a previous report (Calzetti et al., 2018).

Among neutrophils, the blood subset hbN₂ corresponded strongly to the tumor subset hN₂, with a similar Type I interferon signature (IFIT1,2,3, RSAD2), and the blood subsets hbN_1,3,6 were related to tumor clusters hN_1,3 and shared the expression of elevated levels of the canonical neutrophil-associated genes MMP9 and S100A8,9,12 (Figs. 7I,J, Table S3). Yet there was no strong PB correspondence to the CCL3⁺ tumor subset hN₅. Overall, blood and tumor neutrophils strongly differed in their gene expression (719 genes >2-fold differentially expressed at 5% FDR Fig. 7K, Table S2). For example, tumor neutrophils expressed far higher levels of the cytokine CXCL8, and much lower levels of the IL17 receptor, IL17RA, than their blood counterparts (Figs. 7L).

PB cells also contained hbN₆ neutrophils, which were undetected or very rare in the tumor (Figs. 7D,I). The hbN₆ subset was found in all patients and represented 6.0 ± 0.4% (mean±SEM) of blood neutrophils (Fig. S7H). It expressed unique markers, including TNFRSF13B, which encodes the receptor for the cytokines APRIL and BAFF, as well as PDE10A and ADRA1A. TNFRSF13B is associated with B cell lineage differentiation, class switching and antibody production (Rickert et al., 2011), and was expressed by tumor, but not PB, B cells. The function of circulating TNFRSF13B⁺ neutrophils is unknown, though mutations in TNFRSF13B are implicated in a common variable immunodeficiency that can in rare cases present neutropenia with poorly understood etiology (Bogaert et al., 2016; Guffroy et al., 2017). The gene encoding BAFF, which binds TNFRSF13B, was expressed by other blood neutrophils (bN_3–4), suggesting possible signaling between PB neutrophils. In summary, PB neutrophils show some reflection of population structure found in the tumor, but both the tumor and PB host non-overlapping and unappreciated neutrophil subsets. We conclude that to understand TIMs it is necessary to directly assess tumor tissues, as the PB does not recapitulate the tumor myeloid cell landscape.

CONTACT FOR REAGENT AND RESOURCE SHARING

Further information and requests for resources and reagents should be directed to and will be fulfilled by the Lead Contact, Allon M. Klein (allon_klein@hms.harvard.edu).

Raw RNA sequence data from human samples are not available publicly due to patient privacy concerns.

EXPERIMENTAL MODEL AND SUBJECT DETAILS

METHOD DETAILS

QUANTIFICATION AND STATISTICAL ANALYSIS