Identification of Critical Pathways and Hub Genes in LanCL1-Overexpressed Prostate Cancer Cells

Introduction

Prostate cancer remains the most common cause of cancer-related death in men, particularly in developed countries.¹ Up to now, many studies have found out multiple critical genes, cellular processes, and signaling pathways that contribute to disease progression. In previous study, we focused on LanCL1 and explored whether it contributes to prostate cancer disease progression. As a result, we determined that LanCL1 expression correlated with the progression of human prostate cancer and that LanCL1 could also promote prostate cancer proliferation. Moreover, LanCL1 protects prostate cancer cells from reactive oxygen species (ROS) via suppression of JNK pathway.²

Lanthionine synthase C-like protein 1 (also known as P40 or GRP69A) is a mammalian member of the LanC-like protein superfamily. It contains highly differentiated peptide-modifying enzymes found in bacteria and plants (LanCs).³ Multiple studies reveal the role of LanCL1 in relieving oxidative stress in neurons during normal development by binding zinc ion and GSH. In addition, LanCL1 can catalyze the formation of thioether products in neurons, thereby protecting them from oxidative stress.^4–6 In breast cancer, LanCL1 expression is related to breast cancer survival, making it a marker of aging.⁷ In prostate cancer, we have preliminarily explored the role played by LanCL1 in the disease process. However, further mechanisms are still unclear. As we all know, a gene or protein usually exerts its function by affecting the signaling pathways downstream. In this study, we attempt to further study and explore the mechanism underlying to discover new potential targets for disease diagnosis and treatment.

Materials and Methods

Constructs

We generated pPB-CAG-ires-Pac vector as previously described.^2,8 We inserted full-length LanCL1 into the multiple-cloning sites (MCS) of pPB-CAG-ires-Pac to generate pPB-CAG-LanCL1-ires-Pac vector.

Cell Lines and Cell Culture

LNCaP and PC-3 cells (ATCC) were maintained in RPMI1640 supplemented with 10% FBS. LNCaP and PC-3 cells were co-transfected with pPB-CAG-LanCL1 and pCMVPBase to obtain control and LanCL1 stable overexpressing LNCaP/PC-3 cells. We used 2 μg/mL puromycin (Amresco) to screen for 2 weeks for selecting the stable cell lines. Stable cell lines were identified by Western blotting. All cells were added with antibiotic-antifungal solution (100 units/mL penicillin, 0.1 mg/mL streptomycin, and 0.25 mg/mL amphotericin B), and under standard cell culture conditions (5% CO2, 95% humidity) grow at 37°C.

RNA Extraction

According to the manufacturer’s instructions (Invitrogen), total RNA was extracted from the tissue using TRIzol^® reagent and genomic DNA was removed using DNase I (TaKara). RNA quality was then determined by 2100 Bioanalyser (Agilent) and quantified using ND-2000 (NanoDrop Technologies). Only high-quality RNA samples (OD260/280 = 1.8–2.2, OD260/230 ≥ 2.0, RIN ≥ 6.5, 28S: 18S ≥ 1.0, >10 μg) were used to construct the sequencing library.

Library Preparation, and Illumina Hiseq4000 Sequencing

RNA-seq was performed according to the manufacturer’s guidelines (C-10,365; Life Technologies) and previous research.⁹ RNA-seq transcriptome library was prepared following TruSeq^TM RNA sample preparation Kit from Illumina (San Diego, CA) using 5μg of total RNA. Libraries were size selected for cDNA target fragments of 200–300 bp on 2% Low Range Ultra Agarose followed by PCR amplified using Phusion DNA polymerase (NEB) for 15 PCR cycles. After quantified by TBS380, paired-end RNA-seq sequencing library was sequenced with the Illumina HiSeq 4000 (2 × 150 bp read length). RNA-seq was performed on three biological replicates.

Analysis of RNA-Seq and Differentially Expressed Genes (DEG)

RNA-seq was performed according to the manufacturer’s guidelines (C-10,365; Life Technologies) and previous research.⁹ We then used EdgeR, an R package to examine differential expressed genes as previously described.^10,11 DEG was determined using the following criteria: | Fold Change (FC) | ≥ 2 and P value < 0.05. DEG is used for further bioinformatics analysis.

Gene Set Enrichment Analysis (GSEA)

