A novel stemness-related lncRNA signature predicts prognosis, immune infiltration and drug sensitivity of clear cell renal cell carcinoma

Abstract

Background

Clear cell renal cell carcinoma (ccRCC) is a prevalent urogenital malignancy characterized by heterogeneous patterns. Stemness is a pivotal factor in tumor progression, recurrence, and metastasis. Nevertheless, the impact of stemness-related long non-coding RNAs (SRlncRNAs) on the prognosis of ccRCC remains elusive. In this study, we aimed to delve into the SRlncRNAs of ccRCC and develop a signature for risk stratification and prognosis prediction.

Method

Gene-expression and clinical data were downloaded from The Cancer Genome Atlas (TCGA) and Gene Expression Omnibus (GEO) databases. We calculated RNA stemness scores (RNAss) for the samples to evaluate their stemness. SRlncRNAs and stemness-related mRNAs (SRmRNAs) in ccRCC were identified through weighted correlation network analysis (WGCNA), which employed sophisticated statistical methodologies to identify interconnected modules of related genes. Enrichment analysis was performed to explore the potential functions of SRmRNAs. Multiple machine learning algorithms were employed to construct a prognostic signature. Samples from TCGA-KIRC and GSE29609 cohorts were designated as the training and validation cohorts, respectively. Based on their risk scores, samples were stratified into low- and high-risk groups. Prognosis analysis, immune infiltration assessment, drug sensitivity prediction, mutation landscape, and gene set enrichment analysis (GSEA) were conducted to investigate the distinct characteristics of the low- and high-risk groups. Additionally, a web-based calculator was developed to facilitate clinical application. Expression and effects of SRlncRNAs in ccRCC were further corroborated through the utilization of single-cell RNA-seq (scRNA-seq), as well as in vitro and in vivo experiments.

Results

SRlncRNAs and SRmRNAs were identified based on RNAss and WGCNA. The least absolute shrinkage and selection operator (LASSO) in combination with multivariate Cox regression was selected as the optimal approach. Six SRlncRNAs were used to construct the prognostic signature. Samples in the low- and high-risk groups exhibited distinct characteristics in terms of prognosis, GSEA pathways, immune infiltration profiles, drug sensitivity, and mutation status. A nomogram and a web-based calculator were developed to facilitate the clinical application of the model. ScRNA-seq and RT-qPCR demonstrated the differential expression of SRlncRNAs between ccRCC tumors and normal tissues. In vitro and in vivo experiments demonstrated that downregulation of EMX2OS and LINC00944 affected the proliferation, migration, invasion, apoptosis, and metastasis of ccRCC cells.

Conclusion

We uncovered the crucial associations between SRlncRNAs and the prognosis of ccRCC. By leveraging these findings, we developed a novel SRlncRNA-related signature and a user-friendly web calculator. This signature holds great potential in facilitating risk stratification and guiding tailored treatment strategies for ccRCC patients. Both in vitro and in vivo experiments confirmed the role of SRlncRNAs in the progression of ccRCC.

Introduction

Renal cell carcinoma (RCC) is a highly malignant urogenital cancer, accounting for approximately 180,000 deaths globally in 2020 [1]. Among the various types of RCC, clear cell renal cell carcinoma (ccRCC) poses the greatest risk of metastasis compared to papillary and chromophobe renal tumor [2]. For patients with localized ccRCC, radical or partial nephrectomy is considered the primary treatment option [3]. Furthermore, clinical trials have demonstrated that immune checkpoint inhibitors (ICIs) and tyrosine kinase inhibitors (TKIs) can enhance the prognosis of patients with advanced ccRCC [4]. Several clinical and pathological factors, such as TNM stage, tumor necrosis, and microscopic vascular invasion, have been established as being correlated with the prognosis of ccRCC [5]. Presently, numerous prognostic models have emerged, with validated prediction models aiding in refining risk stratification and facilitating the development of individualized treatment plans and follow-up protocols [5].

