Single-Cell Transcriptome Analyses Reveal Endothelial Cell Heterogeneity in Tumors and Changes following Antiangiogenic Treatment

Solid tumors are composed of a complex assortment of tumor cells and various stromal cells, including endothelial cells (EC), fibroblasts, and smooth muscle cells, as well as immune cells, comprising the cellular tumor microenvironment. Stromal cells are intimately involved in tumor growth and progression, and as such have increasingly become the targets of anticancer therapies. For example, the growth of new blood vessels (angiogenesis) is important for tumor growth (1), and a number of antiangiogenic agents are in clinical use or are being tested (2). Additionally, fibroblasts promote a protumor microenvironment via various mechanisms (3), one of which involves deposition of extracellular matrix components (ECM; ref. 4) that contribute to tumor progression (5, 6).

In addition to the various cell types within the tumor microenvironment, there is heterogeneity within each cellular compartment. Some dynamic interactions among the heterogeneous cells are critical to tumor growth. For example, in tumor angiogenesis, dynamic interactions between specialized ECs allow some ECs to sprout and migrate from a blood vessel (so-called tip cells), whereas other cells remain relatively more static and form the shaft behind the sprout (so-called stalk cells). Key angiogenic signaling pathways, such as VEGF and Dll4-Notch, regulate these interactions between ECs and help to shape the cell specialization. While specific blockers of these key angiogenic pathways show antitumor effects and reduce tumor perfusion, their particular effects on the processes of tumor angiogenesis differ. For example, VEGF inhibition reduces vessel density (7, 8), whereas Dll4 blockade results in increased, but nonfunctional tumor vasculature (9, 10).

The heterogeneity of tumor stromal cells has long been recognized, and various approaches have been used to document the different cell subpopulations. Histologic and immunostaining methods have been widely applied, and using tumor angiogenesis as an example, these approaches have shown specific expression of genes in specialized ECs such as apelin (Apln) and Esm1 in tip cells (11). Other approaches to study heterogeneity of tumor stromal cells include laser dissection/capture and flow cytometry (12, 13). However, these approaches require preknowledge of the specific marker genes to be visualized. The effects of VEGF or Dll4-Notch blockade on ECs have been studied by various groups using morphologic analyses, molecular characterization as well as gene profiling (2, 8, 11, 14–16). However, these approaches focused on whole tumor or tissue lysates, or in some cases on isolated EC populations using flow cytometry or laser microdissection, all of which provide pooled ECs instead of providing data at a single-cell level. Thus, there is need to better characterize cellular heterogeneity of tumor stromal cells. Recent technology developments now allow RNA profiling of individual cells derived from tumors or tissue, which provides more information on single-cell properties and heterogeneity. The single-cell transcriptome technique allows dissection of different cell populations in the tumor microenvironment and capture of changes even in rare cell subpopulations that would be missed by bulk RNA-Seq.

Here, we used single-cell transcriptome profiling to not only characterize tumor EC heterogeneity and compare them with normal ECs, but also to evaluate changes in EC subpopulations following blockade of key angiogenic signaling pathways, namely VEGF and Dll4-Notch. In addition, we phenotypically classified other stromal cells in the tumor microenvironment, particularly smooth muscle cells, pericytes, and fibroblasts. To achieve this, we isolated single cells from xenograft tumor-bearing mice that had been treated with the anti-VEGF agent aflibercept (VEGF Trap), an anti-Dll4 Ab (aDll4), a combination of aflibercept and aDll4, or control hFc. In parallel, single cells from similarly treated normal hearts, which are an abundant source of nonangiogenic endothelial cells, were also collected to serve as a normal comparator in this study. Although heart endothelial cells may express specialized transcription factors, angiocrine factors and surface markers (17, 18), these vessels lack fenestrations and other specialized morphologic features and were, therefore, chosen as representative of quiescent vasculature. Single-cell transcriptome sequencing using the 10× Genomics platform allowed analyses of individual cells. Our findings better define the heterogeneity of tumor ECs and other stromal cells. In addition, the results show the roles of VEGF and Dll4-Notch signaling in specifying tumor EC phenotype, and highlight the response of stromal cells to antiangiogenic therapies.

Animal studies were approved by Regeneron's Institutional Animal Care and Use Committee (IACUC). Tumor cell lines were obtained from the ATCC. All cell lines were IMAPCT tested (including mycoplasma) by IDEXX BioResearch and authenticated between 2012 and 2015 using the STR Profile Testing by ATCC. Cells used are less than 5 passages in our experiments. A total of 2 × 10⁶ COLO205 human colon carcinoma or 2.5 × 10⁶ HT-29 human colon carcinoma cells were grown subcutaneously (s.c.) in male CB.17/SCID mice (Taconic). When tumors reached 100–150 mm³, mice remained untreated or were treated by s.c. injection with hFc (control protein, 25 mg/kg), aflibercept (VEGF Trap, ziv-aflibercept, 25 mg/kg), a DII4 blocker (REGN 1035, 10 mg/kg) or a combination of aflibercept plus anti-DII4 antibody. Twenty-four-hour or 72-hour posttreatment heart and tumor tissue was harvested.