GSEA v3.0 was used to analyze the mRNA expression level of biological function annotations and pathways after LanCL1 overexpression. The results indicate the effects of LanCL1 overexpression on various biological functional gene sets. The number of permutations is set to 100. In the analysis, we used the following values as cut-off values: the P-value cut-off value is 0.05 and the false discovery rate (FDR q-val) is <0.25.

Functional Annotation and Pathway Enrichment Analysis

We uploaded the top 2000 DEGs and the top 2000 downgraded DEGs to the database, respectively, for annotation, visualization, and integrated discovery (DAVID) website for GO analysis and KEGG analysis,¹² while specifying a P-value <0.05 for statistical significance.

Integration of Protein–Protein Interaction (PPI) Networks and Module Analysis

We used STRING database, which provides a critical assessment and integration of PPI, to analyze the protein interaction.¹³ To evaluate the interactive relationships among the DEGs above, we uploaded the top 1000 upregulated and top 1000 downregulated DEGs together to STRING. Then, only experimentally validated interactions (total score > 0.4) were selected as meaningful. The PPI network was then used for module screening by Molecular Complex Detection (MCODE) in Cytoscape (scores > 3, nodes > 4).¹⁴

Publicly Available Gene Expression and Clinical Datasets for Prognosis Analysis

We download the gene expression data with corresponding prognostic profiles of 140 prostate cancer patients from online dataset GSE21032 for prognosis analysis.¹⁵ The Kaplan–Meier method using Log-rank test was used to calculate the clinical results between different gene expression groups by Graphpad. Values of P <0.05 were considered statistically significant.

Statistical Analysis

The Kaplan–Meier method using Log-rank test was used to calculate the clinical results by Graphpad. Multivariate Cox regression analyses were used to investigate the association between prognosis and hub genes expression/clinical factors. The FDR in GSEA was tested several times using the Benjamini–Hochberg program to control the FDR.^16,17 Values of P <0.05 were considered statistically significant. All statistical analyses were performed using SPSS v19, Graphpad, and R 3.3.0.

Results

RNA-Seq and Identification of DEGs

To further explore the critical signaling pathways or the regulatory network downstream, we obtained the LanCL1 stably overexpressed LNCaP cell line. Western blotting showed upregulation of LanCL1 in stable cell lines (Figure 1A). We then performed RNA-seq on a stable cell line with LanCL1 overexpression or an empty vector as a control. RNA-seq was performed on three biological replicates.

Based on computer analysis, we used P < 0.05 and Fold Control (FC)≥ 2.0 or ≤0.5 as criteria, and identified 8928 genes as DEGs, of which 4755 were upregulated and 4173 were downregulated. The volcano map of DEGs is shown in Figure 1B, while the heat map of DEGs expression for the first 100 genes is shown in Figure 1C.

GO and KEGG Analyses of DEGs

To further analyze the DEGs at the functional level, we submitted the top 2000 upregulated and top 2000 downregulated DEGs, respectively, online for further analyses using DAVID. In the upregulation group, the GO-BP (Biological Process) analysis of DEGs (Figure 1D) suggested significant enrichment in extracellular matrix organization, collagen fibril organization, proximal/distal pattern formation, cell adhesion, peptidyl-tyrosine phosphorylation, axon guidance, anterior/posterior pattern specification, chemical synaptic transmission, regulation of ion transmembrane transport, ureteric bud morphogenesis, angiogenesis, negative regulation of neuronal apoptosis, ureteric bud development, calcium ion-regulated exocytosis of neurotransmitter, and positive regulation of synapse assembly. We also showed the results of GO-CC (Cellular Component) and GO-MF (Molecular Function) in Supplementary Figures 1A and B. While KEGG pathway analysis showed enrichment mainly in focal adhesion, ECM–receptor interaction, PI3K-AKT signaling pathway, arrhythmogenic right ventricular cardiomyopathy (ARVC), axon guidance, pathways in cancer, Ras signaling pathway, cell adhesion molecules (CAMs), amoebiasis, cAMP signaling pathway, hypertrophic cardiomyopathy (HCM), protein digestion and absorption, proteoglycans in cancer, platelet activation, and Rap1 signaling pathway (Figure 1E).

