Identification of Immune-Related Seven-Long Non-Coding RNA Signature for Overall Survival and Validation of the Effect of LINC01270 in Malignant Phenotypes of Clear Cell Renal Carcinoma

Introduction

Kidney cancer is a serious threat to human health. Approximately 73,750 new kidney cancer cases were diagnosed in the United States in 2020.¹ The 5-year survival rates of patients with localized kidney cancer (ie, stage I/II in clinical) have increased to 90%.² However, the overall prognosis of patients who present with high-stage disease (ie, stage III/IV) is still poor.³ Clear cell renal cell carcinoma (ccRCC) is the most common subtype of kidney cancer.⁴ Patients with localized RCC usually have a better prognosis after surgical removal of part or all of the kidney.⁵ However, late-stage ccRCC (stage III/IV) generally cannot be cured because they usually have metastatic lesions.⁵ Patients with late-stage ccRCC may be treated with surgery (cytoreductive nephrectomy and/or metastasectomy), sometimes combined with anti-angiogenic therapy and immunotherapy or with exclusive systemic treatment.^4,6 To date, several prognostic scores have been proposed for the prognosis of ccRCC patients, such as the stage, size, grade, and necrosis (SSIGN) score⁷ and The University of California Los Angeles Integrated Staging System (UISS) score.⁸ However, they were less useful for individualized risk prediction.³ For the prognosis of advanced RCC, an IMDC prognostic score was also proposed.⁹ Moreover, these prognostic systems depend on evaluating clinical parameters without consideration of molecular markers.¹⁰ Therefore, the identification of molecular-level prognostic biomarkers is warranted for the clinical management of late-stage ccRCC.

In recent years, medical treatment for ccRCC has transitioned to a novel immunotherapy era.³ Identification of immune-related genes contributes to understanding the molecular basis of the effects of immunotherapy in ccRCC. Protein-coding genes that can be immunotherapy targets are under investigation, such as immune checkpoint molecules.¹¹ Long non-coding RNA (lncRNA) is an important type of RNA. LncRNAs were closely associated with cell cycle, proliferation, metastasis, and drug resistance of cancer cells.¹² Functionally, they may act as signal, decoy, scaffold, guide, enhancer RNAs, and short peptides. For example, the classic lncRNAs HOTAIR, upregulated in diverse cancer types, was found to promote invasion and metastasis by regulating the H3K27 methylation of tumor suppressor genes.¹³ Recently, lncRNAs have been linked to PD-L1/PD-1 checkpoint pathways, implying the clinical potential of lncRNAs.¹⁴ Thus, lncRNAs may serve as novel biomarkers in the future.¹² However, the prognostic value of immune-related lncRNAs has not been fully investigated. Thus, we systematically assessed the prognostic potential of immune-related lncRNAs in ccRCC by leveraging the knowledge of MSigDB C7 immunologic signatures.¹⁵

Materials and Methods

Transcriptome Data

