Single cell analysis reveals intra‐tumour heterogeneity, microenvironment and potential diagnosis markers for clear cell renal cell carcinoma

We found specific gene sets and cell clusters that were correlated with the prognosis for clear cell renal cell carcinoma (ccRCC). We also discovered transcription factor regulons and ligand–receptor pairs that might become potential therapeutic targets.


Single cell analysis reveals intra-tumour heterogeneity, microenvironment and potential diagnosis markers for clear cell renal cell carcinoma
We found specific gene sets and cell clusters that were correlated with the prognosis for clear cell renal cell carcinoma (ccRCC). We also discovered transcription factor regulons and ligand-receptor pairs that might become potential therapeutic targets.
Single-cell RNA sequencing (scRNA-seq) has yielded insights into tumour origin or immune composition's effect on clinical outcome in ccRCC, 1,2 but was limited to characterize potential new targets.
Here, we performed scRNA-seq on seven patients' tumours and matched five normal samples (Table S1 and Figure 1A). A total of 37 243 cells were divided into six cell types according to marker genes ( Figure 1B). Potential malignant cells were distinguished via copy number variations (CNVs). 3 Epithelial cells having CNVs were named as 'changed' ( Figure 1B). The top five differentially expressed genes (DEGs) were shown ( Figure 1E). According to the cell distributions and sample origins ( Figure 1C), we found the heterogeneity of tumour epithelial cells was more obvious than other cells. We then used the Harmony algorithm to minimize the batch effect ( Figure 1D,F). Myeloid cells appeared more in tumours than in normal samples, indicating their active role in tumourigenesis ( Figure 1G). Scoring for genes with different biological functions (Table S2) indicated that hypoxia score was highest in tumour and endothelial cells ( Figure S1).
According to canonical makers (Table S3), epithelial cells were divided into six types ( Figure 2A): tumour cell, proximal tubule cell (PT), collecting duct cell (CD), principal cell (PC), podocytes and VCAM1 + PT. The ccRCC was thought to originate from VCAM1 + PT, 1 and our data confirmed the existence of this cluster of cells, and these cells were only found in normal samples ( Figure S2B). Based on CNVs, epithelial cells were divided into three subpopulations: tumour, normal and 'changed'. Normal cells were mostly PT, while 'changed' cells were comprised CD, PC and VCAM1 + PT ( Figure 2B,C and Figure S2A). To reveal the dynamics of transcriptional profiles of epithelial cells, we applied trajectory analysis ( Figure 2D,E). Normal cells are located at the right end of the pseudotime trajectory, while tumour cells are located at the left end, indicating a transition tendency. VCAM1 + PT occupied all the branches, suggesting their pluripotency ( Figure 2I). PAGA also confirmed VCAM1 + PT were most similar to tumour cells ( Figure 2F). Mon-ocle2 defined the branches of tumour cells as two 'states' (state 3 and state 4, Figure 2E). State 3 genes functioned in proliferation, hypoxia and angiogenesis, while state 4 genes functioned in stress and immune response ( Figure 2G,H and Figure S2C,D and Table S4). High expression of state 3 genes was related to improved overall survival in a TCGA cohort containing 531 patients 4 ( Figure 2J). NMF 5 found six tumour epithelial-specific gene programs ( Figure 2K). Highly expressed program 2 was correlated with favourable overall survival, while highly expressed program 4 (Table S5) was correlated with unfavourable survival ( Figure 2L,M). Program 2 functioned in interferon signalling and antigen processing ( Figure 2N), whereas program 4 functioned in stress response, growth and apoptosis ( Figure 2O). SCENIC 6 discovered transcription factor regulons differentially expressed in epithelial cells ( Figure 3A,B). NR1I3 regulon was normal cell specific, while IRX3 regulon was located mainly in tumour cells ( Figure 3C). Co-expressed regulons clustered together by correlation analysis ( Figure S3A). IRX3 level in the tumour was apparently higher than normal samples ( Figure 3D). IRX3 was correlated to ccRCC marker CA9 and EMT gene VIM ( Figure 3E and Figure S3B). High expression of IRX3 was correlated to unfavourable survival in the TCGA cohort ( Figure S3C). We also performed immunofluorescence staining for the seven patients and tissue microarrays containing 340 ccRCC samples. We found IRX3 staining in patient 1t ( Figure 3F) and tissue microarrays   Figure S3E). IRX3 played a role in ccRCC progression.
Myeloid cells had 10 clusters ( Figure 4A and Figure S4A). We generated a gene signature file of the 10 clusters and used CIBERSORTx 7 to infer cell-type proportions in the TCGA cohort. Higher infiltration of C4_C3AR1 + macrophages in the tumour was associated with worse survival, while higher infiltration of C3_CD163 + macrophages or C8_cDC2 was related to improved survival ( Figure 4B). C4_C3AR1 + cells expressed tissue-resident marker CD74 and CD81 and scavenger receptor MSR1, which was an anti-inflammatory M2 marker. M1 and M2 gene signatures (Table S6) scoring found that C1_IL1B + cells had higher M1 score, while C3_CD163 + , C5_MSR1 + cells had higher M2 scores ( Figure 4D). RNA velocity analysis 8 indicated M1 to M2 transition tendency exist ( Figure 4C), and M2-like cells were mainly from tumour samples based on the sample origins ( Figure S4B).
In subgroups of T cells ( Figure 4E), CD8 + T cells had three clusters: 4-1BB+, ex_ef (expressing exhausted makers and effector molecules) and cycling T cells. Higher infiltration of 4-1BB + T cells or NKT0 cells was related to worse survival, while higher infiltration of ex_ef T cells was related to better survival ( Figure 4F). The previous opinion that CD8 + T cells predicted bad outcomes 9 was not accurate, more markers were needed to evaluate the prognosis.
CellPhoneDB inferred interactions between different cells. The strongest interactions appeared between epithelial and myeloid cells ( Figure S4C,D and Table S7). Typical crosstalks were listed ( Figure 4G,I). Intriguingly, C3_C3AR1 crosstalk happened between tumour epithelial and myeloid cells. The exact roles of complement components C3 and C3AR1 in ccRCC have not been described. C3 and C3AR1 did co-localize in tumour samples (Figure 4H). LGALS9_CD47 interactions were also obvious ( Figure 4I), staining validated this crosstalk ( Figure 4J). CD99_PILRA between endothelial and myeloid cells (Figure 4I) suggested a decreasing immune infiltration effect. 10 LGALS9_HAVCR2/TIM3 also reminded a suppressive impact, suggesting endothelial cells' immune-suppression role.
Our studies revealed gene sets and specific macrophage and T cell clusters that can predict prognosis. Potential therapeutic targets such as IRX3 and inhibitory interactions among different cells provided insights for ccRCC therapies.