In downregulation group, in GO-BP (Biological Process) analysis, DEGs were significantly enriched in oxidation-reduction process, homophilic cell adhesion via plasma membrane adhesion molecules, inflammatory response, establishment of skin barrier, nucleosome assembly, phosphatidylserine acyl-chain remodeling, acute-phase response, angiogenesis, phosphatidylinositol acyl-chain remodeling, positive regulation of cell adhesion, phosphatidylethanolamine acyl-chain remodeling, phospholipid metabolic process, negative regulation of viral genome replication, positive regulation of inflammatory response, and transmembrane transport (Figure 1F). The results of GO-CC (Cellular Component) and GO-MF (Molecular Function) were shown in Supplementary Figures 1C and D. KEGG pathway analysis showed enrichment mainly in arachidonic acid metabolism, alcoholism, systemic lupus erythematosus, linoleic acid metabolism, alpha-Linolenic acid metabolism, drug metabolism – cytochrome P450, chemical carcinogenesis, influenza A, glutathione metabolism, mineral absorption, ether lipid metabolism, transcriptional misregulation in cancer, Aldosterone-regulated sodium reabsorption, Hepatitis C, and Bladder cancer (Figure 1G).

GSEA

To determine the effect of LanCL1 overexpression on prostate cancer progression, we then analyzed various gene sets by GSEA. As shown in Figure 2 and Supplementary Figure 2, the GSEA results indicate biological processes including UV response to DN, WNT beta-catenin signaling, Hedgehog signaling, hypoxia, myogenesis, angiogenesis, apical junction, KRAS signaling, complement, IL2-STAT5 signaling, epithelial–mesenchymal transition, allograft rejection, pancreas beta cells, and mitotic spindle were significantly enriched.

PPI Network Analysis

To investigate the protein–protein interaction and find out the hub genes, we then uploaded the top 1000 upregulated DEGs and 1000 downregulated DEGs to the website to screen the information in the STRING database. The top 15 genes ranked by degree were identified as hub genes. These hub genes included TNF, FGF2, BMP4, CXCL8, BDNF, ESR1, IGF1, KDR, ICAM1, MMP2, KNG1, NCAM1, CDH2, KIT, and APOE (Table 1). TNF is the top hub gene with the highest degree of 137. Then, we identified the gene modules in the PPI network through MCODE. The top three modules were selected, and we performed GO-BP and KEGG pathway enrichment analysis (Figure 3). Functional annotations in module 1 genes are richer in G protein-coupled receptor signaling pathways, chemotaxis, positive regulation of cytosolic calcium ion concentration, chemokine signaling pathways, neuroactive ligand–receptor interactions, and Cytokine–cytokine receptor interaction.

Table 1 Functional Roles of the 15 Hub Genes

Enrichment analysis showed that the genes in module 2 were mainly associated with interferon-gamma-mediated signaling pathway, type I interferon signaling pathway, collagen catabolic process, protein digestion and absorption, ECM–receptor interaction, and herpes simplex infection. In module 3, genes were mainly enriched in extracellular matrix organization, nucleosome assembly, integrin-mediated signaling pathway, systemic lupus erythematosus, alcoholism, and hypertrophic cardiomyopathy (HCM).

Clinical Outcome and Hub Genes

To select the hub genes that are closely related to clinical prognosis, we use the data and clinical information of GSE21032.¹⁵ Figure 4 and Table 2 show that 4 of the 15 hub genes strongly associated with clinical outcomes, including APOE, FGF2, ICAM1, and KIT. Lower APOE expression level correlated with better outcome, while higher FGF2, ICAM1, and KIT correlated with better outcome. Then, we reorganized the clinical data of GSE21032 and found that it contains two clinical factors, Gleason score and Tumor stage (Table 3). To investigate the association between prognosis and these clinical two factors, and to evaluate whether the expressions of the 4 hub genes above are independently associated with disease recurrence, we first did Kaplan-Meier analysis. The results suggest that both Gleason score and Tumor stage significantly correlated with the prognosis of the disease (Table 2). Then, we conducted a multivariate COX regression analysis. Results showed that only the expression of ICAM1 in the four genes, as well as Gleason score and T stage, is an independent risk factor for the prognosis (Table 4). At last, to confirm our results, we investigated the expression of the 4 critical hub genes after LanCL1 overexpression in another prostate cancer cell line, PC-3, which revealed the same results as in LNCaP cells (Figure 5). These results suggested that these 4 genes, especially ICAM1, might be potential key genes for further research.