The raw CEL files (GSE53757¹⁶ and GSE66270)¹⁷ from the Affymetrix Human Genome U133 Plus 2.0 Array platform were retrieved from the Gene Expression Omnibus (GEO) database. For the GSE53757 dataset, a total of 144 ccRCC samples (72 tumor samples and 72 normal samples) were downloaded. For the GSE66270 dataset, a total of 28 ccRCC samples (14 tumor samples and 14 normal samples) were downloaded. Moreover, the level 3 gene expression data (HT-Seq read counts) of TCGA-KIRC (n = 603) were retrieved from the Genomic Data Commons (GDC) data portal (https://portal.gdc.cancer.gov/repository). The RNA-Seq data (read count) of late-stage ccRCC (stage III/IV, n = 186) of the ICGC project were downloaded from the Xena data hub (https://xena.ucsc.edu/public). Follow-up clinical data were also retrieved for analysis. The.gtf annotation file from GENCODE (v22) database was used to obtain the lncRNA records. The genes with the identifier “lincRNA” in the annotation file were extracted and analyzed in this study.

Differential Expression Analysis

For microarray expression profiles, the background correction was performed for the signal intensities. For comparison between sample groups, quantile normalization was performed for samples in the same dataset. The expression matrix of lncRNA genes was retained for differential expression analysis. The limma package¹⁸ was used for differential expression analysis of microarray data. For the read count of RNA-Seq data, the DESeq2 package¹⁹ was used. Genes with Benjamini–Hochberg (BH) adjusted P-value <0.05 were statistically significant.

Identification of Immune-Related lncRNAs and the Function Analysis

Differentially expressed protein-coding genes in significant enriched Molecular Signatures Database (MSigDB) C7 immunologic signatures¹⁵ were considered immune-related genes. The lncRNAs with Pearson correlation |r| ≥0.4 and FDR-value <0.05 with BH adjustment were determined as immune-related lncRNAs. A functional analysis of the lncRNA signature was performed based on the lncRNA-associated protein-coding genes. Protein-coding genes correlated with lncRNAs were determined using Pearson correlation analysis. The mRNA-lncRNA pairs with FDR <0.05 were considered. The mRNAs in the “pathway in cancer” category from the PathCards database were included in this correlation analysis. Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathways were included in the enrichment analysis.

Statistical Analysis

Univariate and Multivariate Cox regressions were used to identify prognostic lncRNAs in ccRCC. LncRNAs with BH-adjusted P < 0.05 in univariate analysis were subjected to multivariate analysis to obtain the lncRNA signature. Similar to Wang et al,²⁰ The prognostic index (Pindex) was calculated using the lncRNA expression levels and Cox regression coefficients of the seven lncRNAs:

where n is the number of lncRNAs, is the expression of lncRNA i, and is the estimated regression coefficient of lncRNA i in the multivariate Cox model. The ccRCC patients were divided into two groups based on Pindex value (median as cut-off). Clinical parameters were added into the multivariate Cox hazard model to assess whether the Pindex was an independent factor associated with the prognosis of ccRCC, such as TNM stage, grade, age, and gender.

Clinical Specimens

The ccRCC tissues (n = 20) and the paired normal tissues from 20 patients were obtained in the First Affiliated Hospital with Chongqing Medical University (Chongqing, China). The use of tissues has been abided by the ethics committee.

Cell Lines and Cell Culture

The ccRCC cell lines ACHN, Caki-1, 769-P, Caki-2, and 786-O (RPMI 1640) and human renal tubular epithelial cells HK-2 were purchased from the Cell Bank of the Chinese Academy of Sciences. The HK-2 cells were cultured in Keratinocyte Serum Free Medium (K-SFM), the ACHN cells were cultured in DMEM, the Caki-1 and Caki-2 cells were cultured in McCoy’s 5A, and 786-P and 786-O cells were cultured in RPMI 1640 with 10% fetal bovine serum FBS (PAN-Biotech, Adenbach, Germany) and 100 U/mL penicillin-streptomycin. All cells were cultured at 37°C with 5% CO2.

Transfection

Si-LINC01270 and NC were purchased from Gene Pharma (Shanghai, China). Transfections were performed using GP-transfect-mate (Gene Pharma, China) according to the manufacturer’s instructions. During the transfection process, many negatively charged proteins in the serum may interfere with the adsorption of cationic liposomes on nucleic acids, affecting the transfection efficiency. We use a serum-free medium for transfection. After 6h of transfection, cells were cultured in complete media for 48h.

RNA Isolation and qRT-PCR

According to the manufacturer’s instructions, total RNA was extracted by Trizol (Invitrogen Corporation, Carlsbad, CA, USA), and 0.5 ug RNA was reverse transcribed using a Reverse Transcription Kit (Takara, China). PCR used SYBR Green Master Mix (Takara, China). GAPDH was used as endogenous control, and all primers used for qRT-PCR are listed as follows: LINC01270 Forward: 5’-CTCACGAAAGCGCAGGAATG-3’; Reverse: 5’-GCTCCAAAAGCAGACAAGCC-3’. FIRRE Forward: 5’-CACTGGAGCAATCCTGGCTG −3’; Reverse: 5’-TGCCTAACTACCTCTATTGGCT-3’. RP11-37B2.1 Forward: 5’-GCGATTTGGGGAGTGTTGTG-3’; Reverse: 5’-GTTGGGGCAGTAGGCAACTA-3’. RP11-253I19.3 Forward: 5’-TGTCGGATTGGTTAGCGACC-3’; Reverse: 5’-CTTAAGTGTGGAGCCCTCGG-3’. RP11-438L19.1 Forward: 5’-TGGTTCTGCTCCTGGTAACG-3’; Reverse: 5’- AGGATAACAGGTCTGCCTGC-3’. RP11-504P24.9 Forward: 5’-TGTGTCCCCAAGGAAGGATG-3’; Reverse: 5’-TCCACACAGTGTAGTCAAGCC-3’. CTB-41I6.1 Forward: 5’-GTGTCCCCAAGGAAGGATGAG-3’; Reverse: 5’- GTAGTCAAGCCGACTCTCCA −3’. GAPDH Forward: 5’-CCATGGGGAAGGTGAAGGTC-3’; Reverse: 5’- AGTGATGGCATGGACTGTGG −3’.

Cell Proliferation and Colony Formation Assays

Cell proliferation was measured by Cell Counting Kit-8 reagent (CCK-8) (Bimake, China). Briefly, the cells were cultured in 96-wells plates, and 10 ul CCK-8 solution was added to each well every day. Then, the cells were cultured at 37°C for 1h and were detected at 450 nm on a Microplate Reader (BIO-TEK, USA). For colony formation assays, 1.5×103 cells were plated in 6-wells plates and cultured for 2 weeks. Then, the cell colonies were fixed, stained, and counted.

Cell Migration and Invasion Assays

The cells were suspended in 200 ul serum-free medium and inoculated into the upper chamber (BD Company, USA) uncoated (for migration analysis) or coated with Matrigel. Then, the lower chambers were fixed with 600 ul 20% FBS serum. Migrating and invading cells were fixed and stained after 36h. Chambers were randomly selected from three fields of view to count, and the average count was considered the number of permeable cells.

Results

Identification of ccRCC-Related lncRNAs and Immune-Associated Coding Genes

In the present study, we analyzed the differentially expressed immune-related lncRNAs with prognostic potential by leveraging the knowledge of immune-related gene signatures from the MSigDB database (Figure 1). First, we identified 2342 protein-coding genes (FDR < 0.05 with BH adjustment, Figure 2A) and 893 lncRNAs (FDR < 0.05 with BH adjustment, Figure 2B) differentially expressed between tumor samples and normal controls based on one TCGA dataset (TCGA-KIRC, Table 1) and two GEO datasets (GSE53757 and GSE66270, Table 1). To determine immune-related protein-coding genes, gene set enrichment analysis based on MSigDB C7 immunologic signatures showed that 1105 differentially expressed protein-coding genes were significantly enriched in 51 signatures (FDR < 0.01 with BH adjustment, Figure 2C), including CD8 T cell, CD4 T cell, and NK cell-associated signatures. Here, 1105 genes were considered immune-related genes associated with ccRCC.

Table 1 Patient Clinical Information of ccRCC Datasets

Figure 1 Flowchart of the strategy to identify immune-related prognostic lincRNAs in ccRCC. Abbreviation: MSigDB, Molecular Signature Database.

Figure 2 Identification of ccRCC-related lncRNAs and immune-associated coding genes. (A) Identification of differentially expressed (DE) protein-coding genes between ccRCC tumor sample and normal controls based on TCGA and GEO databases. (B) Identification of DE lncRNA between tumor and normal. (C) Gene set enrichment analysis of DE protein-coding genes based on MSigDB C7 immunologic signatures. Genes and signatures with Benjamini-Hochberg (BH) adjusted P-value < 0.05 were significant.

Identification of Immune-Related lncRNAs

We performed a Pearson correlation analysis to identify differentially expressed lncRNAs with potentially functional relevance to immune-related genes. The hypothesis is that if a lncRNA is functionally associated with a protein-coding gene, it should be co-expressed. The lncRNA-gene pairs with Pearson correlation coefficients |r| ≥0.4 and BH-adjusted P <0.05 were significant co-expression. As a result, a total of 290 differentially expressed lncRNAs were co-expressed with immune-related genes (Figure 3). Hierarchical clustering of the Pearson correlation matrix showed that 290 lncRNA was associated with immune-related genes involved in activities of CD4/CD8 T cells, B cells, and natural killer (NK) cells (Figure 3).

Figure 3 Identification of immune-related lncRNAs in ccRCC. Pearson correlation coefficient matrix of 290 differentially expressed lncRNAs (y-axis) and 1105 immune-related protein-coding genes (x-axis). For better visualization, the cells in the heatmap (correlation of specific gene-lncRNA pair) with BH-adjusted P < 0.05 and Pearson correlation coefficients |r| < 0.4 were considered no significance and were replaced with zero. The Euclidean distance metric with Ward linkage was used for hierarchical clustering of the correlation matrix. The functional classification of the four clusters of immune-related genes was performed based on MSigDB C7 signatures (BH-adjusted P < 0.05 was used as a cut-off value).

7-lncRNA Signature Could Be an Independent and Unfavorable Prognostic Biomarker for Late-Stage ccRCC

Univariate Cox proportional hazard regression based on 225 immune-related lncRNAs was performed to identify lncRNAs associated with overall survival of patients with late-stage ccRCC. Here, to make full use of available data and make the results more robust, the TCGA cohort and the ICGC cohort were used for training and validation, respectively. Results showed that 82 lncRNAs were associated with overall survival of patients with late-stage ccRCC (P < 0.05, Log rank test). Further multivariate analysis using the 82 lncRNAs showed a 7-lncRNA signature (LINC01270, FIRRE, RP11-37B2.1, RP11-253I19.3, RP11-438L19.1, RP11-504P24.9 and CTB-41I6.1, log-rank P = 6.04E-06, Figure 4A). The scoring method (Pindex) was calculated based on the 7-lncRNA signature. To investigate whether the 7-lncRNA model was affected by clinical factors, multivariate regression was performed using the Pindex and clinical parameters as candidate variables. The Kaplan–Meier curves were used to evaluate the survival rates of the 7-lncRNA signature. In this multivariate model, the Pindex values of patients were divided into two groups (median as cut-off value). Results showed that our 7-lncRNA signature could be an independent factor for the overall survival prognosis of late-stage ccRCC in both training (P = 2.14E-04, Figure 4B, Table 2) and validation datasets (P = 0.003, Figure 4C, Table 2). The Pindex-high group patients have significantly poor survival compared with low-Pindex patients (Figure 4B and C).

Table 2 The Univariate and Multivariate Cox Proportional Hazard Model for Stage III/IV ccRCC Patients

Figure 4 The prognostic potential of our 7-lncRNA signature. (A) The expression of the seven lncRNAs in the TCGA cohort. (B and C) Overall survival of patients with late-stage ccRCC in training (TCGA, n = 162) and validation datasets (ICGC, n = 158). The prognostic index (Pindex) was calculated based on the expression and Cox regression coefficients of the seven lncRNAs. The median value of Pindex was used as a cut-off to divide patients into high- and low-Pindex groups.

The Functions of the 7 lncRNAs and the Correlations with Immunotherapy Targets

Immune checkpoint blockade therapies were effective for diverse cancer types, such as ccRCC. The major targets of immune checkpoint blockade therapies are PD-1 (PDCD1 or CD279), PD-L1 (CD274), and CTLA4. We performed correlation analysis between immunotherapeutic targets and the seven lncRNAs. It showed that they have universal associations with immunotherapeutic targets (Figure 5). For example, both RP11−504P24.9 and CTB-41I6.1 were associated with the expression of PD-1 and CTLA4 (P values <0.05, r >0.2). RP11−37B2.1 was associated with the expression of PD-L1 and CTLA4 (P values <0.05, r >0.2, Figure 5). We further explored the potential functions and pathways to investigate how the 7-lncRNA signature is involved in tumor immunity. We showed that 260 genes of “pathway in cancer” were related to the seven lncRNAs (FDR <0.05 and Pearson |r| >0.4). They were involved in many cancer-related signaling pathways, such as MAPK, cell adhesion-related pathways, ERBB, Wnt, apoptosis, T/B cell, and chemokine signaling (Figure 6A). Based on GO enrichment analysis, the cellular components are protein receptors/complexes in the cell surface and extracellular matrix, such as immune-related CD40 receptor complex (Figure 6B). Moreover, based on the enrichment results of molecular functions, the protein receptors/complexes were growth factors, transcriptional factors, kinase, and cytokine/cytokine receptors (Figure 6C). Biological process enrichment results indicated that these protein receptors/complexes were involved in proliferation and apoptosis (Figure 6D). In summary, the seven lncRNAs may play an important role in tumor immune cells and tumor progression through regulating the binding of protein receptors/complexes.

Figure 5 The correlations between PD-1/PD-L1/CTLA4 and the 7 lncRNAs in ccRCC. The three genes PD-1/PD-L1/CTLA4 are the most effective drug targets of immune checkpoint blockade therapy. The Pearson correlation was used in this analysis.

Figure 6 The functions and pathways of the 7 lncRNAs. (A) The KEGG pathways analysis enriched based on the protein-coding genes associated with the 7 lncRNAs. (B–D) The GO cellular components, molecular function, and biological process enrichment results, respectively.

Experimental Validation of the Expression of the Seven lncRNAs in ccRCC Tissues and Cell Lines

qRT-PCR was performed to examine the expression level of the seven lncRNAs in ccRCC tissues (n = 20). We found that the seven lncRNAs had a higher expression level in ccRCC tissues than in adjacent normal tissues (Figure 7A). Based on the lncRNA expression level in ccRCC tissues, we finally selected LINC01270 as the object of subsequent research. qRT-PCR showed that the expression level of LINC01270 was increased in ccRCC cell lines than in human renal tubular epithelial cells HK-2 cells (Figure 7B).

Figure 7 Experimental validation for the expression and function of LINC01270. (A) The differential expression level of seven lncRNAs between renal clear cell carcinoma tissues (n=20) and their paired normal tissues were detected by qRT-PCR. (B) The differential expression levels of LINC01270 between renal clear cell carcinoma cells and human renal tubular epithelial cells HK-2 cells were measured by qRT-PCR. (C) The transfected efficiency of LINC01270 was measured using qRT-PCR. (D and E) The CCK-8 and colony formation assays were used to detect cell proliferation ability. (F–I) The migration and invasion were used to measure the cell migration and invasion ability. *p<0.05; **p<0.01; ***p<0.001.

The Effect of LNC01270 in Proliferation, Invasion, and Migration of ccRCC Cell Lines

We further investigate the function of LINC01270 in ccRCC cell lines. The transfected efficiency of LINC01270 was verified. The expression of LINC01270 was significantly decreased in both ACHN and Caki-1 cells by siRNA transfection (Figure 7C). To explore whether LINC01270 plays an important role in ccRCC progression, CCK-8, colony formation assays, and Migration and Invasion assays were used. CCK-8 assays and colony formation assays showed that the cell proliferation of ccRCC cells was inhibited in si-LINC01270 group compared with the control group (Figure 7D and E). Moreover, results of the transwell suggested that LINC01270 knockdown inhibited migration and invasion of ccRCC cells (Figure 7F–I). Taken together, LINC01270 may act as an oncogene to promote ccRCC progression.

Discussion

This study comprehensively evaluated the prognostic potential of immune-related lncRNAs in ccRCC by leveraging the knowledge of MSigDB C7 immunologic signatures. First, we identified ccRCC-related lncRNAs using differential expression analysis. Next, the immune-related genes were identified by gene set enrichment analysis of immunologic signatures. The lncRNAs co-expressed with immune-related genes were considered immune-related lncRNAs and were used to perform an overall survival analysis. Finally, a 7-lncRNA signature was identified, which could be an unfavorable biomarker for the overall survival of late-stage ccRCC patients. Patients in the high-Pindex group showed poorer survival than those in the low-Pindex group. This signature provided a novel clinical tool for the prognosis of late-stage ccRCC at the molecular level. Importantly, we found that the LINC01270 was upregulated in ccRCC tissues and cell lines. Silence of LINC01270 resulted in inhibition of cell proliferation, invasion, and migration of ccRCC. Although current prediction model has significant limitations in clinical application, prognostic models provide information on which risk factors are associated with disease outcomes and the relative magnitude of their influence on disease outcomes. Moreover, prospective validation of this signature is necessary for future clinical applications.

The polygenic molecular signature is a composite signature that has shown strong ability in tumor diagnosis, prognosis, and classification of pathological subtypes, probably because it organically integrates the predictive roles of multiple key genes. An expression analysis using three independent cohorts showed that the seven lncRNAs were upregulated in ccRCC compared with normal controls. Based on the functional analysis, the 290 immune-related lncRNAs showed strong relevance to immune cell activity, including B cells, NK cells, and CD4/CD8 T cells. The seven prognostic lncRNAs were mainly associated with CD4/CD8 T cells activity. It was reported that RP11-37B2.1 was associated with thrombocytopenia development during anti-TB therapy and suggested that its genetic variants could be potential biomarkers for thrombocytopenia during anti-TB treatment.²¹ FIRRE can serve as a platform for trans-chromosomal co-localization in that the loci located on different chromosomes are brought together around the transcription site of FIRRE to modulate nuclear architecture across chromosomes.^22,23 Another study reported that FIRRE helped position the inactive X chromosome near the nucleolus and was related to the stabilization of H3K27me3 histone modifications.²⁴ Notably, it was revealed that FIRRE was involved in NF-kappaB signaling to positively regulate the expression of inflammatory genes in macrophages following LPS stimulation, suggesting that FIRRE was associated with the innate immune system.²⁵ There is no functional study for other lncRNAs to date. Thus, the 7-lncRNA signature could be a novel prognostic biomarker for the overall survival of ccRCC.

Cancer immune checkpoint blockade therapy is the most valuable cancer immunotherapy by blocking the activation of co-suppressive information molecules by T cells and releasing the anti-cancer immune response in a suppressed state to attack cancer cells. At present, the most effective immune checkpoint blockade drugs mainly target three genes: PD-1, PD-L1, and CTLA4. Their expression levels in tumors are closely related to the efficacy of this immunotherapy. We found that the seven lncRNAs were universally correlated with the three genes. Notably, RP11−504P24.9 and RP11−37B2.1 have a high correlation with CTLA4, implying that they may be of value in predicting the anti-CTLA4 immunotherapy.

Conclusions

In summary, we identified a 7-lncRNA signature (LINC01270, FIRRE, RP11-37B2.1, RP11-253I19.3, RP11-438L19.1, RP11-504P24.9, and CTB-41I6.1) with great potential in overall survival prognosis for patients with late-stage ccRCC, which could be a novel clinical tool at the molecular level.