Study on the role of SLC14A1 gene in biochemical recurrence of prostate cancer

Prostate cancer (PCa) is a common malignant disease among men and biochemical recurrence (BCR) is considered to be a decisive risk factor for clinical recurrence and PCa metastasis. Clarifying the genes related to BCR and its possible pathways is vital for providing diagnosis and treatment methods to delay the progress of BCR. An analysis of data concerning PCa from previous datasets of The Cancer Genome Atlas (TCGA) and the Gene Expression Omnibus (GEO) was performed. Immunohistochemical (IHC) staining were used to evaluate the expression of SLC14A1 in prostate tissues. Kaplan–Meier analysis, Pearson correlation, and single sample Gene Set Enrichment Analysis (ssGSEA) were used to identify the potential pathway and molecular mechanism of the function of SLC14A1 in BCR of PCa. The expression of SLC14A1 is significantly reduced in prostate cancer cells and tissue comparing to normal prostate epithelial cell and para-cancerous tissue. As indicated by Kaplan–Meier analysis, High expression of SLC14A1 could increase the BCR-free survival time of PCa patients. This effect might be related to the interaction with miRNAs (has-miR-508, has-mir-514a2, and has-mir-449a) and the infiltration of B cells. SLC14A1 is a novel important gene associated with BCR of PCa, and further studies of its molecular mechanism may delay the progress of BCR.