Table 2 Disease Free Survival and Clinicopathological Affairs/Hub Genes Expression

Table 3 Summary of the Clinical Characteristics of Patients with Prostate Cancer in GSE21032

Table 4 Multivariate Cox Regression Analyses Investigating the Association Between Prognosis and Hub Genes Expression/Clinical Factors

Discussion

Prostate cancer is still one of the most common malignant tumors in male urinary system all around the world, especially in developed countries. In-depth research on the downstream pathways and regulatory networks of key genes involved in disease progression is of great significance for understanding the disease and finding new therapeutic targets. In our previous study, we revealed the role of LanCL1 as an oncogene in prostate cancer to promote cell proliferation and protect cell from attack of ROS via JNK pathway. However, important downstream genes and regulatory signaling pathways are still not clear.

Here in this study, we overexpressed LanCL1 gene in LNCaP cells and did RNA-Seq analysis. By selecting DEGs and performing GO and KEGG analysis on the top 4000 DEGs, we found that LanCL1 overexpression was significantly associated with the metabolism, axon and synaptic function, and several signaling targets or pathways, such as cAMP signaling pathway, Ras signaling pathway, PI3K-AKT-mTOR signaling pathway, and Rap1 signaling pathway. Previous researches have shown that Rap1 (ras-related protein 1) signaling pathway and Ras signaling pathway are related to JNK and participate in cellular oxidative stress and the production and elimination of ROS in a variety of cellular backgrounds.^18–21 In addition, the cAMP signaling pathway is also involved in the production, metabolism, and clearance of ROS in the body, which has the same important role as the JNK-related pathway.^22,23 These results suggest that LanCL1 plays an important role in assisting cells in processing ROS, which is consistent with our previous findings. Besides, because the oxidative stress state and related signaling pathways in the cells will change after H₂O₂ treatment, there may be significant differences in the signaling pathways between treated and untreated cells. This might be the main reason why JNK pathway was not enriched in the RNA-seq results of untreated cells in this study. JNK pathway is the main mechanism by which LanCL1 protects prostate cancer cells from oxidate stress in H₂O₂ treated cells. Our results here also show that the mechanism of LanCL1 in dealing with ROS is more complicated. It also validates our previous experiments, that is, changing the activity of JNK pathway in rescue experiments only partially restored the ROS-related phenotype.

In addition to oxidative stress such as ROS, we also found that many pathways related to tumorigenesis and progression appeared enriched after LanCL1 overexpression, mainly including WNT beta-catenin signaling, PI3K-AKT signaling pathway, and KRAS signaling. These pathways play important roles in the pathogenesis and progression of various tumors, including prostate cancer.

In prostate cancer, WNT/β-catenin and AR signals inhibit each other. β-catenin activation and PTEN deletion synergistically drive non-AR independent castration-resistant prostate cancer.^24–26 In most cancer types, the PI3K signaling pathway is activated and promotes tumorigenesis by regulating nutrient metabolism, cell proliferation, survival, migration, and angiogenesis.²⁶ The potential mechanism of PI3K/AKT activation is mainly due to the deletion or mutation of its key negative regulatory gene, PTEN.²⁷ In fact, it is known that the mutual interference between PI3K/AKT and other pathways (such as AR, WNT, and ERK signaling pathways) plays a vital role in the disease progression and drug resistance of prostate cancer, making co-targeting the PI3K/AKT signaling pathway and its synergistic pathways critical to improve the anticancer efficacy of PI3K/AKT inhibitors.^28,29 Our results indicate that LanCL1 may contribute to prostate cancer by influencing these important tumor regulatory signal pathways or signal networks, especially its promotion effect on prostate cancer cell proliferation. LanCL1 also may affect one of these pathways. As for the specific regulatory mechanism, further molecular biology experiments are needed to verify it.

