The Cancer Genome Atlas (TCGA) database contains multi-omics datasets, including gene transcript, miRNA, and lncRNA expression data from numerous cancer types. We downloaded all KIRP datasets from TCGA (https://portal.gdc.cancer.gov/) using the TCGAbiolinks package in R software (version 3.5.1), including lncRNA, mRNA, and miRNA expression profiles from KIRP samples and the corresponding clinical follow-up data. Next, we used the ESTIMATE algorithm to output stromal scores (SSs) and immune scores (ISs) and compared the differentially expressed lncRNAs, mRNAs, and miRNAs in samples with high and low SSs and ISs (high >0 and low <0) (9).

We utilized the weighted gene co-expression network analysis (WGCNA) package in R to analyze the co-expression network of selected lncRNAs and mRNAs (10). The gene expression profile matrix was first transformed into a similarity matrix based on the Pearson test between pairwise genes. Next, the similarity matrix was converted into an adjacency matrix, and scale-free gene co-expression networks were constructed by the WGCNA package. The topological overlap matrix (TOM) and dissimilarity TOM (dissTOM) were then obtained by TOM similarity and dissimilarity on the basis of the adjacency matrix. Lastly, after hierarchical clustering analysis based on dissTOM, modules were generated from the Dynamic Tree Cut method for Branch Cutting.

Two independent cohorts of KIRP samples were acquired for this study from patients that underwent surgery at the Air Force Military Medical University Tangdu Hospital (Xi’an, China). From cohort A, 30 paired KIRP tissues and adjacent non-tumor tissues were gathered and used for reverse transcription (RT)-PCR and immunoblotting. From cohort B, 80 KIRP specimens were collected and used for constructing tissue microarrays and performing immunohistochemistry (IHC). This study was approved by the Ethics Committee on Human Research of Tangdu Hospital, and written informed consent was obtained from all patients.

Human full-length CAMK2B cDNA (NM_001220.5) was obtained from VectorBuilder (Guangzhou, China). This fragment was cloned into the pCDH lentiviral expression vector (System Biosciences, Palo Alto, CA, USA) and inserted between the XbaI and EcoRI sites, using the In-Fusion HD Cloning Kit (Takara, Tokyo, Japan). Plasmid constructs expressing CAMK2B shRNA and scramble shRNA were generated using the lentiviral expression plasmid PLKO.1 and obtained from Vectorbuilder. Lentiviral shRNA target sequences were as follows: CCGGAAGCAGGAGATCATTAA (sh1), GACCAGATGTGATTTGTTAAA (sh2), and ATAGAGGATGAAGACGCTAAA (sh3). Stable cell lines were created by antibiotic selection with puromycin for one week, beginning 72 h after transduction.

Total RNA was isolated from KIRP tissues and adjacent non-tumor tissues using TRIzol Reagent (Invitrogen, Thermo Fisher Scientific, Waltham, MA, USA). Isolated RNA was reverse transcribed using the PrimeScript RT Reagent Kit (TaKaRa, Shiga, Japan), and mRNA expression was determined by quantitative (q) qRT-PCR using SYBR Premix Ex Taq II (TaKaRa). Primers used for amplification of human genes were as follows: lncRNA GUSBP11 forward primer 5’-TCCCCTGTCCCGAAGGATTAC-3’ and reverse primer 5’ -TAAGGGACTAACGGCTTCGCT-3’; miR-432-5p forward primer 5’-ACTCAAACACTTCGGACATGG-3’ and reverse primer 5’-CAAAGAGCAACAGAGAGTAGCA-3’; CAMK2B forward primer 5’-CAGTGTACTTCAGTTGGTG-3’ and reverse primer 5’-TCAGGTTTTGCTCTTC TC-3’; GAPDH forward primer 5’-AGGTCGGAGTCAACGGATTTG-3’ and reverse primer 5’-TGTAAACCATGTAGTTGAGGTCA-3’.

Tumor tissue was fixed, embedded in paraffin, and sliced into 5-μm thick sections. IHC staining and immunoblotting were performed using primary antibodies to CAMK2B (1:100 dilution, Proteintech, Chicago, IL, USA), according to standard protocols, as previously described (11). Relevant bioinformatics methods, Culture conditions for cell lines, and the various functional assays for CAMK2B in vivo (animal models with subcutaneous xenografts) and in vitro (migration, proliferation, and antibodies and cell lines information) are described in the Supplementary Materials and Methods (Doc. S1) .

Statistical analyses were performed using SPSS 15.0 for Windows (SPSS, Inc., Chicago, IL, USA). IHC and qRT-PCR experiments were performed at least three times, and statistical quantitative data were evaluated by Student’s t-tests, with P <0.05 considered statistically significant.

We analyzed a patient cohort containing 233 KIRP patients for whom expression data were available in TCGA database. These data were analyzed using the ESTIMATE algorithm, as described in the flowchart in Figure 1A . Infiltrating cells in the tumor tissue were assessed by incorporating two gene signatures within the ESTIMATE algorithm. SSs were designed to capture the presence of stromal cells in tumor tissue, and ISs were calculated to determine the infiltration of immune cells in tumor tissue. Based on ISs, 233 samples were divided into an IS-high (178 samples) group and an IS-low (22 samples) group. We detected no difference in overall survival (OS) between patients in the IS-high and IS-low group ( Figure 1B, a). We then divided patients into SS-high (27 samples) and SS-low (206 samples) groups, and found a significant difference in survival between patients in these groups, with SS-high patients showing decreased survival relative to SS-low patients (P = 1×10^-4, Figure 1B, b). Based on these findings, we compared gene expression profiles in the SS-high and SS-low groups.

Identification of DEGs in KIRP samples with high vs. low SSs using the methods described above and in Figure 1A revealed 1668 differentially expressed mRNAs, 783 differentially expressed lncRNAs, and 58 differentially expressed miRNAs. All DEGs in each independent dataset were analyzed using unsupervised hierarchical clustering analysis. The cluster analysis heat map shows the correlation between expression maps and group conditions ( Figure 2A ). We further performed gene ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment analyses of 1189 upregulated mRNAs in SS-high subjects and identified immune response and cytokine-cytokine receptor interactions among the top pathways. Analysis of 479 downregulated mRNAs yielded oxidation-reduction process and metabolic pathways as top hits ( Supplementary Figure 1 ).

We next performed sample clustering detection of lncRNAs and mRNAs and found no outlier samples, indicating it was not necessary to remove any samples in our co-expression analyses ( Figure 2B ). The most important step in this process is selection of the value for soft-threshold power (β). To determine the relative equilibrium between scale independence and mean connectivity, we analyzed the network topology with soft-threshold power from 1 to 20, eventually confirming a β value of 8 for lncRNAs and mRNAs ( Figure 2C ). Next, we identified co-expression modules by hierarchical clustering and Dynamic Branch Cutting and assigned a unique color as an identifier in each module using WGCNA ( Figure 2D ).

To identify interactions among these co-expression modules, we analyzed the connectivity of eigengenes in a cluster analysis. Seven modules (turquoise, blue, grey, brown, yellow, green, red) were generated, with the grey module representing a gene set that was not assigned to any of the modules. We then generated a heatmap plot of the adjacencies in the hub gene network, wherein red represents positive correlation with high adjacency, and blue represents negative correlation with low adjacency ( Figure 2E ). The eigengene dendrogram and heatmap were used to identify groups of correlated eigengenes, and the dendrogram revealed a significant association between these modules and KIRP clinical traits. Notably, the turquoise module was found to be negatively correlated with pathological stage of the tumor and was chosen as the focus of research for this study ( Figure 2F ). This module contained a total of 1004 differentially expressed RNAs, including 364 mRNAs, 600 lncRNAs, and 40 other noncoding RNAs. GO enrichment analyses of the 364 mRNAs in the turquoise module identified negative regulation of apoptotic process as a top function, with the largest number of enriched genes ( Figure 2G, a), and KEGG-pathway enrichment analyses found metabolic pathway as a key mechanism ( Figure 2G, b).

Through our predefined strategy based on SSs, DEGs in the turquoise module were found to be mainly associated with favorable survival outcomes in KIRP patients. We then analyzed co-expression of ceRNA networks in the turquoise module by WGCNA and built multiple groups, including RP3-416H24.1/hsa-miR-127-5p/ITGB8, SNHG11/hsa-miR-214-3p/ARVCF, SNHG11/hsa-miR-214-3p/HDAC11, SNHG11/hsa-miR-214-3p/KAZN, SNHG11/hsa-miR-214-3p/MRPS25, SNHG11/hsa-miR-214-3p/PRR15L, and GUSBP11/miR-432-5p/CAMK2B ( Figure 3A ).

The prognostic effects of these DEGs were analyzed and presented in a forest plot, and eight lncRNAs, including GUSBP11, IQCH-AS1, MTHFD2P1, RP11-12A2.3, RP3-416H24.1, SNHG11, TMEM51-AS1, and RP11-567M16.6, were found to exhibit pro-survival functions with longer overall survival (OS) ( Figure 3B, a). A total of seven miRNAs, including mir-432-5p, mir-145-5p, mir-143-5p, mir-127-5p, mir-214-3p, and mir-134-5p, were correlated with a poor prognosis and shorter OS ( Figure 3B, b). In addition, 15 mRNAs were correlated with better prognosis, including CAMK2B, MACC1, PKHD1, PTPRD, RAB36, ZNF320, FLRT3, ITGB8, ARVCF, HDAC11, KAZN, MRPS25, P4HTM, PAX8, and DACT2 ( Figure 3B, c). We then plotted Kaplan–Meier curves using the Cox proportional hazard regression model for SS-related lncRNAs, miRNAs, and mRNAs in KIRP and confirmed that all these genes are associated with survival ( Supplementary Figure 2 ). Based on results of correlation strength and survival analyses, CAMK2B was selected as one of the most promising target genes.

We then assessed the prognostic impact of CAMK2B in 33 types of human tumors using GEPIA2 (12), and found that high expression of CAMK2B is associated with improved prognosis in KIRP, liver hepatocellular carcinoma (LIHC), and ovarian cancer (OV, Figure 3C, a). To further confirm this finding, we used Kaplan–Meier plotter and found that that high expression of CAMK2B is associated with longer OS and recurrence-free survival (RFS) in KIRP ( Figure 3C, b).

We next measured expression of the selected ceRNA network in KIRP tissue by qRT-PCR. We detected lower levels of lncRNA GUSBP11 (P=0.0618) and CAMK2B (p=0.0447) in tumor tissue samples than in the normal tissue adjacent to the carcinoma, whereas levels of miR-432-5p were found to be higher in tumor tissue (P=0.0097; Figure 4A ). Upon further exploration, we found that high expression of lncRNA GUSBP11 is associated with no metastasis (P=0.0023; Figure 4B, a), early-stage disease (P=0.0023; Figure 4B, b), and small tumor size (P=0.0249; Figure 4B, c). Similarly, high expression of CAMK2B mRNA is associated with no metastasis (P = 0.002; Figure 4D, a), early-stage disease (P=0.0299; Figure 4D, b), and small tumor size (P=0.0226; Figure 4D, c). In contrast, high expression of mir-432-5p is associated with tumor metastasis (P=0.0037; Figure 4C, a), late-stage disease (P=0.007; Figure 4C, b), and large tumor size (P=0.027; Figure 4C, c).

We further confirmed these findings by analyzing KIRP tumor samples and adjacent non-tumor tissues by IHC using tissue microarray ( Figure 4E, F ). Indeed, we found that protein levels of CAMK2B are lower in tumor tissue than in paired tumor-adjacent tissue (P=0.0227, Figure 4F, a, b). Further, expression of CAMK2B in tumors <3 cm in size is significantly higher than in tumors >3 cm in size (P=0.0005, Figure 4F, c, Confirm whether the insertion of the Ethics Statement section is fine. Note that we have used the statement provided at Submission. If this is not the latest version, please let us know.). High expression of CAMK2B was also found to be associated with lower tumor grade (P=0.0328, Figure 4F, d) and low metastasis (P=0.0375, Figure 4F, e). Thus, these data confirm the expected pattern of CAMK2B expression and suggest it acts as a protective factor in KIRP.

To evaluate the function of CAMK2B, we stably overexpressed and silenced CAMK2B expression in the KIRP SK-RC-39 cell line ( Figure 5A ). We found that upregulation of CAMK2B significantly inhibits cell proliferation ( Figure 5B, a), whereas the reverse occurs in response to CAMK2B silencing ( Figure 5B, b). Moreover, we confirmed the inhibition of proliferation in response to CAMK2B overexpression using sphere formation assays, and showed that CAMK2B silencing has the opposite effect ( Figure 5C ). Migration assays further showed that a lower number of CAMK2B overexpressing cells cross the basement membrane relative to SK-RC-39 control cells, whereas migration is enhanced after CAMK2B silencing ( Figure 5D ).

In nude mouse models, diminished subcutaneous tumor growth was observed in mice injected with SK-RC-39 cells stably overexpressing CAMK2B (P=0.0351, Figure 5E ) as compared to mice injected with SK-RC-39 control cells. Immunohistochemical staining further revealed that SK-RC-39-CAMK2B subcutaneous tumors from these animals show significantly decreased expression of CD34 (P=0.0363), α-SMA (P=0.0451), and PCNA (P=0.0177, Figure 5F ) relative to SK-RC-39 tumors. We also confirmed this negative correlation between expression of CAMK2B and both α-SMA (R=-0.26, Figure 5G, a) and CD34 (R=-0.35, Figure 5G, b) in human KIRP tissue using the GEPIA2 database (12), indicating that elevated CAMK2B levels are associated with decreased infiltrating fibroblasts and vascular endothelial cells. Additionally, the negative correlation between expression of CAMK2B and endothelial cells (R=-0.455) and cancer-associated fibroblast (R=-0.384) was confirmed using Xcell and EPIC. In contrast, only a weak correlation is observed between CAMK2B and immune cells ( Figure 5H ). Thus, the above findings suggest that upregulation of CAMK2B decreases proliferation and inhibits tumor stromal cell infiltration in KIRP.

In order to identify potential target genes downstream of CAMK2B, we screened for CAMK2B co-expressed genes in KIRP from TCGA database using LinkedOmics (13) and Metascape (14). We then constructed a volcano plot to visualize genes that are positively and negatively correlated with CAMK2B and identified CHL1 as the most positively correlated gene and VEGFA as the most negatively correlated gene ( Figure 6A ). Using Metascape, we found that the top 100 genes correlated with CAMK2B are significantly enriched in cellular component organization and biogenesis ( Supplementary Figure 3 ). To confirm the relationship between CAMK2B, VEGFA, and CHL1 in KIRP, we measured the levels of these proteins in frozen tissue samples from KIRP patients by immunoblotting, revealing a positive correlation between CAMK2B and CHL1 and a negative correlation between CAMK2B and VEGF ( Figure 6B ). Further evaluating their expression in KIRP tissue using GEPIA2 confirmed a significant positive correlation between expression levels of CAMK2B and CHL1 (R=0.7433), as well as a negative correlation between expression of CAMK2B and both VEGF (R=-0.3) and TGFβ1 (R=-0.23, Figure 6C ). Finally, we examined expression of several effector molecules in SK-RC-39 cells stably overexpressing CAMK2B or subjected to shRNA-mediated CAMK2B silencing. Notably, we found that expression of CHL1 and E-cadherin are significantly increased, whereas TGFβ1, VEGF, vimentin, and PCNA are significantly downregulated, in SK-RC-39 cells overexpressing CAMK2B. Conversely, when CAMK2B expression in SK-RC-39 cells is silenced using specific shRNA, the opposite expression pattern is observed ( Figure 6D ). Thus, we identify CHL1, VEGF, and TGFβ1 as downstream effectors of CAMK2B, and provide further evidence that enhanced CAMK2B expression inhibits proliferation and stromal cell infiltration.