Stemness, a defining characteristic of cancer cells, plays a pivotal role in tumor recurrence, metastasis, resistance to adjuvant therapy, and metabolic reprogramming [6]. However, the impact of stemness on the prognosis of ccRCC remains unclear. Long non-coding RNAs (lncRNAs), transcripts exceeding 200 nucleotides in length with limited protein-coding potential, can regulate gene expression through cis- or trans-regulatory mechanisms [7]. Prior research has shown that lncRNAs are involved in the proliferation, apoptosis, epithelial-mesenchymal transition, and metastasis of ccRCC [8,9,10,11]. LINC00887 could enhance the stem-like properties of by recruiting SPI1 and activating CD70 transcription [12]. LncRNA HOTAIR, in conjunction with the androgen receptor, synergistically activated the Hedgehog pathway and promoted the stemness of ccRCC cells [13]. Previous studies have identified stemness-related lncRNA (SRlncRNA) signatures in breast cancer and gastric cancer [14, 15]. Despite this growing understanding, the correlations between the expression of SRlncRNAs and ccRCC remain elusive. To address this gap, we aimed to clarify the roles of SRlncRNAs in ccRCC progression and develop an SRlncRNA-related signature to predict the prognosis of ccRCC patients. We developed a web-based calculator to facilitate the clinical application of this signature (https://dupeng1587.shinyapps.io/KIRC_SRlncRNAs_prognosis/).

Materials and methods

Data collection

In this study, gene expression and clinical data were obtained from the Cancer Genome Atlas (TCGA, https://portal.gdc.cancer.gov/) and Gene Expression Omnibus (GEO, https://www.ncbi.nlm.nih.gov/) databases. Our inclusion criteria comprised patients diagnosed with ccRCC and possessing complete mRNA and lncRNA data. Patients with incomplete clinicopathological information were excluded. The training cohort consisted of TCGA-KIRC, encompassing 72 normal samples and 539 tumor samples. Clinicopathological details for each patient were extracted. Additionally, clinical characteristics and genome expression data from 39 ccRCC patients were sourced from GSE29609 dataset [16]. GSE29609 was served as the validation cohort.

RNA stemness scores calculated

RNA stemness score (RNAss) referred to a quantitative metric used to assess the level of stemness associated with RNA molecules. It is indicative of the potential to maintain or influence cellular properties reminiscent of stem cells, such as proliferation, pluripotency, or resistance to differentiation. We downloaded the syn3156503 dataset from the Progenitor Cell Biology Consortium (PCBC, [17]. Utilizing the “gelnet” package, we employed the One-class Logistic Regression (OCLR) machine-learning algorithm to determine the weights of genes and construct a prediction model [18]. Subsequently, RNAss for the TCGA-KIRC samples were computed based on this model. The higher the RNAss, the greater the stemness of samples.

Identification of stemness-related mRNAs and lncRNAs

Weighted correlation network analysis (WGCNA) was a systematic framework designed for constructing gene co-expression networks from high-throughput genomic data [19]. In WGCNA, genes are considered as nodes, and the connections between them represent the strength of their co-expression patterns across different samples. In our study, it was employed to identify stemness-related mRNAs (SRmRNAs) and SRlncRNAs [19]. Using the “WGCNA” package, we generated co-expression networks for both mRNAs and lncRNAs in TCGA-KIRC. An appropriate soft threshold β was determined to satisfy the scale-free network criteria. Weighted adjacency matrices and topological overlap matrices (TOM) were constructed to establish the co-expression networks. Dynamic tree cutting was performed to identify relevant modules, with the target modules being recognized as those most significantly correlated with the RNAss. mRNAs and lncRNAs within the modules exhibiting the highest correlation were defined as SRmRNAs and SRlncRNAs, respectively, and selected for further investigation. After acquiring the SRmRNAs and SRlncRNAs, we explored the differentially expressed genes (DEGs) between carcinoma and adjacent tissue, with the statistical significance threshold at an absolute fold change (FC) > 1.2 and a false discovery rate (FDR) < 0.05.

GO enrichment, KEGG pathway and GSEA enrichment

To gain insights into the potential functions of DEGs, Gene Ontology (GO) analyses were utilized to understand their cellular components (CC), molecular functions (MF), and biological processes (BP) [20]. Furthermore, Kyoto Encyclopedia of Genes and Genomes (KEGG) and Gene Set Enrichment Analysis (GSEA) were performed to explore the pathways in which these genes may participate [21, 22]. These analyses (GO, KEGG, and GSEA) were conducted using the Database for Annotation, Visualization, and Integrated Discovery (DAVID 6.8, https://david.ncifcrf.gov/) and the Molecular Signatures Database (MSigDB v7.5.1, http://www.gsea-msigdb.org/gsea/msigdb/index.jsp). We sorted the top five pathways based on ascending p-value and presented them accordingly.

Establishment and validation of the signature

TCGA-KIRC and GSE29609 were utilized as the training cohort and external validation, respectively. We established overall survival (OS) as the primary endpoint of the study, defined as the interval from diagnosis to death. For patients who were alive, the endpoint was recognized as the last follow-up date. The expression of SRlncRNAs, OS time, and status were integrated using 10 machine learning algorithms, totaling 101 algorithm combinations. The algorithms employed in our study encompassed random survival forest (RSF), survival support vector machine (survival-SVM), ridge regression, elastic net (Enet), the least absolute shrinkage and selection operator (LASSO), supervised principal components (SuperPC), generalized boosted regression modeling (GBM), partial least squares regression for Cox (plsRcox), CoxBoost, and stepwise Cox regression. Five of these algorithms featured a variable filtering function: RSF, LASSO, CoxBoost, backward stepwise Cox, and both versions of stepwise Cox regression. One hundred and one algorithm combinations were applied to fit predictive models based on the leave-one-out cross-validation (LOOCV) framework. The Harrell's concordance index (C-index) was calculated for both datasets separately. The model with the highest average C-index was deemed optimal. The risk score for each case was calculated using this optimal model. The median of all sample scores was set as the cut-off value. Patients were then stratified into low- and high-risk groups based on their risk scores. Time-dependent receiver operating characteristic (ROC) curves were plotted, and the area under the curves (AUCs) at 1, 3, and 5 years were calculated to evaluate the accuracy of the established model. To assess the predictive performance of our model more accurately, we retrospectively collected previous published signatures related to ccRCC. Risk scores of all samples in TCGA-KIRC and GSE29609 cohorts were calculated by the published formula and coefficient of genes. And C-index value of these signature could be obtained.

Co-expression network

The co-occurrences of SRmRNAs and SRlncRNAs were investigated. Given the vast quantity of SRmRNAs, we intersected the SRmRNAs identified in our study with those from a previous study [14]. A co-expression network of SRlncRNAs and SRmRNAs was then constructed based on the correlations between them. The criterion for inclusion in the network was set at an absolute Pearson correlation coefficient greater than 0.4 and a p-value less than 0.00,001.

Construction of a ceRNA network based on SRlncRNAs and SRmRNAs

A Competitive Endogenous RNA (ceRNA) is defined as an RNA molecule, encompassing but not limited to lncRNAs, that operates by competing with other RNAs for shared miRNA binding sites. This competition modulates the accessibility of miRNAs to their canonical target mRNAs, subsequently influencing the post-transcriptional expression of those target genes. In our study, we constructed a ceRNA network centered on SRlncRNAs and SRmRNAs, both of which possess the capability to bind a common miRNA. We utilized the "multiMIR" package to explore potential mRNA-miRNA interactions. These interactions were rigorously validated across three databases: miRecords (http://mirecords.biolead.org/), miRTarBase (http://mirtarbase.mbc.nctu.edu.tw/index.html), and TarBase (http://microrna.gr/tarbase/). Specifically, we selected only those interactions that had been experimentally confirmed using a dual-luciferase reporter gene system. To delve deeper into potential lncRNA-miRNA interactions, we consulted the ENCORI database (https://rnasysu.com/encori/). Subsequently, we identified the miRNAs shared between these two sets of validated interactions. Ultimately, based on these findings, a ceRNA network was constructed.

Immune infiltration and immunotherapy

The Tumor Immune Estimation Resource (TIMER, [23], Estimating Relative Subsets of RNA Transcripts (CIBERSORT, [24], XCELL ([25], MCP-counter (available at [26], quanTIseq ([27], and EPIC ([28] algorithms were employed to compute the abundance of tumor-infiltrating immune cells across all samples. Drawing from a previous study, we performed single-sample gene set enrichment analysis (ssGSEA) on 29 immune signatures, encompassing 16 immune cell types and 13 immune function scores [29]. Subsequently, we analyzed the differences in tumor-infiltrating immune cell populations and immune function scores between low- and high-risk groups. We aimed to shed light on the activation status of immune cells within tumor tissues through the aforementioned analysis [30].

We predicted the response of samples to immunotherapy through the Tumor Immune Dysfunction and Exclusion (TIDE, http://tide.dfci.harvard.edu/) website [31]. The Submap algorithm was applied in a previously published immunotherapy cohort to further assess the accuracy of risk stratification in predicting the efficacy of immunotherapy [32,33,34]. By reviewing the CheckMate cohort, we analyzed the relationship between the expression of SRlncRNAs and the prognosis of ccRCC patients undergoing immunotherapy [35].

Drug sensitivity prediction

The drug sensitivity of the samples was predicted using the “pRRophetic” package [36]. We compared the differences in the 50% inhibitory concentration (IC50) values of the candidate drugs between the low- and high-risk groups, with the understanding that a lower IC50 value indicates higher drug sensitivity. DEGs between the low- and high-risk groups were identified. To further explore potential therapeutic drugs for KIRC based on these DEGs, we utilized the Connectivity Map (CMap, [37], L1000 Fireworks Display (L1000FWD, [38], and the Drug Gene Interaction Database (DGIdb, [39]. The intersection of potential drugs identified from these sources was obtained, and these drugs were considered as SRlncRNA-related therapeutic options. The structures of these potential drugs were retrieved from the DrugBank database [40].

Mutation analysis

Somatic mutation data was downloaded using the "TCGAmutations" package. Tumor Mutation Burden (TMB) was defined as the number of mutations detected per million bases. Genetic variants could play a significant role in the progression of cancer [41]. Our primary focus was on the differences in mutation status, frequencies, and TMB between the low- and high-risk groups.

Construction of nomogram

Each SRlncRNA was categorized as either protective or risk-associated. Both univariate and multivariate Cox proportional hazards analyses were conducted to identify clinical factors and the risk score that are associated with survival time. Based on the results, a nomogram was constructed, integrating the prognostic signatures. This nomogram can be utilized to predict the 1-, 3-, and 5-year OS rates for ccRCC. To facilitate clinical application of this model, we aimed to develop a web-based calculator.

Landscape of scRNA-Seq and spatial transcriptome data

We obtained a Single-cell RNA sequencing (scRNA-seq) dataset, GSE159115, from the GEO database [42]. This dataset comprised 27,669 cells from 7 KIRC patients, including 7 tumor samples and 6 normal samples. The scRNA-seq data of KIRC was independently analyzed using R (version 4.3.0) and Python (version 3.12). The "Seurat" package was utilized for quality control and further analysis. Genes expressed in fewer than three cells and cells expressing less than 200 genes were excluded. The "PercentageFeatureSet" function was applied to calculate the percentage of mitochondrial genes, and cells with a mitochondrial proportion exceeding 20% were removed due to their low quality or extensive mitochondrial contamination. Total transcripts per single cell were normalized to 10,000, followed by log transformation. Using the "FindVariableFeatures" function, we identified the top 2000 features for each single-cell dataset, which varied significantly among cells. Subsequently, the "ScaleData" function was employed to scale the expression levels of genes. Principal Component Analysis (PCA) was performed using the "RunPCA" function to achieve linear dimensional reduction. The "JackStrawPlot" and "ElbowPlot" functions were utilized to determine the optimal number of principal components (PCs) for downstream analyses. Cells were further clustered using the "FindClusters" function. The "RunUMAP" function was applied for dimensionality reduction and visualization of clusters. The "FindAllMarkers" function was used to identify specific biomarkers of different clusters. We set the parameters "only.pos = TRUE" and "min.pct = 0.25" to select highly expressed markers. Clusters in the single-cell dataset were annotated based on known marker genes, corresponding original articles, and the Tumor Immune Single-cell Hub 2 (TISCH2, http://tisch.comp-genomics.org/gallery/?cancer=KIRC&species =) database [43].

We thoroughly investigated the correlation between SRlncRNAs and the progression of ccRCC at the single-cell level, employing the following four strategies. Firstly, we conducted a comparative analysis of SRlncRNA expression between normal epithelial cells and malignant cells at the single-cell resolution. Based on these expression levels, cells were categorized into high-expression and low-expression groups for SRlncRNAs. Secondly, pseudotime analysis aimed to order cells along a continuous trajectory that reflects their biological progression. We used “monocle” package to simulate the trajectory of different cell states and differentiation paths. We employed pseudotime analysis to examine the differentiation disparities between cells with high and low SRlncRNA expression. Thirdly, we utilized the “infercnvpy” function to calculate Copy Number Variations (CNV) values for diverse cell types at the single-cell level, subsequently analyzing the differences in CNV values between cells with high and low SRlncRNA expression. Lastly, we leveraged the "CytoTRACE" package to quantify cell differentiation potential. Through this, we compared the variations in CytoTRACE scores across different cell populations.

Spatial transcriptome sequencing data from a ccRCC patient (PD45816, female, stage IIIa, and grade 4) was downloaded from a previous study [44]. Cell clusters and annotation references were based on the cellxgene tool (https://cellxgene.cziscience.com/). The H&E Stained images showed the interface between normal and tumor tissues. We mainly concentrated on the spatial expression and distribution of SRlncRNAs. These analyses were conducted by “anndata” package (python 3.12).

Cell culture

The Caki-1, 786-O, 769-P renal cancer cell lines, along with the HK-2 renal tubular epithelial cell line, were procured from MeisenCTCC (Hangzhou, China). These cell lines were cultured in RPMI-1640 medium (Gibco, NY, USA), supplemented with 10% fetal bovine serum (Gibco, NY, USA), and maintained in a 5% CO2 atmosphere.

Human tissue specimens

One hundred and two surgically resected ccRCC tissues were collected from Peking University Cancer Hospital. All specimens were frozen in liquid nitrogen. Our study was approved by the Ethics Committee of Peking University Cancer Hospital on July 16th, 2020 (Approval No. 2020KT85). All patients provided informed consent. The inclusion criteria were as follows: patients must not have received systemic therapy prior to surgery, be 18 years old or older, and have normal cardiac, renal, and liver functions. Patients with distant metastasis were excluded from this study.

Quantitative real-time PCR (RT-qPCR)

Total RNA was extracted from cell lines and ccRCC tissues using RNAiso Plus 9108 (Takara, Japan). The quantity and quality of the RNA were assessed using an ultraviolet spectrophotometer (Thermo, USA). Subsequently, the total RNA was reverse transcribed into cDNA using the PrimeScript RT Reagent Kit with gDNA Eraser (Takara, Japan). RT-qPCR was performed on a Roche LightCycler 480 system (Roche, Switzerland). GAPDH was used as a reference gene to normalize the RT-qPCR results. The primer sequences for SRlncRNAs and GAPDH are presented in Supplementary Table 1. The expression levels of SRlncRNAs were calculated using the 2^−ΔΔCt method, with GAPDH serving as the normalizing factor. The amplification protocol was as follows: initial denaturation at 95 °C for 30 s, followed by 40 cycles of PCR at 95 °C for 5 s and 60 °C for 30 s, and a final annealing step at 95 °C for 5 s and 65 °C for 1 min.

Western blotting

Cells were lysed with radioimmunoprecipitation assay buffer (RIPA, Beyotime Biotechnology, China) containing PMSF on ice. After centrifugation, the supernatant was collected. Protein concentration was measured using a BCA protein assay kit (Applygen Technologies, China). Protein samples (50 μg) were separated by sodium dodecyl sulfate–polyacrylamide gel electrophoresis (SDS-PAGE, 8%) and transferred onto polyvinylidene fluoride (PVDF) membranes (Millipore, USA). Following blocking with 5% skim milk, the membranes were incubated with primary antibodies at 4 °C overnight. The membranes were then incubated with horseradish peroxidase-conjugated secondary antibodies. Images were captured using a ChemiDocTM CRS + Molecular Imager (Bio-Rad, Hercules, CA, USA). GAPDH was used as an internal reference protein. The antibodies used in this study in Supplementary Table 2.

Organoid culture and growth assay

Resected tumor specimens from a ccRCC patient were cut into small fragments of 0.5 mm, washed using cold basal culture medium, filtered through a 100 µm pore mesh, and then centrifuged at 1500 rpm for 3 min to collect the samples. Then they were digested in DMEM (sourced from Gibco, NY, USA) enriched with IV collagenase (Sigma, Germany) at a temperature of 37 °C for a duration of 15 min. To halt the digestion process, PBS was introduced.

Following centrifugation and filtration, the tissue was then embedded within Cultrex Reduced Growth Factor Basement Membrane Extract (Trevigen, USA) and nurtured in a basal medium specifically for human ccRCC (Biozellen, USA), which encompassed advanced DMEM/F12 medium, N21-MAX and N2-MAX supplements, glutamine, nicotinamide, HEPES, A83-01, N-acetylcysteine, FGF-7, FGF-10, BMP7, EGF, Noggin, Wnt3α, R-Spondin 1, Gastrin I, Forskolin, and Y-27632. This mixture of organoid and medium was incubated in a 24-well cell incubator for a period of 3 days. To recover and propagate the organoids, Cultrex Organoids Harvesting Solution (3700–100–01, R&D Systems, USA) was utilized. Following matrix degradation, the sample underwent centrifugation at 1000 rpm for 5 min and was subsequently washed with PBS. The precipitated cells were then cultured in the aforementioned proliferation medium.

Transient transfection of cell lines and organoids with shRNA

We designed and synthesized shRNAs through RiboBio company (Guangzhou, China), with the sequences provided in Supplementary Table 3. Subsequently, these shRNAs were transfected into the cell lines using DharmaFECT4 (Thermo Scientific, USA) following the manufacturer's instructions. Twelve hours post-transfection, fresh complete medium was added to the cultures. Forty-eight hours after transfection, RT-qPCR was conducted again to evaluate the efficiencies of gene silencing.

Prior to transfection, the organoids were dissociated into single cells. The procedures for transfecting shRNAs into organoids mirrored those employed for cell lines [45]. Single cells derived from the organoids harvested to validate the transfection efficacy and to continue culturing the organoids. Briefly, they were well mixed with viral particles for 24 h at 37 °C in a culture incubator. The supernatant was removed, and the pellet was resuspended in 50 µl of ice-cold matrigel. The medium was incubated at 37 °C for 1 h to solidify the material. The organoid medium was added to start culture.

Cell counting Kit-8 and colony formation assays

Cell proliferation was assessed in vitro using the Cell Counting Kit-8 (CCK-8) and plate colony formation assays. Human ccRCC cell lines were harvested during their logarithmic growth phase and seeded into 96-well plates at a density of approximately 2,000 cells per well in 100 μl of culture media. These plates were then incubated at 37 °C in a 5% CO₂ atmosphere for 96 h. At 0, 24, 48, 72, and 96 h, 10 μl of CCK-8 reagent (APExBIO) was added to each well. Following an additional hour of incubation, the absorbance/optical density (OD) value of each well was measured at 450 nm.

For the colony formation assay, single-cell suspensions of the cell lines were prepared and seeded into a 6-well plate at a density of 1000 cells per well. The plates were incubated at 37 °C in a 5% CO₂ incubator. Once visible colonies emerged, the cultures were terminated, and the cells were fixed with 4% paraformaldehyde (Maokang Biotechnology Co., Ltd, Shanghai, China) for 30 min. Following fixation, the colonies were stained with 0.1% crystal violet for 30 min and subsequently counted.

Cell migration and invasion assays

Cell migration and invasion assays were performed following the manufacturer's instructions (Corning, 3422). Cells were harvested, centrifuged, and resuspended in serum-free DMEM (Gibco, NY, USA). For migration assays, 100 μl of the cell suspension containing 5*10⁴ cells was pipetted into the upper chamber of the Transwell insert. For invasion assays, an equivalent volume of the cell suspension was added to the Matrigel-coated upper chamber. The lower chamber was filled with 600 μl of DMEM media supplemented with 10% FBS. The Transwell plate was then placed in a cell culture incubator and incubated for 48 h. After the incubation period, the migrated or invaded cells were fixed with 4% paraformaldehyde for 20 min and subsequently stained with crystal violet for 20 min.

Wound healing assay

Cells were seeded into 6-well plates at a density of 5*10⁵ cells per well and incubated overnight in cell culture medium at 37 °C. An artificial wound was created using a P200 pipette tip, and the detached cells and debris were thoroughly washed away with PBS. Following a 24-h incubation period, images of cell migration in the vicinity of the wound were captured using a microscope.

Flow cytometric analysis of apoptosis and cell cycle

Cells in the logarithmic growth phase were harvested, digested, and suspended in cold PBS, followed by transfer to binding buffer. For apoptosis analysis, the cell suspensions were incubated with 5 µl of annexin V-FITC (KeyGen Biotech, China) and 10 µl of propidium iodide (PI, KeyGen Biotech, China) in the dark for 15 min at room temperature. Subsequently, apoptosis was detected using a FACSCalibur flow cytometer (BD Biosciences). For cell cycle analysis, the cell suspensions were fixed overnight in 70% ethanol at 4 °C. Following centrifugation of the fixed cells, they were incubated with PI in the dark for 30 min at room temperature. The cell cycle was then analyzed using a FACSCalibur flow cytometer (BD Biosciences).

Sphere-formation assay

The DMEM/F12 medium adding with B27, EGF, and bEGF were used to resuspend ccRCC cells. Two thousand of cells per well were seeded into Ultra-Low Attachment plate (Corning Incorporated, USA). We changed the medium every 2 days. After 2 weeks, we counted the spheres and acquired the pictures under a microscope (Olympus).

Immunofluorescence staining

Organoids were preserved in PBS and immobilized in 4% paraformaldehyde sourced from Maokang Biotechnology Co., Ltd, Shanghai, China, at a temperature of 4 °C for a duration of 90 min. They were subsequently blocked using 5% bovine serum albumin (BSA) for 60 min at room temperature. Following this, the organoids were incubated at room temperature for 3 h and allowed to settle due to gravitational force. Attention was directed towards the expression of Ki67 within the organoids to assess their proliferation. The samples were incubated with an anti-Ki67 antibody at 4 °C overnight. After being washed with PBS, they were then incubated with the secondary antibodies at 37 °C for 30 min. The nuclei were counterstained using DAPI for a period of 15 min. Ultimately, the samples were examined and analyzed under a fluorescence microscope.

Xenograft mouse models

All animal experimental procedures were approved by the Ethics Committee of Peking University First Hospital (Approval No. 2023-549). Five-week old BALB/c nude mice were purchased from Beijing Vital River Laboratory for tumor xenograft experiments. Caki-1 and 769-P cells (3*10⁶) with specific shRNA and negative control shRNA (shNC) were diluted and mixed with 100 µL PBS. The suspension was subcutaneously injected into mice. Major and minor axes of the tumors were measured every 3 days, and the volume was calculated by the formula: (major axe * minor axe²)/2. After 30 days, the mice were sacrificed, and the tumors were isolated and weighed.

Lung metastasis models

Tail vein injection on nude mice was performed to induce the lung metastasis models. One million Caki-1 with specific shRNA and shNC were injected into the caudal vein. After 3 weeks of tail vein injection, the mice were intraperitoneally injected with 100 mg/kg D-luciferin (Goldbio, USA). The images were captured to detect lung metastasis using the IVIS Lumina Series III imaging system (Revvity, USA). And lung tissues were harvested for gross examination.

Statistical analysis

In this study, a Student's t-test or Mann–Whitney U test was performed to compare continuous variables adhering to normal or abnormal distributions, respectively. Data analyses and graph generation were carried out using R software (version 4.3.0) and Python (version 3.12). We utilized the “limma” package to identify DEGs. The R packages “regplot” and “DynNom” were employed to construct and visualize nomograms. A P-value of less than 0.05 was considered statistically significant, unless otherwise specified.

Results

Selection of SRlncRNAs and SRmRNAs

The flowchart of this study is presented in Supplementary Fig. 1. The cohort syn3156503 comprised 78 stem cell lines and 151 non-stem cell lines. The OCLR algorithm was employed to establish the prediction model of mRNAsi. Subsequently, the RNAss of all samples in the TCGA-KIRC cohort was calculated. Initially, we performed sample clustering based on the lncRNA expression matrix to detect outliers, and no outliers were identified (Fig. 1A). In the WGCNA procedure, the soft threshold β was set to 6 (Fig. 1B). An appropriate power value was selected to construct the co-expression network, resulting in the identification of 15 modules, each represented by a distinct color (Fig. 1C). Correlations between the modules and RNAss were calculated and presented in a heatmap (Fig. 1D). Among all module-RNAss relationships, the MEgrey module exhibited the highest correlation. The eigengene from the MEgrey module was extracted, comprising 2773 lncRNAs. Within the MEgrey module, the correlation coefficient between module membership (MM) and gene significance (GS) was 0.49. Four hundred and fifty-two lncRNAs with MM > 0.2 and GS > 0.1 were selected and designated as SRlncRNAs (Supplementary Fig. 2A; Supplementary Table 4). Similar processes were applied to mRNAs. Specifically, the MEorange, MElightcyan, MEblack, and MEtan modules showed the highest correlations with RNAss (all ≥ 0.5). These modules collectively contained 1907 mRNAs (Fig. 1E–H). SRmRNAs were defined as those with MM > 0.5 and GS > 0.3, resulting in the identification of 1005 SRmRNAs (Supplementary Fig. 2B-E; Supplementary Table 5).

Fig. 1

SRlncRNAs and SRmRNAs were identified by WGCNA method. SRlncRNAs. A Sample dendrogram and trait heatmap. B Analyses of scale-free fit index and mean connectivity for various soft-thresholding powers (β). C Cluster dendrogram of modules identified using the Dynamic Tree Cut algorithm. D Heatmap of the correlation between the MEgrey module and RNAss. SRmRNAs. E Sample dendrogram and trait heatmap. F Analyses of scale-free fit index and mean connectivity for various soft-thresholding powers (β). G Cluster dendrogram of modules identified using the Dynamic Tree Cut algorithm. H Heatmap of the correlation between the MEorange/MElightcyan/MEblack/MEtan modules and RNAss

Enrichment analysis of SRmRNAs

We identified 787 DEGs from the SRmRNAs, comprising 413 upregulated and 374 downregulated mRNAs (Fig. 2A, Supplementary Table 6). Additionally, 66 upregulated and 99 downregulated SRlncRNAs were discovered (Fig. 2B, Supplementary Table 7). GO and KEGG analyses were conducted. GO analysis revealed that the DEGs were primarily enriched in BP components, such as positive regulation of neurogenesis, gland development, and positive regulation of neuron differentiation. Within the CC category, the DEGs were mainly enriched in nuclear chromatin, transcription factor complex, and nuclear transcription factor complex, among others. In the MF category, the DEGs were predominantly enriched in DNA-binding transcription factor activity and DNA-binding transcription repressor activity (Fig. 2C). The KEGG pathway analysis demonstrated that the DEGs were mostly enriched in pathways related to Epstein-Barr virus infection, transcriptional misregulation in cancer, cellular senescence, and the TGF-beta signaling pathway (Fig. 2D).

Fig. 2

Potential functions of differentially expression stemness-related genes. A-B Volcano plots showed DEGs of SRmRNAs and SRlncRNAs between tumor and normal samples. A SRmRNAs. B SRlncRNAs. C GO and D KEGG enrichment analyses for SRmRNAs

Construction of SRlncRNA signature and ceRNA network

In TCGA-KIRC cohort, one hundred and one algorithm combinations were utilized to establish an SRlncRNA signature by the LOOCV framework (Fig. 3A). The C-index of the two cohorts was employed to evaluate the predictive capabilities of the algorithms. LASSO regression combined with stepwise Cox regression exhibited optimal fitting, achieving the highest average C-index of 0.716. In this model, LASSO regression was employed to select survival-related SRlncRNAs (Fig. 3B, C). The optimal λ was achieved when the partial likelihood deviance attained its minimum value, utilizing the LOOCV framework. Ten SRlncRNAs emerged as candidates for the prognostic signature, including AC138150.2, AL109741.1, Z95115.1, LINC01711, LINC00944, AC245128.3, USP27X-AS1, AC002070.1, EMX2OS, and AL162586.1. All of these SRlncRNAs were significantly associated with the prognosis of KIRC (Supplementary Figs. 3A-J). We further conducted multivariate Cox proportional hazards analyses to determine the coefficient of each SRlncRNA (Fig. 3D). Among the SRlncRNAs included in the signature, LINC01711, LINC00944, AC245128.3, and AL162586.1 were identified as hazardous factors (HR > 1, P < 0.05), while AC002070.1 and EMX2OS were classified as protective factors (HR < 1, P < 0.05). The formula used to calculate the risk score is as follows:

$${\text{Riskscore}} = {\text{LINC}}0{1711}*0.{13673} + {\text{LINC}}00{944}*0.{21378} + {\text{AC245128}}.{3}*0.{3}0{175} - {\text{AC}}00{2}0{7}0.{1}*0.{24289} - {\text{EMX2OS}}*0.0{424}0 + {\text{AL162586}}.{1}*0.{26419}$$

Fig. 3

Establishment and validation of the prediction models. A A total of 101 types of algorithm combination and their average C-index. B-C LASSO and the matching λ value selected 10 SRlncRNAs to construct the signature. The optimal λ value was determined by minimizing partial likelihood deviance using LOOCV framework. D Multivariate cox regression for expression of SRlncRNAs

After calculating the risk score using the aforementioned formula, the median value was chosen to stratify patients into low- or high-risk groups. We summarized the risk group, M stage, T stage, grade, gender, age, survival status, and expression of SRlncRNAs in a heatmap (Supplementary Fig. 4). Survival analysis revealed that patients in the high-risk group had a poorer prognosis compared to those in the low-risk group (P < 0.001, Fig. 4A). A heatmap was constructed to depict the risk types and SRlncRNA expression levels across all samples (Fig. 4B). As the risk score increased, the survival rate of patients decreased (Fig. 4C). The AUCs for 1-, 3-, and 5-year survival rates were 0.762, 0.761, and 0.792, respectively (Fig. 4D). We retrospectively calculated the C-index of 142 published ccRCC signatures (Supplementary Table 8). It showed that SRlncRNA signature ranked as the highest among all published signatures in TCGA cohort (Fig. 4E).

Fig. 4

Development and evaluation of SRlncRNA signature related to OS in TCGA cohort. A Survival analysis of low- and high-risk groups. B Heatmap of all samples. C Scatter plot about risk score, survival time and status. D ROC curves of the signature to predict 1-, 3- and 5-year survival. E Comparisons of the C-index between SRlncRNA signature and 142 previously published signatures in TCGA cohort

Data from the GSE29609 cohort were used to validate the signature. Thirty patients with complete information were extracted, and the risk score for each patient was calculated to perform risk stratification. Patients with scores above or below the median value were classified into the high- or low-risk group, respectively. Patients in the high-risk group had significantly shorter OS compared to those in the low-risk group (P < 0.001, Fig. 5A). An expression heatmap of each patient was plotted in ascending order of risk score (Fig. 5B). As the risk score increased, the proportion of deceased patients also increased (Fig. 5C). ROC curves demonstrated that the AUC values of the signature for predicting 1-, 3-, and 5-year survival were 0.800, 0.868, and 0.851, respectively (Fig. 5D). Compared to other 99 published signatures, SRlncRNA signature manifested good predictive value in GSE29609 cohort (Fig. 5E).

Fig. 5

Validation of SRlncRNA signature related to OS inGSE29609 cohort. A Survival analysis of low- and high-risk groups. B Heatmap of all samples. C Scatter plot about risk score, survival time and status. D ROC curves of the signature to predict 1-, 3- and 5-year survival. E Comparisons of the C-index between SRlncRNA signature and 99 previously published signatures in GSE29609 cohort

To delve into the relationships between SRlncRNAs and SRmRNAs, we established a SRlncRNA-SRmRNA co-expression network (Supplementary Fig. 5A). Additionally, a Sankey diagram was generated to illustrate the risk categories associated with SRlncRNAs (Supplementary Fig. 5B). Furthermore, we constructed a ceRNA regulatory network based on SRlncRNAs and SRmRNAs (Supplementary Fig. 5C). This ceRNA network encompassed EMX2OS, four SRmRNAs, and five miRNAs. We identified five ceRNA interactions that were pertinent to stemness features (Supplementary Table 9).

Functional analysis

We further analyzed the potential biological mechanisms leading to the difference between low- and high-risk group. GSEA was performed to identify the significantly enriched KEGG and GO pathways. The results showed that alpha linolenic acid metabolism, cytokine receptor interaction, endometrial cancer, homologous recombination and intestinal immune network for IGA reproduction were mostly enriched in the high-risk group. And primary immunodeficiency, prostate cancer, proximal tubule bicarbonate reclamation, tight junction and vasopressin regulated water reabsorption were enriched in the low-risk group (Fig. 6A). In the same way, the enrichment of GO pathways was shown in Fig. 6B.

Fig. 6

The enriched KEGG A and GO (B) pathways were revealed in the high- and low-risk groups by GSEA software

Immune infiltration

Samples in the low- and high-risk groups exhibited distinct characteristics of immune infiltration. The TIMER algorithm revealed that neutrophils and myeloid dendritic cells were more abundant in the high-risk group (Fig. 7A). The CIBERSORT algorithm demonstrated that T-regulatory cells (Tregs), activated NK cells, CD8 + T cells, follicular helper T cells, M0 macrophages, memory B cells, and activated CD4 + memory T cells were significantly elevated in the high-risk group. Conversely, monocytes, resting CD4 + memory T cells, resting NK cells, activated mast cells, M2 macrophages, and naive B cells were significantly higher in the low-risk group (Fig. 7B). The ssGSEA results indicated that samples in the high-risk group had significantly higher scores for immune functions and immune cell types (Fig. 7C, D). The outcomes of other immune infiltration algorithms are presented in Supplementary Figs. 6A-E. Furthermore, the expression of immune checkpoint-related genes was significantly upregulated in the high-risk group, including PDCD1 and CTLA4 (Supplementary Fig. 6F). We employed the TIDE algorithm to predict the response to immune checkpoint inhibitors in both high- and low-risk groups (Fig. 7E). TIDE score was significantly lower in the high-risk group, indicating a better response to immunotherapy in this group. We utilized the SubMap algorithm to compare the different risk groups with a published immunotherapy cohort comprising 47 melanoma samples (Fig. 7F). The results confirmed that the high-risk group exhibited a more favorable immunotherapy response. In the CheckMate cohort, the higher the expression level of EMX2OS, the better the prognosis for ccRCC patients undergoing immunotherapy (Fig. 7G, H).

Fig. 7

Immune microenvironment and immunotherapy efficacy of samples in low- and high-risk groups. (A-B) Different immune infiltration characteristics in low- and high-risk groups. A TIMER algorithm. B CIBERSORT algorithm. C-D Different enriched scores of immune functions and immune cells in low- and high-risk groups. C Immune functions. D Immune cells. E TIDE scores of high- and low-risk groups. F Submap analysis showed the sensitivity of different risk groups to immunotherapy. In Submap analysis, a smaller p-value suggested a greater similarity between paired expression profiles. G-H Association between the expression level of EMX2OS and the prognosis of ccRCC patients undergoing immunotherapy in CheckMate cohort. G OS. H Progression free survival

Drug sensitivity analysis

The drug sensitivities of all samples were predicted, revealing that samples in the low- and high-risk groups displayed varying sensitivities to the candidate drugs. Specifically, samples in the low-risk group were more sensitive to bicalutamide, lapatinib, and thapsigargin (Fig. 8A). Conversely, samples in the high-risk group exhibited greater sensitivity to JNK.Inhibitor.VIII, RO-3306, and roscovitine. The top 100 DEGs between the low- and high-risk groups were identified. To explore potential therapies based on SRlncRNAs, the L1000FWD, CMap, and DGIdb databases were utilized. The intersection of these analyses yielded seven drugs: mebendazole, daunorubicin, loperamide, paclitaxel, testosterone, bepridil, and nadolol (Fig. 8B). Furthermore, we delved into the structures of these drugs (Fig. 8C).

Fig. 8

Potential therapeutic drugs for ccRCC based on SRlncRNA signature. A Different drug sensitivities of samples in low- and high-risk groups. The lower the IC50 value, the higher the sensitivity of the drug. B Intersection of the candidate drugs in L1000FWD, CMap and DGIdb database. C Structures of the candidate drugs in the intersection

Mutation status analysis

Mutation analysis revealed that the top five mutated genes in the TCGA-KIRC cohort were VHL, PBRM1, TTN, SETD2, and BAP1 (Supplementary Fig. 7A). Among these, missense mutations were the most prevalent variant classification, with single nucleotide polymorphisms (SNPs) constituting the most common variant type. The single nucleotide variation (SNV) class exhibiting the highest mutation frequency was "C > T" (Supplementary Fig. 7B, C). Co-mutation analyses indicated that the concurrent occurrence of mutations among the top ten mutated genes was significantly more frequent than expected by chance. Notably, SETD2-PBRM1 was the most frequently co-occurring mutation pair, whereas BAP1-PBRM1 displayed mutual exclusivity (Supplementary Fig. 7D). Both the low- and high-risk groups shared the same top four genes with the highest mutation frequencies: VHL, PBRM1, TTN, and SETD2 (Fig. 9A, B). Characteristics related to variant classification, variant type, SNV class, and variants per sample were similar between the low- and high-risk groups (P < 0.05, Fig. 9C, D). However, compared to the low-risk group, patients in the high-risk group exhibited a significantly higher frequency of mutations in BAP1 and SETD2 (Fig. 9E). Additionally, TMB was calculated for all samples, revealing that TMB was significantly higher in the high-risk group compared to the low-risk group (P < 0.01, Fig. 9F).

Fig. 9

Mutation analysis of the samples of ccRCC. A-B Waterfall plots summarized mutation frequencies and types of the top 10 most commonly mutated genes in different groups. A Low-risk group. B High-risk group. C-D Summary of the mutation status in different groups. C Low-risk group. D High-risk group. E Mutational differences between the high- and low-risk groups. F TMB of the high- and low-risk groups. **P < 0.01

Construction of the prognostic prediction model

Samples with complete demographic, clinicopathological, and prognosis records were selected to develop a nomogram. Univariate analysis identified age, grade, stage, T stage, M stage, and risk score as potential predictors of OS (Fig. 10A). Multivariate Cox regression analysis revealed that age (HR 1.032, 95% CI 1.016–1.048, P < 0.001), grade (Grade I&II-ref; Grade III-HR 1.491, 95% CI 0.994–2.235, P = 0.053; Grade IV-HR 2.152, 95% CI 1.306–3.547, P = 0.003), T stage (T1-ref; T2-HR 1.084, 95% CI 0.607–1.938, P = 0.784; T3-HR 1.762, 95% CI 1.137–2.732, P = 0.011; T4-HR 2.317, 95% CI 0.970–5.535, P = 0.059), M stage (M0-ref; M1-HR 2.592, 95% CI 1.734–3.876, P < 0.001), and risk score (HR 1.010, 95% CI 1.005–1.015, P < 0.001) were independently associated with prognosis (Fig. 10B). Based on these factors, a nomogram for predicting OS was constructed (Fig. 10C). By summing the respective points for each factor, the total score for each patient can be calculated. This nomogram integrates clinical and pathological characteristics with the novel SRlncRNA-related prognostic signature, achieving a C-statistic of 0.770 (95% CI 0.751–0.790). Moreover, we have developed a web-based calculator for clinical application (Supplementary Fig. 8, available at https://dupeng1587.shinyapps.io/KIRC_SRlncRNAs_prognosis/).

Fig. 10

Construction of the prognostic prediction model. A Forest plot of univariate cox regression. B Forest plot of multivariate cox regression. C Nomogram to predict prognosis of ccRCC based on clinic-pathological factors and SRlncRNAs

Single cell landscape and spatial sequencing of SRlncRNAs

We analyzed the scRNA-seq landscape of ccRCC to clarify the expression patterns of SRlncRNAs and their potential roles in tumor progression. The overall cellular composition remained consistent among the specimens after filtering steps (Supplementary Fig. 9A, B). Following data normalization, we identified the top 2000 variable features for downstream analysis (Supplementary Fig. 9C). The dimensionality of the data was determined using the JackStraw and ElbowPlot methods (Supplementary Fig. 9D, E). Subsequently, thirty-three clusters were identified and visualized using the Uniform Manifold Approximation and Projection (UMAP) method (Fig. 11A; Supplementary Fig. 9F). Based on previous reports and marker gene expression, we classified the cells into eight distinct types: malignant cells, monocytes/macrophages, epithelial cells, plasma cells, endothelial cells, CD8 T cells, erythroblasts, and pericytes (Fig. 11B; Supplementary Fig. 9G). The expression of EMX2OS and LINC00944 was recorded in GSE159115 cohort. Depending on the expression levels of these SRlncRNAs, we further stratified malignant cells into high and low expression groups for both EMX2OS and LINC00944 (Fig. 11C, D). When comparing the expression of EMX2OS and LINC00944 between epithelial and malignant cells, we observed that the expression of EMX2OS was significantly lower in malignant cells compared to epithelial cells, whereas the expression of LINC00944 was significantly higher (P < 0.05, Fig. 11E–G). Pseudotime series analysis illuminated the developmental trajectory of cells from epithelial to malignant states (Fig. 11H, I). Both epithelial and malignant cells exhibited three distinct states during this transition (Supplementary Fig. 10A). Notably, during this developmental process, the expressions of EMX2OS and LINC00944 underwent significant changes (Fig. 11J, K).It suggested their potential involvement in mediating the transition from normal to malignant phenotypes. Upon further subdivision of malignant cells into EMX2OS high- and low-expressed groups, as well as LINC00944 high- and low-expressed groups, the cells maintained their presence in three distinct states, with only subtle alterations noted in their developmental trajectory (Fig. 11L–O; Supplementary Fig. 10B, C).

Fig. 11

Single cell landscape of ccRCC in GSE159115 cohort. A Thirty-three clusters were identified by UMAP method. B Cells were classified into CD8 T cells, plasma cells, monocyte/macrophage, endothelial cells, epithelial cells, malignant cells, erythroblasts and pericytes. C-D Malignant cells were divided into 2 clusters based on the expression of EMX2OS and LINC00944. C EMX2OS. D LINC00944. E–G Expression of SRlncRNAs in different cell types. E EMX2OS. F LINC00944. H–K Pseudotime trajectory analysis of epithelial and malignant cells in GSE159115 cohort. H Pseudotime was colored in gradients from epithelial cells to malignant cells. I Pseudotime of epithelial and malignant cells. J-K Correlations between the pseudotime trajectory and the expression of EMX2OS or LINC00944. L-M Pseudotime trajectory analysis of epithelial, EMX2OS-low-expressed, and EMX2OS-high-expressed malignant cells. N–O Pseudotime trajectory analysis of epithelial, LINC00944-low-expressed, and LINC00944-high-expressed malignant cells

We computed the CNV scores for each cell type, revealing distinct characteristics in malignant cells (Fig. 12A, B). Further analysis using the Leiden algorithm classified the cells into 14 distinct clusters (Supplementary Fig. 10D, E). Notably, the CNV score of malignant cells was significantly elevated compared to other cell types (Fig. 12C). Within the malignant cell population, cells with high LINC00944 expression exhibited a higher CNV score than those with low LINC00944 expression (Fig. 12D). We also delved into the associations between SRlncRNA expression and CytoTRACE scores (Fig. 13A–D). Among all cell types, malignant cells displayed the highest CytoTRACE scores (Fig. 13E), and the expressions of both EMX2OS and LINC00944 were significantly correlated with CytoTRACE scores.

Fig. 12

CNV landscape of ccRCC in GSE159115 cohort. A CNV heatmap of all types of cells. B CNV scores of all cell types in UMAP figure. C Different CNV scores in the 8 types of cells. D Different CNV scores in epithelial, EMX2OS-low-expressed, EMX2OS-high-expressed, LINC00944-low-expressed, and LINC00944-high-expressed malignant cells

Fig. 13

CytoTRACE landscape of ccRCC in GSE159115 cohort. A CytoTRACE identified the stemness in different cell types. B-D Phenotype plots showed the clusters of different cell types. E CytoTRACE scores of different cell types

To gain deeper insights into the spatial expression patterns of EMX2OS and LINC00944 in ccRCC, we conducted a comprehensive analysis of spatial transcriptome sequencing data obtained from a representative tumor-normal interface (Fig. 14A, D). Cell type annotation was shown in Fig. 14C, F. The expression of EMX2OS and LINC00944 is unevenly distributed across the boundary between tumo and normal tissues in ccRCC (Fig. 14B, E).

Fig. 14

Spatial transcriptome sequencing data of 2 slices from a ccRCC patient. A Spatial mapping of a representative tumor-normal interface in slice 1. B Expression of EMX2OS (ENSG00000229847) and LINC00944 (ENSG00000256128) in slice 1. C Annotation of the spots in slice 1. D Spatial mapping of a representative tumor-normal interface in slice 2. E Expression of EMX2OS (ENSG00000229847) and LINC00944 (ENSG00000256128) in slice 2. F Annotation of the spots in slice 2

Expression of EMX2OS and LINC00944 in ccRCC and their impacts on prognosis

We performed RT-qPCR to assess the expression levels of EMX2OS and LINC00944 in renal cancer cell lines and tissue specimens. Our results indicated that LINC00944 was significantly upregulated in Caki-1, 786-O, and 769-P renal cancer cell lines compared to the normal HK-2 cell line (all P < 0.05; Supplementary Fig. 11A). Conversely, EMX2OS was downregulated in renal cancer cell lines (P < 0.05). A comparable trend was observed in 102 tissue specimens (Supplementary Fig. 11B), where the expressions of EMX2OS and LINC00944 were significantly distinct between normal and renal cancer samples (all P < 0.05). Survival analysis further confirmed that all six SRlncRNAs were significantly correlated with the prognosis of ccRCC (P < 0.05; Supplementary Fig. 11C, D). Through the analysis of TCGA Pan-Cancer cohort, we further elucidated the expression patterns of EMX2OS and LINC00944 across various cancer types, as well as their associations with prognosis (Supplementary Fig. 12A-D).

Downregulation of EMX2OS and LINC00944 in ccRCC influenced proliferation, migration and invasion in vitro

We utilized targeting shRNAs to specifically knockdown the expressions of EMX2OS and LINC00944 in ccRCC cell lines, achieving effective knockdown efficiencies for both targets (Fig. 15A). To evaluate the proliferative abilities of these cells, CCK-8 and colony formation assays were performed. Our results demonstrated that knocking down EMX2OS significantly promoted the proliferation of ccRCC cells, whereas knocking down LINC00944 significantly inhibited their proliferation (Fig. 15B, C). Additionally, we conducted Transwell assays to investigate the effects of EMX2OS and LINC00944 on cell migration and invasion. Knockdown of EMX2OS enhanced cell migration and invasion, whereas knocking down LINC00944 suppressed these processes (Fig. 15D, E; P < 0.01 for all comparisons). Furthermore, wound healing assays were utilized to assess the migratory capacities of ccRCC cells, revealing an increase in migratory capacity following EMX2OS knockdown and a decrease following LINC00944 knockdown (Fig. 16A; P < 0.001 for all comparisons).

Fig. 15

Downregulation of EMX2OS promoted ccRCC cell proliferation, migration and invasion, while downregulation of LINC00944 suppressed them. A Assessment of relative expression levels of EMX2OS and LINC00944 in Caki-1 and 769-P cell lines, following transfection with shNC and shRNA constructs, was conducted by RT-qPCR. B Cell proliferation ability was evaluated by CCK-8 assay in Caki-1 and 769-P cell lines. C The ability of cell proliferation was evaluated by colony-forming unit assay in Caki-1 and 769-P cell lines. D–E The abilities of cell migration and invasion were evaluated by transwell assays in Caki-1 and 769-P cell lines. D Migration ability. E Invasion ability. Error bars showed the SD from three independent experiments. **P < 0.01, ***P < 0.001

Fig. 16

Downregulation of EMX2OS promoted ccRCC cell migration, suppressed apoptosis, and increased S phase entry, while downregulation of LINC00944 appeared the opposite effects. A The ability of cell migration was promoted and suppressed after knockdown of EMX2OS and LINC00944, respectively. B Apoptosis rates of Caki-1 and 769-P cells were evaluated by flow cytometry after transfection. C Cell cycle analysis of Caki-1 and 769-P cells were conducted by flow cytometry after transfection. *P < 0.05. **P < 0.01. ***P < 0.001

Downregulation of EMX2OS and LINC00944 influenced apoptosis and cell cycle of ccRCC cells

We analyzed apoptosis and cell cycle progression using flow cytometry. Knockdown of EMX2OS was found to inhibit apoptosis in ccRCC cell lines, whereas knockdown of LINC00944 promoted apoptosis in ccRCC cells (Fig. 16B; P < 0.05 for all comparisons). Specifically, upon EMX2OS knockdown, ccRCC cells in the G0/G1 phase decreased, with a corresponding increase in cells in the S phase. Conversely, knockdown of LINC00944 led to an increase in ccRCC cells in the G0/G1 phase and a decrease in cells in the S phase (Fig. 16C).

Downregulation of EMX2OS and LINC00944 enhanced the stemness of ccRCC

We explored the changes in the stemness and its related genes after the silencing of EMX2OS and LINC00944. Downregulation of EMX2OS upregulated the mRNA and protein levels of CD133, EPCAM, and SOX2, and promoted sphere formation (Fig. 17A, B). Conversely, downregulation of LINC00944 downregulated the expression of stemness-related genes, and inhibited sphere formation (Fig. 17A, B). Organoid model derived from a ccRCC patient was successfully established. We transfected shRNAs of EMX2OS and LINC00944 into ccRCC organoids (Fig. 17C). Knockdown of EMX2OS promoted the growth of ccRCC organoids, whereas silencing LINC00944 inhibited the growth of ccRCC organoids (Fig. 17D). Immunofluorescence demonstrated that upon silencing SRlncRNA, the expression of KI67 in the organoids was altered, which was consistent with the changes observed in the growth of the organoids (Fig. 17E).

Fig. 17

Downregulation of EMX2OS enhanced the stemness of ccRCC, whereas Downregulation of LINC00944 reduced the stemness of ccRCC. A Expression of stemness markers (CD133, EPCAM, and SOX2) were upregulated and downregulated after knockdown of EMX2OS and LINC00944, respectively. B Downregulation of EMX2OS and LINC00944 affected spheroid formation ability of ccRCC cells. C The transfection efficiency of shNC and shRNA in ccRCC organoids. D The growth status of ccRCC organoids before and after transfection with shNC and shRNA. E The expression levels of KI67 in ccRCC organoids after transfection with shNC and shRNA. F The subcutaneous tumorigenesis of CAKI1 and 769P cells in nude mice was observed following transfection with shNC and shRNA. Additionally, the variations in tumor volume and weight over a 30-day period post-treatment were documented. G The lung metastatic tumor formation of CAKI1 cells following transfection with shNC and shRNA. *P < 0.05. **P < 0.01. ***P < 0.001

To further assessed the effects of EMX2OS and LINC00944 on ccRCC growth and metastasis in vivo, we generated xenograft mouse and lung metastasis models. Xenograft tumor model indicated that both tumor size and weight were significantly increased in sh-EMX2OS group, while decreasing in sh-LINC00944 group (Fig. 17F). We observed the similar trends in lung metastasis models. There were more lung metastasis lesions in sh-EMX2OS group, and less lesions in sh-LINC00944 group (Fig. 17G).

Discussion

A previous study utilizing the SEER database indicated that the prognosis of ccRCC has improved in recent years [46]. However, further reducing the mortality rate of ccRCC continues to pose a significant challenge. The acquisition of stemness potential is a crucial process when a non-stem cell acquires malignant traits [47]. It is widely recognized that stemness is intimately linked to cancer pathogenesis, progression, and therapy resistance [48]. The stemness phenotype exhibits considerable heterogeneity within and across different tumor types [49]. In this study, we have established an SRlncRNA-related signature for predicting the stemness and prognosis of ccRCC. Our findings demonstrate that this signature performs exceptionally well. Prior studies have shown that SRlncRNA-related signatures can effectively predict the prognosis of various tumors, including gastric cancer, colorectal cancer, breast cancer, and lung adenocarcinoma [14, 15, 50, 51]. To the best of our knowledge, our study presents the first SRlncRNA-related signature specifically designed to predict the prognosis of ccRCC.

In previous studies, multiple signatures have been developed to predict the prognosis of ccRCC. Qu et al. conducted a multicenter cohort study involving 2401 patients and established a four-lncRNA signature, termed RCClnc4, comprising ENSG00000255774, ENSG00000248323, ENSG00000260911, and ENSG00000231666 [52]. This signature demonstrated potential in predicting the survival of patients with stage I-III ccRCC, achieving a C-index of 0.720. Additionally, Gui et al. constructed an 8-lncRNA signature associated with autophagy [53], which revealed the ability to facilitate precise risk stratification and prognostic prediction, yielding an AUC of 0.764. Our stemness-related signature for ccRCC was established based on LASSO combined with stepwise Cox regression. LASSO encouraged simpler models by shrinking some coefficients precisely to zero, thereby efficiently performing feature selection and diminishing the model's complexity. This signature achieved an average C-index of 0.716, with AUC predictive values for 1-, 3-, and 5-year survival standing at 0.800, 0.868, and 0.851, respectively. These indices underscore the satisfactory predictive performance of our signature.

The loss of differentiation and the acquisition of stemness are culpable for tumor progression. Cancer stem cells (CSCs), a crucial subpopulation of cancer cells, possess the capacity to divide and differentiate into mature cells [54]. CSCs play a pivotal role in treatment resistance, tumor recurrence, and metastasis. Their presence in ccRCC has been confirmed, offering a promising new target for anticancer drugs [55]. Furthermore, the CSC theory may provide insight into explaining the tumor heterogeneity observed in ccRCC [49]. Malta et al. identified the biological mechanisms underlying the dedifferentiated oncogenic state and developed two independent stemness indices, mDNAsi and mRNAsi, utilizing the OCLR machine learning algorithm [18]. These stemness indices positively correlate with active biological processes associated with CSCs and the dedifferentiated oncogenic state [18]. Wang et al. classified ccRCC into three subtypes based on the expression of stem cell-related genes: stem cell activated, stem cell dormant, and stem cell excluded [56]. Notably, the stem cell activated subtype, characterized by the highest concentration of CSCs, exhibits the highest level of immune cell infiltration, a distinct mutation landscape, and the poorest prognosis. Additionally, this subtype demonstrates the highest sensitivity to immunotherapy, particularly CTLA-4 blockade, and to the vascular endothelial growth factor receptor (VEGFR) inhibitor sunitini [56]. LncRNAs can influence the expression of downstream genes through various mechanisms, including functioning as a ceRNA, activating or inhibiting transcription factors, and other pathways. We will further explore the impact of the expression of SRlncRNAs in ccRCC on immunotherapy efficacy and drug sensitivity [57].

GO analysis revealed that stemness-related genes were predominantly enriched in the positive regulation of neurogenesis, nuclear chromatin organization, and DNA-binding transcription activator activity, RNA polymerase II-specific. As is widely recognized, mTOR plays a crucial role in cell growth and proliferation, and tumor growth and angiogenesis can be suppressed by mTOR inhibitors. Neurogenesis, nuclear chromatin organization, and DNA-binding transcription are all vital processes in the progression of ccRCC. The results of KEGG enrichment pathways encompassed transcriptional misregulation, cellular senescence, and the TGF-beta signaling pathway. These findings align with previous studies indicating that stemness levels are intimately linked to senescence-related transcription. Boguslawska et al. suggested that TGF-beta1 can coordinate the expression of a microRNA network that regulates cellular adhesion in ccRCC [58]. Nishida et al. demonstrated that decreased TGF-beta glycan expression enhances the metastatic potential of ccRCC [59]. Our study further explores the relationship between stemness and the TGF-beta signaling pathway in patients with ccRCC.

This SRlncRNA-related signature comprises six lncRNAs, which exhibit differential expression patterns between normal and ccRCC samples. Our research has further validated the expression of EMX2OS and LINC00944 (also known as DMDRMR) in ccRCC cell lines and specimens, along with their impacts on the prognosis of ccRCC. Prior studies have reported associations between LINC00944 and EMX2OS with the progression of ccRCC. Chen et al. demonstrated that LINC00944 can regulate TYMP expression and suppress Akt phosphorylation in ccRCC cells [60]. Meanwhile, EMX2OS overexpression has been shown to inhibit the growth of Wilms' tumor by targeting miR-654-3p [61]. Furthermore, EMX2OS restricts the stemness, epithelial-mesenchymal transition, and metastasis of Wilms' tumor cells through direct interaction with the IGF2BP1 protein. LINC00944, under the alias of DMDRMR, signifies a DNA methylation-deregulated and RNA m6A reader-cooperating lncRNA. Gu et al. reported that DMDRMR promotes the growth and metastasis of ccRCC by binding to IGF2BP3, an m6A reader that stabilizes its target genes, including CDK4, COL6A1, LAMA5, and FN1. In an m6A-dependent manner, DMDRMR specifically enhances the activity of IGF2BP3 on its target genes [62]. Through the ceRNA mechanism, DMDRMR can promote angiogenesis in ccRCC by acting as a sponge for miR-378a-5p, thereby increasing the expression of EZH2 and SMURF1, both of which repress DAB2IP expression. This DMDRMR/miR-378a-5p/DAB2IP axis activates the VEGFA/VEGFR2 pathway, ultimately promoting angiogenesis in ccRCC [63]. Located in the enhancer region of the EMX2 gene, expression of EMX2OS showed a strong correlation with EMX2 in gastric cancer [64]. Our study confirmed that EMX2OS and LINC00944 exerted influence on the proliferation, migration, invasion, apoptosis, cell cycle, and stemness of ccRCC cells. These findings hinted at the potential of SRlncRNAs to impact the progression of ccRCC from multiple facets. Despite making considerable efforts in this regard, there remained a substantial gap between computational predictions and the realm of clinical trials [65].

Our study has shed light on the potential roles of SRlncRNAs in the progression of ccRCC, albeit with notable limitations that warrant attention. Firstly, the signature was constructed solely based on bioinformatic methodologies, inherently introducing the possibility of biases. Small sample size in the validation cohort could lead to an increased risk of both false-positive and false-negative results. Consequently, external validation through real-world analyses is imperative to corroborate our findings. Secondly, the nomenclature of lncRNAs may vary across different databases, posing challenges in the standardization and comparison of studies. Thirdly, most of the analytical processes in this study relied on retrospective cohorts, which might introduce potential biases. The findings of this study necessitated further validation through prospective research. Lastly, to firmly establish the mechanistic underpinnings of SRlncRNAs in driving ccRCC progression, further mechanism explorations are indispensable.

In conclusion, we have identified the expression patterns and potential functional roles of SRlncRNAs in ccRCC. We delved into the potential contributions of SRlncRNAs to the prognosis, drug sensitivity, and immune infiltration of ccRCC. Furthermore, we established a SRlncRNA-related signature to predict risk stratification, which could potentially serve as a valuable tool in informing clinical decision-making for ccRCC patients.