In the PPI network analysis, we identified the top 15 central genes with the highest degree of interaction. TNF (tumor necrosis factor) is the top hub gene, which is mainly secreted by macrophages. TNF can bind to its receptors TNFRSF1A/TNFR1 and TNFRSF1B/TNFBR to exert its function. TNF could involve in the development of tumors by regulating inflammatory and immune responses. Among the 15 hub genes, we also identified 4 hub genes that strongly associated with clinical outcomes in prostate cancer, including APOE, FGF2, ICAM1, and KIT. We used data from a classic and authoritative prostate cancer RNA-seq database GSE21032 to analyze. These 4 genes might play critical roles in the development of prostate cancer and are closely related to the prognosis of the disease. Then, we conducted a multivariate COX regression analysis and found that only the expression of ICAM1 is an independent risk factor for the prognosis. This result emphasizes the important potential role of ICAM1 in prostate cancer initiation and development. ICAM1 (Intercellular adhesion molecule 1) is a cell surface glycoprotein expressed on endothelial cells and immune system cells. It binds to CD11a/CD18 or CD11b/CD18 type integrins and is also used by rhinovirus as a receptor.^30,31 Previous studies have shown that ICAM1 may facilitate breast and lung.^32–34 Researchers reported that CCN3-ICAM1 and MiRNA-296-3p-ICAM1 axes could influence cell migration in prostate cancer cells.^35,36 However, the role and mechanism of ICAM1 in prostate cancer are still largely unknown. ICAM1 may be further studied as key a gene to find more clues about the initiation and progression of prostate cancer disease. As for APOE, FGF2, and KIT, many previous studies have investigated their roles in prostate cancer, and our results further confirmed these findings.

Our study contained one limitation. In our study, we analyzed the data of RNA-seq to identify the DEGs after LanCL1 overexpression and figured out the critical genes and pathways that might play key roles after LanCL1 overexpression in prostate cancer using bioinformatics analysis. The mechanism and validation of LanCL1 overexpression in prostate cancer still need further research in clinical and molecular biology experiments.

Conclusions

In conclusion, our results revealed the critical genes and pathways that might play key roles after LanCL1 overexpression in prostate cancer. We also provide 4 critical hub genes that could play critical roles in the development of prostate cancer and are closely related to the prognosis of the disease, which may help in developing prostate cancer treatment strategies.

Figure 1 Differentially expressed genes (DEGs) and functional enrichment results of differentially expressed genes. DEG was determined using the following criteria: | Fold Change (FC) | ≥ 2 and P value <0.05. For enrichment analysis, a P-value <0.05 is taken to indicate statistical significance. (A) Western blotting result of LanCL1 overexpression in LNCaP cells. (B) Volcano plot for differentially expressed genes. (C) Heat map of the top 100 differentially expressed genes (50 upregulated genes and 50 downregulated genes). Red: upregulation; purple: downregulation. (D) The GO- BP (Biological Process) enrichment terms of upregulated genes. (E) The KEGG pathway analysis of upregulated genes. (F) The GO- BP (Biological Process) enrichment terms of downregulated genes. (G) The KEGG pathway analysis of downregulated genes.

Figure 2 Gene set enrichment analysis (GSEA) results of LanCL1 overexpression in prostate cancer cells (Top1–9). In the analysis, we used the following values as cut-off values: the P-value cut-off value is 0.05 and the false discovery rate (FDR q-val) is <0.25.

Figure 3 Top three modules from the protein–protein interaction (PPI) network. (A) PPI network of module 1. (B) GO-BP (Biological Process) and KEGG analyses of module 1. (C) PPI network of module 2. (D) GO-BP (Biological Process) and KEGG analyses of module 2. (E) PPI network of module 3. (F) GO-BP (Biological Process) and KEGG analyses of module 3.

Figure 4 Prognosis analyses of the hub genes. (A–D) Four hub genes that significantly related to prognosis by analyzing the data from GSE21032. The Kaplan-Meier method using Log-rank test was used to calculate the clinical results. Values of P <0.05 were considered statistically significant. Abbreviations: APOE, apolipoprotein E; FGF2, fibroblast growth factor 2; ICAM1, intercellular adhesion molecule 1; KIT, KIT proto-oncogene.

Figure 5 The expression levels of the 4 hub genes in PC-3 LanCL1 overexpression cells. (A–D) The mRNA expression levels of Apolipoprotein E (APOE, (A)), Fibroblast Growth Factor 2 (FGF2, (B)), Intercellular Adhesion Molecule 1 (ICAM1, (C)) and KIT Proto-Oncogene (KIT, (D)) in LanCL1 overexpressed PC-3 cells. N=5 and a P-value <0.05 is taken to indicate statistical significance. **P<0.01, ***P<0.001.