Tissues were minced, enzymatically digested for 10 (COLO205 tumor) or 15 minutes (HT-29 tumor, heart) using 2.5 mg/mL collagenase Type II (Worthington), 2.5 mg/mL collagenase Type IV (Gibco), and 0.5 mg/mL DNase (Sigma) in PBS/1% BSA at 37°C, treated with DMEM containing 10% FBS and filtered (70-μm nylon filter; Falcon). Single-cell suspensions were subjected to 1X PharMLyse (BD Biosciences) for 4 (tumor) or 10 (heart) minutes, washed, and subjected to dead cell removal kit (130-090-101, Miltenyi; manufacturer's instruction). Tumor single-cell suspensions underwent a depletion step using the following FITC–conjugated Abs: COLO205 cells were incubated with anti-CD24 (11-0247-42, eBioscience) and anti–E-cadherin (324104, Biolegend) and HT-29 cells were incubated with anti-CD44 (338804, Biolegend) and anti-EpCam (324204, Biolegend) antibodies using 10 μL anti-FITC-beads for 1 × 10⁶ total cells. The remaining cells were used for single-cell sequencing.

Cellular suspensions (∼6000 cells) were loaded on a Chromium Single Cell Instrument (10X Genomics) to generate single-cell GEMs. Single-cell RNA-seq libraries were prepared using version 1 Chromium Single-cell 3′ Library, Gel Bead & Mutiplex Kit (10X Genomics). Sequencing was performed on Illumina NextSeq500 containing transcript length of 59 bp. The Cell Ranger Suite version 1.1.0 was used to perform sample demultiplexing, barcode processing and single-cell gene UMI (unique molecular index) counting (http://software.10xgenomics.com/single-cell/overview/welcome). Data are accessible through GEO (GSE110501).

Single cells were filtered for downstream analysis by the following criteria: UMI (unique molecular identifier) count is within the range between 3,000 and 40,000, and mitochondria percentage is less than 10 percent of the total UMI count. Gene expression (in UMI) is scale normalized then transformed in log₂ (UMI+1).

All treatment experiments were duplicated. We used the dataset from one experiment for discovery and the other for result validation. For principal component analysis (PCA), the analysis was run on normalized and transformed UMI counts on variably expressed genes, which were defined as genes expressed with at least 10 UMI counts in at least three percent of all cells. These PC outputs were loaded as input to generate t-SNE plots. 3D t-SNE plot was drawn using Rtsne package with default parameters. 2D t-SNE plot was generated with Seurat package. Two hundred eighty-four to 667 highly variable genes were selected in different analyses with the first 10 or 15 PCs applied to find clusters.

To identify genes that specifically expressed in a cell population or cluster, one-way ANOVA test and F-test for multi-group comparison implemented in ArrayStudio (www.omicsoft.com/array-studio/) were performed. Genes are ranked by the FDR-BH corrected P values. Consensus clustering (https://genepattern.broadinstitute.org/) was used to assign ECs into tip-like, stalk-like, or transition subpopulations based on the 60-gene expression signature.

Z-score was calculated on log₂(UMI + 1) for each of the 88 genes in the cell-cycle gene panel. The ceiling of a positive z-scores was set at 2.5; the ceiling for a negative z-score was set at −2.5. Then, z-score values were combined for the 88 cell-cycle genes for each single cell. The cutoff of combined z-score value for calling a cell under active cell cycle was determined by k-means clustering (K = 2) with 200 resampling iterations. P values were calculated by two-tailed Fisher exact test for enrichment of cells in active cell cycle between drug-treated and control-treated groups. For differentially expressed genes, one-way ANOVA (assuming equal variance) was used on normalized and log₂ transformed UMI counts. Gene Ontology (GO) and pathway enrichment (using KEGG database) analysis were performed in DAVID v6.8 (https://david.ncifcrf.gov).

For FFPE sections, tissues were immersion-fixed in 10% neutral-buffered formalin, dehydrated through graded ethanol and xylenes, embedded in paraffin and cut into 4-μm (tumor) or 6-μm (heart) sections. The sections were air dried, deparaffinized, incubated in H₂O₂, subjected to target retrieval procedure. For ISH, sections were subjected to protease plus treatment before the required six hybridization/amplification steps with the respective probes according to the manufacturer's recommendations. Positive RNA detection was revealed with 3,3′-diaminobenzidine (DAB) or alkaline phosphatase (AP). Slides were then counterstained with hematoxylin, dehydrated, and coverslipped. For IHC, sections were stained using an anti-CD31 Ab (Abcam, ab28364) or anti-Lyve1 Ab (in-house Regeneron Ab) and protein positivity was revealed with DAB.

Before analyzing tumor endothelial and stromal cells, we first characterized the heterogeneity of cells in a normal organ from adult mice, namely heart where the response of ECs to antiangiogenic treatment is well established (19). Single cells were collected and RNA-seq libraries were made using the 10x Genomics platform. Cell-clustering analysis with Seurat (20) using 307 highly variable genes revealed four distinct cell populations (Fig. 1A), which represent ECs (clusters 1 and 2 share similar marker genes), fibroblasts (cluster 0), leukocytes (cluster 4), and muscle cells (cluster 3) based on expression of known marker genes, such as Pecam1, Col1a1, Lyz2, and Myl9, respectively. The population classified as muscle cells was mainly composed of smooth muscle cells, as distinguished by selective expression of the smooth muscle actin gene (Acta2). Using the current protocol, very few cardiac muscle cells were captured, likely due to difficulty of single-cell disassociation and/or unusual cell size.

We further interrogated the EC population with a total of 838 isolated heart ECs from one sequencing batch. Initial cell clustering analyses generated five possible subpopulations (Supplementary Fig. S1A). However, a heatmap with the combined top 25 cluster-discriminative genes showed that clusters 0 and 1 were very similar based on their sharing of gene signature and were thus merged into one cluster/subpopulation (Supplementary Fig. S1B). Multi-group ANOVA test generated a list of 90 most variable genes that could distinguish the four endothelial subpopulations (Fig. 1B and Supplementary Table S1). The presence and lineage connections of the four subpopulations within heart ECs were further supported by a PCA using the whole transcriptome gene signature followed by 3D t-distributed stochastic neighbor embedding (t-SNE) projection (21; Fig. 1B). For example, lymphatic ECs (EC4) shared an overlapping gene signature with venous ECs (EC3), and were closer to venous EC in the whole transcriptome 3D t-SNE map. This similarity is consistent with lymphatic ECs originating from venous ECs during development (22). These four subpopulations and subpopulation-discriminative gene signatures were confirmed in other independent experiments (Supplementary Fig. S1C).

As the 838 ECs for the above analysis were pooled from 24-hour treatments (see below), we also looked at the distribution of treatment conditions among the subpopulations. We did not observe any skewed distribution of cells with a particular treatment in any subpopulation. Thus, the four subpopulations identified were not treatment related.

Subpopulation-distinctive genes that have preferential expression in ECs compared with other types of cells in the heart microenvironment are listed in Supplementary Table S2. Guided by annotated EC lineage-specific markers in adult mice (23–26), three out of the four EC subpopulations were identified as lymphatic (EC4), venule (EC3), and arteriole (EC1) derived ECs (Fig. 1C; Supplementary Fig. S2A–S2O). For example, endomucin (Emcn) was highly expressed in venule ECs, Flt4, and Lyve1 were more restricted to lymphatic endothelium, and expression of Sox17 was enriched in arterioles. The fourth subpopulation (EC2), which was the largest and exclusively expressed Cd300lg, appears to represent capillary endothelium (Fig. 1C; ref. 27). Dll4 was detected in heart capillary cells, but was not significantly enriched compared with other endothelial subpopulations (Supplementary Fig. S2E).

Certain ephrin ligands and their Eph receptors have been described as markers for arteries and veins during mouse embryonic development, respectively (28, 29). We examined all ephrin/Eph ligand and receptor genes in the four EC subpopulations and found that their expressions were either below detection level or not significantly associated with any endothelial subpopulations in adult heart with some examples shown in Supplementary Fig. S2. Instead, we identified other novel markers distinctive for subpopulations of heart vasculature. For example, Ramp3 was specifically expressed in capillaries, Mmrn1 was restricted to lymphatics, and Vwf was significantly enriched in venous ECs (Fig. 1C). We further validated two novel EC subpopulation-restricted markers identified in this study by histologic methods (RNAScope), which clearly showed the venule-specific expression for Bace2 in a subset of CD31 positive EC structures and lymphatic-specific expression for Reln in Lyve1-positive lymphatic EC structures (Fig. 1D). Taken together, these single-cell transcriptome analyses provide a better resolution of EC heterogeneity in mouse heart.

Given the insights gained from exploring endothelial cellular heterogeneity in the heart, we designed a systematic experimental approach to study stromal cell heterogeneity as well as responses to antiangiogenic agents in tumors and heart (Fig. 2A). To obtain tumor stromal cells, immunodeficient mice were engrafted with either COLO205 or HT-29 tumor cells. Once tumors were established, mice were either sacrificed or treated with a single injection of control protein (human hFc, 25 mg/kg), aflibercept (afl, 25 mg/kg), aDll4 (10 mg/kg), or a combination of aflibercept plus aDll4 Ab for 24 or 72 hours. Heart tissue from the same mice was collected to study treatment effects on normal ECs. Each experiment was performed in duplicate for cross validation.

In xenograft tumors, host-derived mouse stromal cells could be differentiated from human tumor cells based on preponderance of RNA-seq reads that mapped to the mouse rather than the human reference genome. In initial experiments, the majority of isolated single cells from COLO205 tumors were tumor cells, whereas only a small fraction of cells were stromal cells (Supplementary Fig. S3A, left). To enrich stromal cells, we used antibodies to deplete most of the human tumor cells (Supplementary Fig. S3A, right). The depletion procedure also removed some of the leukocytes, allowing our analyses to focus on ECs and fibroblasts. For all the mouse single cells included in our analyses, a median of 77,920 reads were mapped to exonic regions and an average of 1,142 genes were detected per cell.

Unsupervised analysis identified three major cell populations in the enriched tumor stromal cell populations, namely fibroblasts-like cells, leukocytes, and ECs (Supplementary Fig. S3B). Using one-way ANOVA test, tumor stromal cell type-specific genes were identified and were compared with cell markers as seen in the heart (Fig. 2B). With similar total read counts per cell, average numbers of expressed genes captured in these three cell types were significantly different, with leukocytes having the lowest number of genes detected and fibroblast-like cells having the highest number of genes (Fig. 2C). Although the leukocyte population was smaller than the fibroblast-like or EC populations following tumor cell depletion, we were still able to include the remaining leukocytes in some of our analyses.

Focusing initially on untreated COLO205 xenograft tumors, two EC subpopulations were revealed by unsupervised hierarchical and Seurat clustering tool (using a total of 183 ECs; Supplementary Fig. S4A). Genes distinctive to each cluster were identified by one-way ANOVA test. One cluster, which preferentially expressed known tip cell genes such as Kcne3, Nid2, and Dll4, may represent endothelial tip-like cells, whereas the other cluster, which dominantly expressed known stalk EC genes Vwf and Selp, may represent nonsprouting cells or stalk-like cells (Supplementary Fig. S4B). We next performed a similar analysis on another set of 256 ECs from control hFc-treated COLO205 xenograft tumors and obtained similar results. By joining the two ANOVA test results, we derived a 60-gene signature for best classifying tip-like and stalk-like EC subpopulations (Table 1). Besides known marker genes such as Kcne3 for tip-like cells and Selp for stalk-like cells, we were able to identify additional potential tip-like cell markers such as Ramp3, Ednrb, and Cldn5 as well as stalk-like cell markers such as Ackr1 and Tmem252. Ramp3 and Ackr1 showed EC-restricted expression and were further validated by histologic approaches (RNAScope; Fig. 3A).

Using this 60-gene signature, ECs from COLO205 tumors treated for 24 hours or 72 hours (all treatments combined) were consistently clustered into three distinct subpopulations, namely tip-like, stalk-like, and an intermediate cell population, which we termed transition cells (Fig. 3B, Supplementary Fig. S4C). Transition cells expressed both tip-like and stalk-like cell markers, apparently representing a population at a stage in-between those two phenotypes. Notably, ECs collected at 72-hour posttreatment (Supplementary Fig. S4C) showed more variability within tip-like as well as stalk-like cell populations compared with the 24-hour-treated ECs. For example, some stalk-like cells did not express Vwf, Selp, and Ackr1, whereas expression of Adm and Ankrd37 became mutually exclusive from Ptn and Ramp3 in tip-like cells. This increase in diversity did not seem to be a result of treatment, as treatment conditions were not significantly associated with cell distribution.

In order to determine if our findings were more broadly applicable to other tumor models, we conducted similar experiments and analyses on single cells isolated from a different xenograft tumor, HT-29. Similar to COLO205 tumors, two EC subpopulations were identified in HT-29 tumors based on the 40 EC cells analyzed (Supplementary Fig. S4D). The subpopulation gene signatures substantially overlapped between COLO205 and HT-29 xenografts. At least 49 genes out of the 60-gene signature derived from COLO205 tumors were associated with tip-like or stalk-like cells in HT-29 tumors, suggesting that the tip-like versus stalk-like cell gene signature is largely valid in different tumor models. However, there was little concordance between the genes defining the four EC subpopulations in the heart and those defining the tumor EC subpopulations. Also, no lymphatic ECs were detected in a total of more than 2,900 ECs from either COLO205 or HT-29 tumors. Supplementary Table S3 lists 19 genes out of the 60 tip-stalk gene signature from COLO205 tumors that also showed subpopulation preference in heart ECs without consistent associations (Fig. 3C). Tumor endothelial tip-like cell markers Kcne3 and Kcnj8 were rarely detected in heart ECs.

After developing a gene signature to distinguish tip-like, transition and stalk-like cells in tumors, we assessed how different antiangiogenic treatments affect the different subpopulations. COLO205 tumor-bearing mice were treated with control human Fc protein (hFc), the VEGF blocker aflibercept (afl), aDll4, or a combination of aflibercept and aDll4 (combo) for 24 hours or 72 hours, and single cells were analyzed. The relative number of tip-like, transition and stalk-like ECs isolated following the different treatments were quite different. For example, the relative number of tip-like cells increased upon aDll4 treatment at 24-hour posttreatment. In contrast, aflibercept treatment strongly reduced the proportion of endothelial tip-like cells at both 24-hour and 72-hour posttreatment, whereas blockade of both pathways resulted in a somewhat less dramatic decrease in tip-like cells (Figs. 3B and 4A). Although the number and types of cells captured could be affected by many factors, this pattern was reproduced in duplicate experiments, suggesting that the changes were due to the respective treatments.

The changes in EC subpopulations could be a result of changes in cell proliferation caused by aDll4 or aflibercept treatments. To address this question, we examined EC proliferation using a panel of cell-cycle–related genes (30). Aflibercept treatment, alone or in combination, caused decreased cell-cycle activity in all three EC subpopulations compared with control hFc treatment at both 24- and 72-hour time points. In contrast, aDll4 treatment promoted cell proliferation across all EC subpopulations beyond what was observed in the control Fc-treated group at 24-hour posttreatment (Fig. 4B, top). The growth stimulation effect of aDll4 treatment diminished at 72-hour posttreatment (Fig. 4B, bottom). In order to quantify cell-cycle activity, we calculated a z-score for each cell-cycle gene and assigned a combined z-score to each cell to represent the overall cell-cycle activity in each cell. ECs were assigned to active or inactive cell-cycle groups based on the combined z-score. Supplementary Table S4 summarizes the cell-cycle activity within each treatment group per endothelial subpopulation. EC cycle activity was significantly decreased by aflibercept (alone or in combination) in all subpopulations at both 24 hours and 72 hours, whereas it was increased by aDll4 treatment at 24-hour posttreatment.

We also looked at differentially expressed genes between drug and control-treated cells within tip-like, transition, and stalk-like EC subpopulations. Due to reduced cell counts upon certain treatments, these analyses were limited to groups with at least 50 ECs. In general, the drug treatments had greater impact on tip-like cells than stalk-like cells. In tip-like cells, aflibercept and aDll4 single agent led to expression changes in similar sets of genes highlighting upregulation of hypoxia response and glycolysis and downregulation of insulin growth factor signaling (Supplementary Table S5; Fig. 4C; Supplementary Fig. S5). Hilpda, which encodes a hypoxia-inducible lipid droplet associated protein, was significantly upregulated in tip-like cells at 24-hour and 72-hour posttreatment, supporting the previous finding that these treatments led to decreased oxygen delivery in the tumor microenvironment (9). On the other hand, cyclin D1 (Ccnd1) was significantly downregulated in stalk and transition cells upon aflibercept or combination treatment at 24 hours, confirming that these treatments resulted in inhibition of EC proliferation. At 24-hour posttreatment, a number of tip genes defined in this study, including Adm, Ankrd37, Cotl1, Rgcc, and Kcne3, were significantly upregulated in tip-like cells upon aDll4 treatment. There are also drug-specific gene signatures such as upregulation of Tgfb1 only observed in aflibercept-treated tip-like cells. Igfbp7, Ybx1, and Plvap were consistently downregulated genes among all treatment conditions.

There have been reports that loss of Notch signaling is associated with an increase in expression of Vegf receptors (Vegfr2/Kdr and Vegfr3/Flt4), for example, in tip cells, but we did not observe such an increase at either 24-hour or 72-hour posttreatment. In contrast, placental growth factor (Pgf) was persistently upregulated in tip-like cells following drug treatments. Unlike Notch4 and Flt4, which were mainly restricted to tip and transition cells, Notch1 and Kdr were expressed in all ECs.

We next examined changes in heart ECs upon drug treatment. Although the EC subpopulations did not show an obvious shift as seen in tumor ECs, cell-cycle activity increased upon aDll4 treatment at 72-hour posttreatment (Fig. 5A, bottom). The ECs undergoing cell-cycle activity were more restricted to the capillary subpopulation. These changes were not observed at the 24-hour time point (Fig. 5A, top), suggesting that the effects of Dll4 blockade on heart ECs took longer than 24 hours to manifest as changes in gene expression under our experimental condition. Because aDll4-induced capillary proliferation was a similar phenomenon seen in tumor tip cells, we examined if aDll4 treatment could invoke tip cell properties in heart ECs by looking at differentially expressed genes in ECs between aDll4 and control hFc-treated subpopulations (Supplementary Table S6). Interestingly, four out of the top 16 differentially expressed genes ordered by FDR-adjusted P values (P < 1e-5) were tip cell marker genes defined above including Adm, Lxn, Fscn1, and Col4a1. These four genes were significantly upregulated in heart EC2 cells (capillaries) upon aDll4 treatment at 72 hours (Fig. 5B), and at the 24-hour time point showed a trend toward upregulation. In contrast, aflibercept treatment led to a slight decrease in expression of these four genes (Fig. 5B). Taken together, these data suggest that the reported changes in heart tissue after prolonged Dll4 blockade (19, 31) already start on a molecular level 24-hour posttreatment.

In addition to ECs, we also obtained other stromal cells including tumor-associated fibroblasts (TAF) from COLO205 and HT-29 tumors. To interrogate stromal fibroblast subpopulations, TAFs from COLO205 were grouped into six clusters by Seurat, which formed into three distinct subpopulations, specifically fibroblast-like cells, smooth muscle-like cells, and pericytes (Fig. 6A; Supplementary Fig. S6A–S6C). Figure 6B shows expression of known cell type markers in the three subpopulations. Cells in cluster 1, 2, 4 express Dcn, suggesting they are fibroblast-like cells; cells in cluster 0 and 3 highly express smooth muscle actin (Acta2) and Tagln, indicating that they are smooth muscle-like cells. And the cells in cluster 5 express known pericyte markers Rgs5 and Des, suggesting they are pericytes. Fap, a gene expressed by activated fibroblasts, was found to be expressed by cells in all six Seurat clusters or all three subpopulations. Two genes likely involved in angiogenesis, Notch3 and Angpt2 (32), were found to be highly expressed in pericytes (Fig. 6B). Col10A1, which is frequently expressed by TAFs of many cancer types and by myofibroblasts during wound healing (33, 34), was expressed by a subset of smooth muscle-like cells in cluster 0 and some pericytes (Fig. 6B). The analysis in TAFs from HT-29 tumor resulted in the same cellular components and similar marker gene expression pattern.

Unexpectedly, we found that the EC gene Tek/Tie2, which has been used as an EC-specific Cre-driver (35) and has also been reported to be expressed in a tumor-associated macrophage/monocyte subpopulation (36–38), was found to be expressed in some fibroblast-like cells from COLO205 tumors (Fig. 6B). The detection rate of Tek/Tie2 in TAFs ranged from 3.4% to 11% in the fibroblast-like population, compared with 28% to 48% in ECs. In contrast, from the approximately 2,000 leukocytes captured from all COLO205 xenograft experiments, of which macrophages/monocytes made up approximately one-third of the population, Tek/Tie2 expression was only detected in a total of eight cells (approx. 1% of tumor monocytes/macrophages). In addition, the expression level of Tek/Tie2 was much lower in leukocytes (1 to 5 UMI) compared with TAFs (up to 30 UMI) and ECs (up to 100 UMI). In parallel, none of the 303 leukocytes isolated from HT-29 tumors expressed Tek/Tie2. These data suggest that either Tie2-expressing monocytes are a very rare population present only in some tumors, or that Tie2 expression in monocytes is lower than the detection limit in our current experiments.

Similar to tumors, cardiac fibroblasts, a collective term, can also be classified into three distinct cell types, which correspond with fibroblast-like cells, smooth muscle-like cells and pericytes (Supplementary Fig. S7A–S7C). We analyzed the major differences between tumor fibroblasts and heart fibroblasts by comparing gene expression profiles of the smooth muscle-like or fibroblast-like cells from COLO205 tumors to those from heart. With three-fold cutoff, smooth muscle-like cells from heart and tumor showed 279 differentially expressed genes. Similarly, fibroblast-like cells from heart and tumor showed 217 differentially expressed genes. Between the differentially expressed genes in the two cell types, there were 107 overlapping genes, which could be placed in several functional categories. Among them, more than two-thirds were upregulated genes and more than half of the upregulated genes produce secreted molecules (Supplementary Table S7).

Tumor stromal fibroblasts showed several distinct molecular features compared to similar cells from heart (Fig. 6C). First, both smooth muscle-like and fibroblast-like cells derived from tumors (TAF) expressed high levels of chemokines that were rarely seen in those derived from the heart. Second, genes encoding extracellular matrix proteins and TGF-ß pathway genes were highly upregulated in TAFs compared with stromal cells in heart. For example, Spp1 and Fn1 were almost exclusively expressed in TAFs. Third, genes involved in protein processing, transport and lysosome activity such as cathepsins were also upregulated in TAFs. Fourth, genes known as cell lineage differentiation markers were downregulated in TAFs as indicated by a lower expression level of smooth muscle-like–specific genes such as Acta2 and Tagln and the fibroblast-specific gene Dcn. Interestingly, Vegfd (Figf) was expressed at much higher level in the fibroblast-like compartment of tumors compared with the equivalent cells in heart.

Finally, we compared genes expressed in pericytes between tumor and heart. Among pericyte-defining genes, signaling molecules like Notch3 and Pdgfrb were expressed in both tumor and normal pericytes. However, the angiogenic factor Angpt2 was only detected in tumor pericytes (Fig. 6B; Supplementary Figs. S6A–S6C and S7A–S7C), consistent with its role in vascular remodeling and active angiogenesis. The relative ratio of cell populations differed in tumor and heart: in heart, isolated smooth muscle-like and fibroblast-like cells were found at a 1:4 cell count ratio, whereas in tumors the ratio was approximately 1:1. These data suggest an enhanced role of smooth muscle-like cells, especially pericytes, in the tumor microenvironment.

Because we observed effects on ECs in COLO205 tumors after treatment with aflibercept or aDll4, we examined differentially expressed genes in fibroblast subpopulations under drug-treated and control-treated conditions. Genes involved in transcription and translation were downregulated upon both antiangiogenic treatments (Supplementary Table S8). As expected, we observed a significant upregulation of several hypoxia-responsive genes, such as Ddit4, Angptl4, Pgk1, and Vegfa, in TAFs after treatment with aflibercept, aDll4, or combination over control hFc treatment at both 24 hours and 72 hours (39, 40).

Using expression levels of these hypoxia genes as a surrogate for hypoxia severity, we found that combination treatment led to the strongest hypoxic effect in TAFs, whereas the effect of aflibercept at the same time points was somewhat less. With all treatments, the hypoxia signature was more pronounced at 72 hours than at 24-hour posttreatment. Notably, in addition to Vegfa, Inhba and Il4a were also upregulated in smooth muscle-like cells, with fold changes correlating with hypoxia severity (Fig. 6D). Downregulation of Igfbp7 also occurred in TAFs, as was observed in ECs upon drug treatments (Fig. 6D and Supplementary Table S8). Thus, our data showed transcriptional changes in TAFs in response to the onset of hypoxia, and the induced genes including pro-angiogenic genes like Vegfa.

In this study, we used single-cell transcriptome profiling to characterize the heterogeneity of ECs and fibroblasts in tumors as well as in heart tissue. COLO205 xenograft tumors were chosen to build on our previous studies, where we used bulk RNA analysis to identify gene expression changes reflecting vascular response to VEGF blockage (8). From nearly 2,900 ECs collected from COLO205 xenograft tumors, cells could be classified into three subpopulations, which appear to correspond to previously described endothelial tip- and stalk-cells, as well as an intermediate we called transitional cells. Marker genes that can distinguish these endothelial subpopulations were identified; in addition to known markers such as Kcne3 and Dll4 for tip cells, and Vwf and Selp for stalk cells, many of the genes are novel markers, including Ramp3 and Ednrb for tip-like cells, and Ackr1 and Spint2 for stalk-like cells. Notably, over 80 percent of these marker genes were reproduced in another xenograft tumor model, HT-29.

In addition to tumor ECs, we characterized over 5,300 ECs from normal adult mouse heart, which provides an abundant source of quiescent arteriole, venous, and capillary cells. In subcutaneous tumor models used in our studies, the tumor vasculature can originate from several tissues surrounding the site of tumor implantation including dermis, subcutaneous space or even underlying muscle, making it difficult to utilize a single tissue source for comparison. Previous studies have documented that heart ECs resemble muscle ECs, but are distinct from other organs including liver ECs (17, 18). The distinctive features of heart ECs reported in other studies were not the basis of distinguishing heart from tumor ECs in our studies. Heart ECs were classified into four subpopulations, corresponding to capillary, arteriole, venule, and lymphatic ECs. It is worth mentioning that EphrinB2 (Efnb2) and its receptor Ephb4 have been previously reported as arteriole and venule markers, respectively, during mouse embryogenesis. However, Efnb2 was not well expressed by ECs in adult heart and Ephb4 was not venule-specific by our single-cell profiling or by histologic analyses using RNA-scope. It is possible that expression of these markers is tissue specific, or dependent on developmental stage, or that the expression levels are below what can be reliably detected with current single-cell transcriptome methods.

We also interrogated heterogeneity of other tumor stromal components, namely TAF cells, and documented several subpopulations, which appear to correspond to fibroblasts-like, smooth muscle-like cells and pericytes. In the past these tumor stromal cells have been difficult to clearly classify because of the upregulation of certain “lineage marker” genes such as aSMA, which are specific in quiescent normal tissues but expressed by several cell types in tumors. To overcome this cross-expression, we utilized a panel of genes to classify these cells.

We also examined changes in EC and fibroblast subpopulations following treatment with antiangiogenic agents that block VEGF and Dll4-Notch signaling. In tumors, treatment with aflibercept rapidly inhibits EC proliferation across all cell subpopulations, particularly tip-like cells. In contrast, treatment with aDll4 Ab promotes EC proliferation. These findings are in agreement with our previous studies (8, 9), where we observed changes in EC numbers upon aflibercept or aDll4 treatment. Measured by cell-cycle activity, the effect of aDll4 was more pronounced at 24 hours than at 72-hour posttreatment, whereas the antiangiogenic effect of aflibercept was prominent at both 24-hour and 72-hour posttreatment. Measured by the extent of gene changes, tip-like cells were affected the most by antiangiogenic treatments, with similar patterns of transcriptional program change upon aflibercept or aDll4 treatment. Several tip-like cell markers were upregulated under either drug treatment, suggesting that expression of many tip cell genes could be driven by the hypoxic tumor microenvironment rather than being drug specific. However, some drug-specific changes, such as Tgfb1 upregulation upon aflibercept treatment in tip-like cells, might be worth further investigation. In heart, aflibercept treatment had minimal impact on EC gene expression, whereas aDll4 Ab treatment promoted cell proliferation as well as tip cell properties specifically in the capillary subpopulation. Plvap, a key molecule for controlling vascular permeability and endothelial diaphragms (41), was highly expressed by all tumor ECs, and its expression was significantly downregulated after both aflibercept and aDll4 Ab treatment at both 24 hours and 72 hours. Pgf, a molecule that can enhance the sensitivity of VEGF receptors to their ligands was upregulated in ECs under hypoxic stress upon antiangiogenic drug treatments.

Treatment with antiangiogenic agents also changed gene expression patterns in TAFs. One of the expected gene signatures relates to hypoxic response. Several molecules involved in angiogenesis are produced by TAFs following antiangiogenic therapies. First, Vegfa was upregulated under all drug treatment conditions and the fold of change is correlated with hypoxia severity. Second, Angpt2, which is expressed in ECs and can act as a Tie2 agonist and limits the effect of VEGF inhibition (32), was found to be also expressed by tumor pericytes but not those in heart. Third, Vegfd was detected at high levels in tumor fibroblast-like cells. Thus, other stromal components can contribute to tumor angiogenesis, and potentially resistance to antiangiogenic drugs, by secreting feedback molecules upon local oxygen deprivation.

As an emerging technology, there are important limitations in current single-cell transcriptome methods. One key consideration is that genes expressed at lower levels are likely detected stochastically (42). This nonhomogeneous detection of some genes necessitates using a panel of genes to define subpopulations of cells. It also demands caution when trying to interpret EC heterogeneity. For example, previous studies have reported that Tek/Tie2 expression is heterogeneous in ECs (43), and thus obvious questions from our studies would be in which subpopulation(s) is Tek/Tie2 expressed and what genes are coexpressed in Tie2-positive ECs. Indeed, in our studies, Tek/Tie2 was detected in only a subset of ECs. However, Tek/Tie2 expression was detected at approximately equal frequencies in both tip or stalk ECs (Supplementary Fig. S4B), thus Tek/Tie2 does not appear to be exclusive to any subpopulation that could be delineated in these studies.

Taken together, our findings provide definitions on the heterogeneity of ECs and other stromal cells in tumors. Further, we provide evidence on the roles of VEGF and Dll4-Notch signaling in specifying tumor EC phenotype, and highlight the response of stromal cells to antiangiogenic therapies.

Y. Wei has ownership interest (including patents)in Regeneron Pharmaceuticals, Inc. No potential conflicts of interest were disclosed by the other authors.

Conception and design: A. Eichten, A.A. Parveen, G. Thurston

Development of methodology: A. Eichten, A.A. Parveen, M. Ni, Y. Wei

Acquisition of data (provided animals, acquired and managed patients, provided facilities, etc.): A. Eichten, A.A. Parveen, C. Adler, T. Nevins, M. Ni, Y. Wei

Analysis and interpretation of data (e.g., statistical analysis, biostatistics, computational analysis): Q. Zhao, A. Eichten, Y. Huang, W. Wang, Y. Ding, G. Thurston

Writing, review, and/or revision of the manuscript: Q. Zhao, A. Eichten, W. Wang, G. Thurston

Administrative, technical, or material support (i.e., reporting or organizing data, constructing databases): C. Adler, A. Adler, T. Nevins

Other (performed all the experiments to generate the data shown in the article): A.A. Parveen

We would like to thank Calvin Lin for his input and discussions regarding experimental design and single-cell data collection, Harish Adoni for technical assistance, and Olulanu Aina for helping to assess the RNAScope data.

The costs of publication of this article were defrayed in part by the payment of page charges. This article must therefore be hereby marked advertisement in accordance with 18 U.S.C. Section 1734 solely to indicate this fact.