TECs from ccRCC tumors have a gene expression profile distinct from NECs isolated from renal cortex.
To characterize the expression phenotype of TECs from ccRCC and to compare to their normal counterpart (NEC) isolated from normal adjacent renal cortex, we performed scRNA-seq on TECs and NECs purified by flow cytometry from four different donor nephrectomies. The four TEC/NEC sample pairs were batched into two separate scRNA-seq libraries. We spiked in human umbilical vein endothelial cells (HUVECs) in all samples allowing us to study batch-to-batch variation between libraries. TECs/NECs, HUVECs and an aliquot of the entire tumor/normal adjacent tissue (NAT) were surface-labeled with unique hashtag oligonucleotide (HTO)-labeled antibodies targeting both human CD298 and β2 microglobulin prior to pooling. The oligonucleotide barcodes attached to the HTO-labeled antibodies remained associated with the cells during library generation. We utilized the readout of these HTO-antibodies as an independent measure of cellular identity at the analysis stage (see below).
Uniform Manifold Approximation and Projection (UMAP) projections showed distinct cell types clustered according to gene expression patterns (Fig. 1A). Using the HTO-barcodes as an orthogonal measure of cell identity, we saw clear separation of isolated ECs (HTO-antibody1) from non-EC cell populations (HTO-antibody2) within tumor/NAT digests and from HUVECs (HTO-antibody3) (Suppl. Figure 1A-C). In addition, The HTO-antibody3 labeled HUVEC populations were overlapped between different libraries indicating minimum batch effect, suggesting that there was no observed bias in representation of cells as a function of donor or sequencing library (Suppl. Figure 1B, D, E). Furthermore, expression of the ACTB housekeeping gene was uniform across the different clusters, consistent with comparable transcript counts/sample over the different cell types (Suppl. Figure 1F). Overall, after quality control and filtering, our dataset consisted of 14,982 total cells (library1: 3,630 cells and library2: 11,352 cells), of which, 7,512 cells were TECs or NECs (library1: 2,467 ECs and library2: 5,045 ECs).
Further exploration of the UMAP projection revealed distinct cell type clusters (ECs, immune cells, tumor cells, normal kidney epithelial cells, and the control HUVEC clusters, Fig. 1A) that expressed expected cell-defining markers (Suppl. Figure 1G) confirming the cell type assignment. For example, ECs expressed PECAM1 (CD31) and CDH5 (VE-Cadherin). Though HUVECs expressed EC markers, their overall expression pattern was distinct from primary NECs and TECs. Focused analysis of the primary EC cluster showed a clear separation of TECs and NECs as well as a transitional population with an intermediate phenotype (Fig. 1B). Pseudotime trajectory analysis indicated a possible developmental trajectory from NECs towards TECs (Fig. 1C). We then classified the ECs according to known phenotypes described as a function of their anatomic location (32) (Table S2). This analysis revealed that NECs from the kidney have a heterogeneous mixture of venous, arterial, capillary, and lymphatic vasculature phenotypes (33–35), whereas TECs predominantly exhibited a venous EC phenotype (Fig. 1D). Cell type classification using a set of angiogenic EC markers (GO angiogenic pathway 0001525 (36, 37)) also revealed that TECs were more likely to exhibit an angiogenic phenotype than NECs (Fig. 1E and F). Taken together, these analyses confirmed that TECs have a phenotype that is distinct from NECs.
To explore the basis of this distinct TEC phenotype, we performed differentially expressed gene (DEG) analysis. A total of 1,636 genes were differentially expressed when comparing TECs with NECs (Table S3). The top 3 DEGs highly expressed in TECs versus NECs were IGFBP3, ENPP2 and INSR, whereas MT2A, PLAT and IFGBP5 were the top 3 genes preferentially expressed in NECs (Fig. 1G). Pathway analysis of the DEGs identified extracellular matrix (ECM) organization, laminin interactions, collagen fibrils assembly, non-integrin membrane-ECM interactions, integrin cell surface interactions and IGF regulation as functional pathways upregulated in TECs (Fig. 1H). Exploration of individual matrix metalloproteinase (MMPs) and a disintegrin and metalloproteinase with thrombospondin motifs (ADAMTS) family members also identified several genes with increased expression in TECs (Suppl. Figure 2A & B). These included enzymes with known roles in promoting angiogenesis such as MMP9 (38), MMP14 (39) and ADAMTS13 (40). Conversely, in comparison to NECs, TECs had reduced expression of pathways involved in immune regulation such as signaling by interleukins and interferon gamma signaling (Fig. 1I).
Identification of IGFBP3 as a marker gene for TECs
UMAP projections clearly showed preferential expression of the top 3 DEGs in TECs (IGFBP3, ENPP2 and INSR) and NECs (IGFBP5, PLAT and MT2A) (Fig. 2A). Previous studies have proposed IGFBP3, ACKR1 and PLVAP as universal TEC marker genes (15). While IGFBP3 was clearly expressed in TECs, we found that ACKR1 was instead expressed in a transitional subpopulation of ECs derived mainly from normal kidney tissue with a smaller contribution of ECs isolated from the tumor. PLVAP had higher expression in TECs, but NECs also expressed this gene and therefore it did not appear to be a TEC-specific marker. The cell-adhesion molecule VCAM1 has reported expression on a subset of pancreatic TECs (41), though in our dataset, it was localized to the same transitional subpopulation of ECs that expressed ACKR1. These UMAP findings confirmed the results of the DEG analysis but also identified differences from previous studies which could be due to the increased resolution of our enriched EC dataset and ccRCC-specific patterns.
Given the consistency of IGFBP3 as a specific marker of TECs in our dataset and in previous studies (15), we validated its anatomic localization as well as that of its paralog IGFBP5, using RNA in situ hybridization (RNA-ISH) on formalin-fixed paraffin-embedded (FFPE) sections of RCC and NAT tissues from our 4 donors. IGFBP3 and IGFBP5 were two of the highest DEGs between TECs and NECs (Fig. 2B, Table S3). Correlating with the scRNA-seq data, RNA-ISH for IGFBP3 showed weak expression in rare peritubular capillaries and some tubules, whereas it showed strong expression in the RCC vasculature and patchy expression in the tumor cells themselves (Fig. 2C). By contrast, IGFBP5 showed strong staining in glomerular and peritubular capillary ECs in normal kidney tissues as well as some arterial smooth muscle cells (not shown). IGFBP5 showed only minimal and patchy expression in RCC tissues. When we examined all the cell populations in our dataset, we confirmed minimal IGFBP3 expression in RCC tumor cells as well as IGFBP5 expression in pericytes, which have an overlapping gene expression signature with arterial smooth muscle cells (Suppl. Figure 3A-C). Interestingly, HUVECs were discordant with both TEC and NEC phenotypes demonstrating low levels of both IGFBP3 and IGFBP5 expression. Having validated IGFBP3 and IGBFP5 as markers of RCC-associated TECs and NECs respectively, we examined the expression profiles of these genes in tumor and normal tissues in The Cancer Genome Atlas (TCGA) dataset. Consistent with our data, IGFBP3 showed the highest expression in RCC tumors, whereas one of the highest sites of IGFBP5 expression was normal kidney (Suppl. Figure 3D, E). Increased expression of IGFBP3 in RCC (and also in lung and colorectal cancer) was associated with inferior overall survival, while IGFBP5 expression had no apparent impact on prognosis in RCC (Suppl. Figure 3F, G).
Integration of published RCC scRNA-seq data highlights the advantage of EC enrichment
Next, we compared the expression data from this study with previously published scRNA-seq datasets generated from ccRCC. We identified 10 published studies with ccRCC scRNA-seq data, of which only 5 included cells likely to be ECs that were defined by measurable coexpression of PECAM1 and CDH5. None of the prior studies specifically enriched for ECs resulting in an average of 217 ECs/sample (range 2-518/cells sample, Table S4). By contrast, the EC count/sample in this study was 3.7-fold higher than the average published dataset at 939 cells/sample, thus representing the deepest reference dataset for RCC ECs. Of the 5 published studies with EC expression data, 2 did not include paired NAT along with the RCC samples (GSE152938 and GSE171306) and study GSE139555 had only 2 ECs per sample on average. The remaining two published datasets from Li et. al. (10.17632) (42) and Zhang et. al. (GSE159115) (12) included paired TECs and NECs and had adequate EC sampling (207 and 245 cells/sample respectively). We co-clustered these two datasets together with our data to generate a unified UMAP projection (Fig. 3A). 98.4% of all ECs from our samples belonged to cell clusters shared with the published datasets (hashed outline). By contrast, only 66.0% of cells from the 2 published studies could be assigned to the composite cell cluster. The non-overlapping cells showed distinct UMAP localizations even between the two published studies, which may be due to methodologic differences, pre-analytic variation, and/or batch effects. Focusing on the shared cluster of ECs across all 3 datasets, we analyzed DEGs in TECs compared to their matched NECs (Fig. 3B). This identified 505 DEGs that were shared across all 3 datasets. Pathway analyses of this gene list identified significant alterations in ECM organization, neutrophil and platelet degranulation, regulation of IGF transport and uptake by IGFBPs and others (Fig. 3C). Assignment of TECs and NECs to the shared UMAP projection showed representation of almost all EC clusters identified in our samples and in the published studies, though there were some minor differences in represented cell populations (Fig. 3D). Of note, the reciprocal expression of IGFBP3 and IGFBP5 in TECs and NECs was reproduced in the combined data (Fig. 3E). This analysis demonstrates that the data from our samples recapitulates major EC clusters seen in previous RCC studies, but also provides a much deeper profiling of matched TEC and NEC expression signatures. The advantage of deeper profiling of isolated ECs is evident since IGBFP5 was not identified as a marker of NEC in either of the two previous studies and IGFBP3 was not identified as a TEC marker in one of them (12, 42).
Congruence of TEC phenotypes between kidney and liver cancer
It is not clear if TECs from different types of tumors adopt a unique organ-specific phenotype or a common phenotype that is necessary for tumor growth regardless of tumor location. The second possibility is intriguing since it could lead to therapeutic strategies that target TECs across tumor types. Supporting this concept, a recent study suggested that tip cells from TECs across tumor types adopt a congruent phenotype (43). The deep profiling of ECs enabled by our dataset allowed us to ask a similar question in comparing ccRCC TECs to those from a recently published study on hepatocellular carcinoma (HCC) (20). As we did previously for the RCC scRNA-seq datasets, we combined the scRNA-seq data from hepatocellular carcinoma with our own data to create a unified UMAP projection (Fig. 4A). 92.4% of our RCC TECs overlapped with the HCC TECs and conversely, 75.8% of TECs from the HCC dataset overlapped with TEC clusters also present in this RCC study. While there were TEC populations that were unique to both cancer types, this analysis confirmed that the majority of TECs from both kidney and liver tumors were clustered together and exhibited a congruent expression phenotype (Fig. 4A). By contrast, only 37.3% of NECs from normal kidney overlapped with NEC clusters in liver and reciprocally, 51.7% of NECs from normal liver overlapped with NEC clusters in kidney. This finding suggests that normal tissues have NECs with distinct and organ-specific phenotypes. We then asked which genes were differentially expressed between TECs and NECs in both liver and kidney cancer (Fig. 4B). 561 DEGs were shared between both cancer types and these genes were enriched for many of the same pathways previously seen in RCC TECs alone (Fig. 4C). Mapping onto the combined UMAP projection identified unique NEC and TEC populations in both tissues (Fig. 4D). IGFBP3 expression remained a marker of both kidney cancer and liver TECs, however IGFBP5 expression was mainly seen in the kidney and not in liver NECs (Fig. 4E). These analyses showed that while NECs from different organs retained distinct phenotypes that may be important for that organ’s unique functions, by contrast, the majority of TECs from two different tumor types shared a congruent gene expression profile with upregulated pathways related to ECM organization, neutrophil/platelet degranulation, IGF transport and uptake by IGFBPs and others.
TECs and NECs maintain their distinct gene expression patterns in ex vivo cell culture
A prerequisite to understanding the functional properties of unique cell types is to establish an in vitro culture system. Few previous studies of TECs have examined their properties in vitro, due to their slow growth and difficulty in maintaining viability. We were able to isolate primary NECs and TECs from 5 donors for successful in vitro culture for up to 2–3 passages in human VEGF-containing media (34, 44). TECs were larger, stellate, and adopted a lightly overlapping cobblestoned growth pattern in vitro while NECs were smaller, spindle shaped, and adopted a whorled arrangement in culture (Fig. 5A). We then generated bulk RNA-seq data from ECs at passage 3 as well as primary TECs and NECs freshly isolated (not cultured) from the same donors. We detected 1,002 DEGs between cultured TECs and NECs compared to 3,028 DEGs between primary isolated TECs and NECs. Comparison of the overlap of the DEGs between TECs and NECs showed that, although culturing ECs in vitro reduced the overall number of DEGs, 783 (78.1%) of cultured EC DEGs were still shared with the primary purified ECs (Fig. 5B).
To explore the core programs retained by TECs and NECs in culture, we performed pathway enrichment analysis of the shared DEGs. The major upregulated pathways that were conserved in vitro related to neutrophil degranulation, ECM organization, and regulation of interactions with immune/lymphoid cells (Fig. 5C). The major downregulated pathways preserved in TECs in vitro related to GPCR ligand binding, ERBB4 and ERBB2 signaling (Fig. 5D). Next, using quantitative RT-PCR, we validated our RNA-seq findings using key TEC and NEC marker genes previously identified (Fig. 5E). These experiments demonstrated that primary cultures of TECs and NECs could be established in vitro and retain key expression programs of their freshly isolated counterparts.
TECs are resistant to VEGF withdrawal and exhibit distinct binding preferences for CD45+ leukocytes
Having established conditions for primary culture of TECs and NECs, we next explored their phenotypic and functional differences in vitro. Previous studies demonstrated that the survival of NECs from kidney in culture requires supplementation with VEGF in the culture medium (45). We first assessed whether TECs and NECs exhibited differential sensitivity to this important trophic factor by removing it from the culture. Within 48 hours of VEGF removal, NECs from 4 donors showed a sharp reduction in viability as assessed by DAPI flow staining (Fig. 6A). In contrast, TECs from those same 4 donors showed no loss in viability under the same conditions of VEGF withdrawal. We also performed live cell imaging of TECs and NECs over 48 hours (Fig. 6B). TECs showed a comparable increase in % confluency over the 48-hour assay period with or without the presence of VEGF in the culture medium. In contrast, the % confluency of NECs increased much more slowly without VEGF compared to the presence of VEGF. This differential sensitivity of TECs and NECs to VEGF may be related to the higher expression of KDR and other receptors for autocrine or serum-derived trophic factors both in the isolated and ex vivo cultured ECs (Suppl. Figure 4).
TECs represent the first barrier to entry of immune cells into the tumor mass and are an integral component of the tumor microenvironment. Immune regulatory pathways were altered in our analysis of DEGs both in freshly isolated (Fig. 1) and cultured TECs (Fig. 5). This observation led us to empirically test whether CD45+ leukocytes exhibited different adherence to TECs versus NECs. We isolated CD45+ leukocytes from the tumor mass, NAT or peripheral blood mononuclear cells (PBMCs). We applied fluorescent labeled leukocytes to autologous labeled TECs or NECs in the IncuCyte imaging platform (Figure. 6C). After 24 hours of co-culture, we enumerated bound CD45+ leukocytes and normalized that count versus the number of TECs or NECs in the culture (Fig. 6D). Regardless of their origin (tumor, NAT, PBMC), more CD45+ leukocytes adhered to TECs compared to NECs (Fig. 6E). We then used flow cytometry to determine if specific CD45+ leukocyte subsets were responsible for this difference (Suppl. Figure 6A). We determined the relative fold binding of specific leukocyte subsets as a function of their origin (Fig. 6F, Suppl. Figure 6B). Regardless of their tissue of origin, T-cells and monocytes but not B cells and NK cells were more likely to interact with TECs. Although they were few in number, CD4+CD25+ T regulatory cells were much more likely to adhere to TECs versus NECs. These results may provide an explanation for the preferential recruitment of T-cells and monocytes into the RCC tumor microenvironment that has been reported previously (46). Taken together these studies demonstrated the utility of in vitro culture of TECs and NECs to study their distinct functional and immunological properties.