Prostate cancer (PCa) is a common malignant disease and causes second-ranked cancer-related death among men 1 . The biochemical recurrence (BCR) of PCa is an important indicator of the prognosis of patients, and it often indicates the recurrence of the tumor. According to the guidelines of the European Association of Urology and the American Urological Association, BCR is defined as consecutive prostate-specific antigen (PSA) values no less than 0.2 ng/mL following radical prostatectomy (RP) 2,3 . It is worth pointing out that BCR does not mean clinical recurrence, but it is a decisive risk factor for PCa-specific mortality and overall mortality 4 . As the disease progressed without a second treatment after the diagnosis of BCR, 30% of patients have the median survival period of 5-8 years and around 32-45% among these patients would suffer PCa-specific mortality within 15 years 5 . Also, the clinically localized PCa, which is a significant cancer-specific survival benefit in patients after RP, still has nearly 30% PCa patients progressing into BCR after surgery 6,7 .
Similar to other types of tumors, the molecular mechanism of PCa development and patient prognosis is still unclear. Therefore, digging out the related genes of BCR in PCa can be beneficial to disease monitoring and treatment after RP. The human solute carrier family 14 member 1 (SLC14A1) gene, which is located on chromosome 18q12.1-21.1, contains about 30 kb nucleotides 8 . It is a solute carrier family gene that attracted people's attention by genome wide association studies (GWAS) and was considered to be related to the field of urothelial cancer 9 . The recent studies only indicate that SLC14A1 expresses differently between PCa and benign prostate tissue. SLC14A1 was found down-regulated in PCa tissue at around 2.88 fold and the castration will increase the expression of SLC14A1 by 3.05 fold, indicating that the expression of SLC14A1 gene in the prostate could be regulated by androgen, but the specific relationship of SLC14A1 and PCa has not been reported previously 10 .
In our study, we aimed to: (I) identify the relationship of SLC14A1 and PCa and verify the correlation between SLC14A1 and BCR using data obtained from the Gene Expression Omnibus (GEO) and The Cancer Genome Atlas (TCGA) database; (II) predict the correlation between SLC14A1 and miRNAs, and find the miRNAs which are related to BCR using Pearson correlation analysis and Kaplan-Meier survival curves; (III) analyze Determination and verification of the cut-off value of SLC14A1 gene expression. A total of 429 samples with BCR follow-up time in TCGA (experimental group) were included in the study, and 106 samples with BCR follow-up time in GEO (validation group) were included in the study. The cut-off value SLC14A1 gene expression was calculated by "surv_cutpoint" using the "survminer" package. In the experimental group, patients were be divided into 2 groups based on their SLC14A1 gene expression, which are high expression group and low expression group. Differential analysis was performed using the "edgeR" package. The same cutoff value was also used in the validation group for verification. Kaplan-Meier analysis with log-rank test was used to assess the difference in BCR-free survival between the high expression group and low expression groups with the help of the "survival" package in both experimental group and validation group.
Immunohistochemical staining. Tissue chip containing 50 prostate cancer patients (including cancer and paracancerous tissues) was purchased from Guangzhou Wo Zhao Biotechnology Co., Ltd. Immunohistochemical staining (IHC) was performed to evaluate the expression of SLC14A1 of the prostate cancer tissue comparing with the prostate paracancerous tissue. Briefly, slide mounted sections were brought to room temperature and dried for 30 min. Slides were then placed in an antigen retrieval solution, which contained 10 mM sodium citrate in water (pH 6. Differential analysis of miRNAs based after grouping. Differential analysis was performed based on the SLC14A1 gene expression using the "edgeR" package, all the miRNAs profiles were downloaded using the "TCGA biolinks" package. Statistically different miRNAs were screened out for correlation analysis (Pearson correlation) to assessed the difference in miRNA expression between tissues with low and high SLC14A1 expression. In addition, Kaplan-Meier analysis was performed, too. All included genes have FDR < 0.01.

Functional enrichment analysis and gene set enrichment analysis.
Functional enrichment analysis and gene set enrichment analysis was conducted to identify the biological processes, cellular components, and molecular functions of SLC14A1 related to BCR using the "clusterprofiler" package, based on TCGA differential analysis and their multiples of difference.
Differences in clinical phenotypes. The clinical data were also obtained from GSE54460, the clinical phenotypes comparison between the high expression and low expression group was performed to identify the relationship of SLC14A1 and clinical phenotypes.

Differences in immune cells after grouping. One of the biggest advantages of ssGSEA (single sample
Gene Set Enrichment Analysis) to quantify immune cell infiltration is that researchers could customize and quantify the types of immune infiltrating cells with the help of gene markers. The information of 24 immune cells was obtained from the currently recognized and used most immune cell markers 16 . Identify the immune cells of TCGA-PRAD and GSE46602 with ssGSEA, and Wilcoxon test were performed to compare the differences between groups of immune cells in different groups. In addition, multivariate cox regression analysis was also performed to screen out independent risk factors for BCR events.
Statistical analysis. All

Results
Four expression profiles of SCL14A1 in normal and PCa tissue were obtained from GEO, and they were GSE32448, GSE46602, GSE69223, and GSE70768, respectively. The expression of SCL14A1 was significantly higher in normal prostate tissue compared to tumor tissue, the p values were 7.9 e−09 , 1.5 e−11 , 2.7 e−05 , and < 2.2 e−16 respectively (Fig. 1a). When taking the expression of SLC14A1 in different cancer types and their related paracencerous tissue into consideration (Fig. 1b). Among many cancer tissues and normal tissue samples, SCL14A1 expressed differently and was often lower in common tumors. In addition, we found that the expression of SLC14A was significantly reduced in PCa tissue via IHC staining in PCa tissue chip. The positive composite score (Assignment of positive IHC results score) in the prostate cancer tissue was significantly reduced compared to the paracancerous tissue, p < 0.0001 (Fig. 1c,d).
The clinical data of 429 patients were obtained from TCGA. According to the cut-off value of SLC14A1 gene expression, which was calculated by "surv_cutpoint" of survminer R package, the patients were divided into 2 groups. Additional 106 patients whose data were downloaded from GEO were enrolled in our study for external verification. Kaplan-Meier survival curves also showed that highly expressed SLC14A1 tends to have a lower recurrence rate of PCa (Fig. 1e). Internal verification showed that the high expression group got less incidence of BCR, which indicated that SLC14A1 was a protective factor for PCa recurrence (Fig. 1f). www.nature.com/scientificreports/ A total of 429 samples with BCR follow-up time in TCGA (experimental group) were enrolled in our study, and Table 1 shows the differences between the high SLC14A1 expression group and the low expression group. The Gleason score of the high group was less than the low group (0.94 vs. 1.04, p < 0.001), which indicated that high expression of SLC14A1 might be a benefit of PCa prognosis. In addition, the low expression group had a much higher rate among stage pathologic N1 PCa compared with the high expression group (11.3 vs. 23.8, p = 0.002).
There were 6 biological pathways with the top enrichment score being selected by GSEA of the high expression group and low expression group (Fig. 2a). They were aminoacyl-tRNA biosynthesis, cell cycle, Fanconi anemia pathway, oxidative phosphorylation, protein export, and the ribosome.
In addition, the mRNAs from the comparison of the high expression group and low expression group were significantly involved in cell-substrate adhesion, O-glycan processing, and maintenance of gastrointestinal epithelium in GO. And the most different pathways in KEGG were the PI3K-Akt signaling pathway, cytokine-cytokine receptor interaction, and hippo signaling pathway (Fig. 2b,c) [17][18][19] .
The same cut-off point was used to distinguish the high expression group and the low expression group. A total of 22 miRNAs expressions were statistically different in the high expression group and the low expression group (log 2 foldchange > 2 and p < 0.05). And the ± top 3 miRNAs which were mostly related to SLC14A1 were selected by Pearson correlation (Fig. 2d). The 6 selected miRNAs were carried out to perform BCR-free Kaplan-Meier plots. Three miRNAs showed statistical differences comparing the high expression group and low expression group (Fig. 2e).
To identify the relationship between SLC14A1 and immune cells, differential analysis was performed using the Wilcoxon test. Violin illustration showed that B cells, cytotoxic cells, mast cells, neutrophils, NK cells, Tcm, Tem, Th1 cells, and T Reg were statistically different (Fig. 3a,b). Nevertheless, when taking them into univariate cox analysis, only B cells showed correlation with BCR events (Fig. 3c,d).
The correlation between SLC14A1 and B cells and Kaplan-Meier survival curves of B cells were shown in Fig. 3e-h. The Pearson correlation coefficient of the B cells and SLC14A1 were 0.41 and 0.43, respectively. Although the Pearson correlation coefficient only indicated that they were moderately related, Kaplan-Meier survival curves showed the significant statistic difference between B cells high group and B cells low group. To further validate independent risk factors for BCR events, multivariate cox regression analysis was performed ( Table 2). Low SLC14A1 expression was more likely to cause BCR events, and B cells, pathologic_N, and gleason score were all independent risk factors for BCR events.

Discussion
The incidence of PCa in men is very high, and RP is an important treatment modality. However, BCR after RP needs to arouse people's attention, because it will reduce the survival time of PCa patients after surgery. It is important to figure out the occurrence of BCR, and early detection of BCR will promote the treatment efficacy, which might include targeted radiotherapy or surgery 20 . Our research found that SLC14A1 was associated with the BCR of PCa, although some genes, such as cripto-1 (CR-1), abnormal spindle microtubule assembly (ASPM), C-X-C motif chemokine ligand 12 (CXCL12), epithelial cell transforming sequence 2 (Ect2), a four-long www.nature.com/scientificreports/ analysis. Each small point in the graph corresponds to a pathway, and the color is sorted by p value from small to large according to red, orange, yellow, green, blue, indigo, and purple. The smaller the p value is, the more color tends to red, and the larger the point is, the more genes in the pathway are. (d) Correlation analysis between ± top 3 differential expression miRNAs and SLC14A1. (e) Kaplan-Meier survival curves of the ± top 3 differential expression miRNAs that had a statistical significance. www.nature.com/scientificreports/ non-coding RNA (lncRNA) signature (RP11-108P20.4, RP11-757G1.6, RP11-347I19.8, and LINC01123) and pleomorphic adenoma gene like-2 (PLAGL2), were indicated to be related to PCa prognosis [21][22][23][24][25][26] , the studies of SLC14A1 in PCa were limited. SLC14A1 encodes type-B urea transporter (UT-B), which facilitates the rapid and passive cross-membrane movement of urea 27 . The low expression could induce urea accumulation which might influence the PCa cells. In the UT-B knock-out mice, the urea concentration reached about 9 times that of the wild type which caused severe apoptosis and DNA damage in the urothelial cells 28 . A high concentration of urea could cause cell cycle www.nature.com/scientificreports/ delay in the G2/M phase and G0/G1 phase, and lead to apoptosis and death of cells 29 . In addition, urea accumulation will alter arginine metabolism, thereby increasing the level of inducible NO synthase (iNOS) in cells 28 . Hypoxia-inducible factor-1 (HIF-1) is stabilized by a high concentration of nitric oxide (NO) catalyzed by iNOS, and HIF-1 is recognized to be related to tumors 30 . Therefore, in our study, the low expression of SLC14A1 seemed to be associated with BCR in PCa, the molecular mechanism might be multiple pathways induced by a high concentration of urea. Some researchers found that when transfecting cancer cells with miR-508 mimics, the cancer cells significantly reduce their ability of cell proliferation, migration, and invasion. And when cancer cells were transfected with miR-508 inhibitor, the effect was opposite to the previous 31,32 . There is no clear report on the relationship between miR-514a2 and cancers. But as for miR-449a, has a recognized tumor suppressor effect, especially in PCa. Gupta et al. found that c-Myc, a key factor that promotes cell cycle regulation, will be downregulated with miR-449a upregulated 33 . Similarly, Noonan et al. found that miR-449a was a miRNA component in the Rb pathway, and its tumor suppressor-like effect partly depends on the Rb status in PCa cells 34 . So miR-449a might be a key factor in PCa cells cycle arrest and senescence, which corresponding to our GESA results. Also, loss of miR-449a would cause PrLZ overexpression and promoted PCa metastasis 35 . And Cumar et al. found that SIRT1, a multifaceted NAD + -dependent protein deacetylase, will upregulate with the loss of miR-449a and promoted the invasion of PCa cells 36 .
Studies have shown that enrichment-based pathway diagrams on androgen-regulated proteomics datasets revealed significant imbalances in aminoacyl-tRNA synthetase, indicating increased protein biosynthesis which was a sign of PCa progress 37 . Cell cycle, which is closely related to the proliferation of PCa cells. The loss of controlling of the cell cycle monitoring mechanism can also occur in any link of cell proliferation damage induction, DNA repair, cell death, etc., finally cause the cell's genome instability. In addition, the biallelic mutation of Fanconi Anemia (FA) will not only lead to the phenotype of FA but also cause the instability of the genome and additional mutations in somatic cells, which leads to the susceptibility of many different types of cancer. And FA mutated had been found in PCa 38 . And the other pathways which were enriched in the SLC14A1 high expression group, such as oxidative phosphorylation, and ribosome, were all important pathways in the development and recurrence of PCa [39][40][41] .
Gleason grading is currently a widely used method of histological grading of prostate adenocarcinoma. Because Gleason classification is well related to biological behavior and prognosis, it is gradually recognized and used more and more widely, becoming an important reference indicator for formulating PCa treatment plans. A smaller Gleason score often indicates a better prognosis. In our study, patients with high expression of SLC14A1 got a lower Gleason score, which promoted survival of PCa patients and can predict the incidence of BCR.
B cells mainly exist in the lymph fluid circulating in the lymphatic vessels and are an important cellular component of the body's immune response function. The release of CXCL 13 from tumor cells could lead to B cells infiltration, which was considered as an important role in PCa progression 42 . In the PCa mouse model, the presence of immunosuppressive B cell subsets is associated with the accelerated recurrence of castrationresistant PCa 43 . In the pathological section of PCa, a large amount of B cell infiltration is related to the failure of chemotherapy of PCa 43,44 . These findings suggested that SLC14A1 might promote PCa immunology and a SLC14A1-targeted therapy could be created to delay tumor progression using the potential interaction of SLC14A1 and B cells.

Conclusion
Here, we identified a novel gene which related to the BCR of PCa and high expression of SLC14A1 could reduce the occurrence of BCR. This effect of SLC14A1 may be related to the interaction with miRNAs (has-miR-508, has-mir-514a2, and has-mir-449a) and the infiltration of B cells. However, further basic research and clinical trials are needed to further prove our point and determine the molecular pathway that how SLC14A1 can reduce BCR in prostate cancer.

Data availability
All our data is pulled from public databases. The datasets GSE32448, GSE46602, GSE69223, GSE70768 and GSE54460 used and analyzed in the current study are all from the GEO database(https:// www. ncbi. nlm. nih. gov/