TCGA dataset for Prostate Adenocarcinoma (PRAD), including a total of 460 patients with prognostic information, was obtained from PCaDB (Li et al., 2021a; Li et al., 2021b) (http://bioinfo.jialab-ucr.org/PCaDB/), a comprehensive and interactive database for transcriptomes from PCa population cohorts. In addition, the complete clinical data of 460 patients were downloaded from TCGA (https://portal.gdc.cancer.gov/) for subsequent analyses. The GSE54460 dataset was downloaded from the GEO database (http://www.ncbi.nlm.nih.gov/geo/). R software version: 4.1.0, website tool “catRAPID omics v2.0,” and website tool “OmicShare tools” were used for all analyses in the study.

Thirteen m5C-related genes were collected from literature mining (Takai et al., 2014; Archer et al., 2016; Blanco et al., 2016; Yamashita et al., 2017; Li Y. et al., 2018; Cheng et al., 2018; García et al., 2018; Chen et al., 2019; Gao et al., 2019; Janin et al., 2019; Carella et al., 2020; He et al., 2020; Mei et al., 2020; Sato et al., 2020). Then, expression data of the 13 m5C-related genes and all lncRNAs were retrieved from the TCGA dataset. Pearson correlation analysis (Schober et al., 2018) was used to explore the correlation between lncRNAs and the 13 m5C-related genes; the criteria were |Pearson R| > 0.4 and p < 0.05, and eventually 678 lncRNAs were filtered out. Next, logistic regression analysis (Tolles and Meurer, 2016) and univariate Cox regression analysis (Cox, 1972) were implemented on these 678 lncRNAs for further screening, and eventually, 102 overlapped lncRNAs with p-value < 0.01 were obtained.

We applied the LASSO regression model (R package “glmnet”) to shrink the candidate genes and, meanwhile, to establish the predictive model. Ultimately, 17 lncRNAs and their respective coefficients were calculated, and the minimum criteria chose the penalty parameter (λ). The m5C-lnc score was calculated by the formula: m 5 C − lnc Score = ∑ i = 1 N ( C o e f i × E x p r e s s i o n   l e v e l   o f   l n c R N A i ) where N (17) denotes the number of lncRNAs included in the predictive model, Coef _i refers to a specific lncRNA's coefficient, and Expression level of lncRNA _i represents the relative expression level of a specific lncRNA. The TCGA PCa patients were separated into low- and high-score subgroups according to the median score (–1.713). The BCR-free time was compared between the two subgroups via Kaplan–Meier (KM) survival analysis (Kaplan and Meier, 1958) in the “survminer” package. The “survivalROC” package was employed to perform 1-, 3-, and 5-year ROC (receiver operating characteristic) curve analyses (Mandrekar, 2010) for assessing the predictive power of the prognostic signature, and we compared the AUC (area under the ROC curve) of m5C-lnc score with the AUCs of other clinicopathological features to determine its clinical value. A cohort from the GEO database (GSE54460) was included for validation of the predictive model. The m5C-lnc score for each patient in the GSE54460 cohort was also calculated by the same formula above. Likewise, patients were separated into two subgroups—low- and high-score subgroups—for KM analysis. Finally, the two subgroups were then compared to validate the predictive model, too. We constructed the nomogram and calibration plots based on our previous study (Zhong et al., 2021).

The “GDCRNAtools” package was employed to construct the ceRNA network (Salmena et al., 2011), and the website tool “catRAPID omics v2.0” (http://service.tartaglialab.com/page/catrapid_omics2_group) was utilized to construct the lncRNA–protein interaction network. The website tool “OmicShare tools” (https://www.omicshare.com/tools/) was then used to explore the functional enrichments of the prognostic lncRNAs.

PCa cell lines, DU145 and PC-3, were both obtained from the National Collection of Authenticated Cell Cultures. Both cell lines were cultured in Dulbecco's modified Eagle's medium. Culture media contained 10% fetal bovine serum and 1% double antibiotics (penicillin and streptomycin). All cell lines were cultivated at 37°C and 5% CO₂. Total RNAs were extracted from PC-3 and DU145 cells using Trizol reagent (15596018, Takara). Total RNAs were then reverse-transcribed into cDNA using TransScript All-in-one First-Strand cDNA Synthesis SuperMix for qPCR (AT341-01, TransGen). RT-qPCR (real-time quantitative PCR) was carried out using the PerfectStart Green (AQ601-02, TransGen) in Applied Biosystems 7500 Real-Time PCR System. Finally, the relative expression of MAFG-AS1 was calculated based on the internal reference glyceraldehyde 3-phosphate dehydrogenase (GAPDH). All experiments were carried out with three replicates. The primers of MAFG-AS1 and GAPDH are listed in Supplementary Table S3.

GenePharm Company synthesized small interfering RNAs (siRNAs) targeting MAFG-AS1. RT-qPCR validated the transfection efficiency after the transfection of siRNAs along with siRNA-Mate (GenePharm) for 72 h. CCK-8 (Cell Counting Kit-8, MA0218-5, Meilunbio) viability assay and colony formation assay were applied to examine the proliferative ability of PCa cell lines after knocking down MAFG-AS1. Transwell assay was conducted to examine the invasiveness of MAFG-AS1 downregulated cells. Detailed procedures of the above assays are available in our previous study (Zhong et al., 2021). All experiments were implemented with three triplicates. siRNAs targeting sites in MAFG-AS1 are shown in Supplementary Table S3. The original scans of the plate colony formation assays and Transwell assays are shown in Supplementary Figure S3.

All bioinformatics analyses were performed by R software version 4.1.0 (The R Project for Statistical Computing, Vienna, Austria). Pearson correlation analysis was carried out to analyze correlations between m5C-related genes and lncRNAs. The “survival” and “survminer” packages were applied for KM plots and Cox regression analysis. GraphPad Prism 7.0 (GraphPad, La Jolla, CA, United States) was utilized to analyze the results of RT-qPCR and cell functional assays. All statistical results were showed as mean ± SD (standard deviation) with two-sided test, and the results with a p-value < 0.05 were regarded as statistically significant ones.

The workflow chart is presented in Figure 1A. Firstly, we searched and identified 13 m5C- regulators from published articles, and then TCGA and GEO datasets of prostate adenocarcinoma were downloaded for further analysis. After extracting the expression matrixes of these m5C regulators and all lncRNAs in the TCGA dataset and excluding lncRNAs with FPKM < 1 (Fragments Per kilobase Million), we performed a Pearson correlation analysis between these m5C-related genes and 2,590 lncRNAs to determine whether a lncRNA was correlated with the m5C modification (|Pearson R| > 0.4 and p < 0.05), and consequently we included 678 qualified lncRNAs into logistic regression analysis and univariate Cox regression analysis for further screening (Supplementary Tables S1 and S2). Subsequently, 102 lncRNAs, whose p-values were both less than 0.01 in the two analyses mentioned above, were chosen for a LASSO regression analysis to build a m5C-related prognostic model to forecast the BCR in PCa (Supplementary Figure S1A). Eventually, we constructed a model that contains 17 m5C-related lncRNAs (Supplementary Figure S1B), and the correlations between the 17 lncRNAs and the 13 m5C regulators in the TCGA dataset are shown in Figure 1B.

Given the LASSO regression analysis, we focused on these 17 lncRNAs' clinical value. The univariate Cox analysis results of the 17 lncRNAs are shown in Figure 2A; 7 lncRNAs (UBXN10-AS1, AC017100.1, Z93930.2, LINC01135, AP001486.2, AC004943.2, and SP2-AS1) tend to be protective factors with HR (hazard ratio) < 1, while another 10 lncRNAs' HRs are more than 1 indicating their high-risk relation with the BCR of PCa. We also separated patients into high-expression and low-expression subgroups based on the median of each lncRNA's expression value. We performed KM survival analysis to detect a significant difference in BCR status between two subgroups. As a result, all 17 KM survival analyses revealed that the high-expression subgroups belonging to the 7 protective lncRNAs mentioned above harbored favorable BCR status. In contrast, the high-expression subgroups of another 10 high-risk lncRNAs gained poorer BCR status, and meanwhile, all the KM analyses hit a p-value < 0.05, as shown in Supplementary Figure S2. Then we extracted 460 patients' clinical data, including BCR status, Gleason score (GS), PSA, pathological T stage, pathological N stage, and metastatic stage from the TCGA-PRAD dataset to inspect the 17 lncRNAs' signatures, and the heatmap is displayed in Figure 2B. Moreover, we implemented a further analysis to explore the correlation between 17 m5C-related modulators and clinical features after classifying these patients into several subgroups, as shown in Figure 2C and Supplementary Figure S1E. The expression of all lncRNAs except AC004943.2 in tumors exhibited significance, compared with that in adjacent normal tissue. Regarding the T stage, N stage, and GS subgroups, five, four, and two lncRNAs did not show significant difference, respectively. Consequently, most of these lncRNAs revealed their predictive value.

Following the analyses above, we then constructed the m5C-related lncRNA prognostic model named m5C-lnc score with these 17 lncRNAs, and their coefficients calculated by LASSO regression analysis are shown in Figure 3A. Patients in the TCGA cohort were divided into low-score and high-score subgroups based on the m5C-lnc score attained by the formula above and its median. KM survival curves demonstrated that patients in the low-score subgroup had better clinical outcomes (longer BCR-free time), as shown in Figure 3B. m5C-lnc score and BCR status distribution of TCGA-PRAD dataset are depicted in Figure 3C. Subsequently, ROC analysis showed that the predictive m5C-lnc score model whose 1-, 3-, and 5-year AUCs were all about 0.8 (12-month AUC = 0.785, 36-month AUC = 0.825, and 60-month AUC = 0.815; Figure 3D) had excellent sensitivity and specificity, exhibiting better power to predict BCR status in the TCGA cohort; moreover, compared with the AUCs of other prognostic factors such as PSA, age, GS, and pathological T stage (AUC of PSA = 0.612, AUC of age = 0.51, AUC of GS = 0.728, and AUC of pathological T stage = 0.685), the AUC of the m5C-lnc score hit the highest score (0.793; Figure 3E), suggesting the model's better sensitivity, specificity, and clinical predictive value. The decision curve analysis (DCA) of the TCGA cohort indicated that the clinical benefit rate of the m5C-lnc score model was significantly higher than that of other clinical characteristics based on GS and T stage, respectively (Supplementary Figure S1C). Moreover, we performed univariate Cox regression with m5C-lnc score and other clinical pathological features including age, PSA, GS, and pathological T stage, respectively. As a result, all factors except age were significantly associated with patients' BCR status, as shown in Supplementary Figure S4A; PSA, GS, m5C-lnc score, and pathological T stage were then included in further multivariate Cox regression, and consequently, m5C-lnc score presented as the primary risk factor with the highest HR of about 3, as shown in Supplementary Figure S4B.

To validate the predictive value of the established m5C-related lncRNA signature, we applied the GSE54460 dataset to verify the results; likewise, patients in the GSE54460 dataset were also separated into low-score and high-score groups by the median m5C-lnc score. Consequently, the survival analysis results were similar to those in the TCGA dataset: patients with higher m5C-lnc scores had a lower BCR-free rate and a shorter BCR-free time (Figure 3F). Likewise, m5C-lnc score and BCR status distribution are illustrated in Figure 3G. The ROC analysis indicated that these m5C-related regulators' signature also had a promising prognostic value for patients in the GSE54460 dataset given that the 1-, 3-, and 5-year AUCs were all around 0.7 (12-month AUC = 0.748, 36-month AUC = 0.68, and 60-month AUC = 0.655; Figure 3H). Nomogram is always used to quantitatively predict the prognosis of patients in clinical practice by calculating the relative points on a set of scales. In our study, we illustrated our predictive model in the form of a nomogram as shown in Figure 3I to predict 1-, 3-, and 5-year BCR-free probabilities, respectively, for patients with PCa according to the “Total Points.” A total point considers six variables including PSA, age, T-stage, GS, m5C-lnc score, and resection status. With the scale bar “Points” to be a reference, the primary value of each variable could correspond to a point on the “Points” scale, and all the six relative points are added together to make the total points. Ultimately, the probabilities of 1-, 3-, or 5-year BCR-free time for each patient could be quantified respectively by the points corresponding to the scale bar “Total Points.” Besides, we also carried out a calibration curve analysis to demonstrate that the estimated 3- and 5-year BCR-free time were in accordance with the actual scenario (Supplementary Figure S1D).

Next, we explored how these lncRNAs function in biological processes by using the R package called “GDCRNAtools” to construct a ceRNA network and visualizing it with the help of the software Cytoscape (Shannon et al., 2003). Consequently, 9 lncRNAs (UBXN10-AS1, AC138956.2, ZNF32-AS2, AC017100.1, AC004943.2, SP2-AS1, Z93930.2, AP001486.2, and LINC01135), 106 miRNAs, and 1,710 mRNAs were included in the ceRNA network with their lncRNA–miRNA–mRNA pathways interacting mutually, as shown in Figure 4A. For the rest of lncRNAs (AC012510.1, AC012065.3, AL117332.1, AC132192.2, AP001160.2, AC129510.1, AC084018.2, and MAFG-AS1), we utilized an online website tool called catRAPID omics v2.0 (Agostini et al., 2013) to inspect the lncRNA–protein interaction. After inputting each lncRNA's information, the tool automatically predicted the candidate proteins and presented the result list. Given the build-in ranking system (0 star to 3 stars), we selected 27 candidate proteins with 2.5 and more stars to construct the interaction network. Likewise, the lncRNA–protein interaction network was visualized by Cytoscape, as shown in Figure 4B. Detailed data can be found in Supplementary Tables S4 and S5.

Furthermore, we used these 17 lncRNAs to perform functional enrichment analysis with an online tool called OmicShare tools. The summary of GO (Gene Ontology) enrichment and details about its three parts are shown in Figures 5A–D, and the KEGG (Kyoto Encyclopedia of Genes and Genomes) pathway enrichment is also demonstrated in Figure 5E. Accordingly, we found that the top five GO Biological Processes that these genes were primarily enriched in are heterocycle metabolic process, nucleobase-containing compound metabolic process, organic cyclic compound metabolic process, cellular nitrogen compound metabolic process, and cellular aromatic compound metabolic process (Figure 5B); neomycin, kanamycin, and gentamicin biosynthesis, carbon metabolism, and RNA transport were significantly enriched in KEGG Pathway (Figure 5E). These data may bring us some insights into the potential functions of these m5C-related lncRNAs in PCa. Complete information are shown in Supplementary Table S6.

With the univariate Cox regression analysis mentioned above, 10 lncRNAs' HRs are greater than 1, indicating their potential promoting PCa. Among them, the abundance of MAFG-AS1 is the highest (Figure 2C), and its HR ranks second (Figure 2A). Therefore, we chose to implement the experimental validation with MAFG-AS1. Based on the data provided on LNCipedia (Volders et al., 2019), a website with comprehensive information of human lncRNAs, MAFG-AS1, known as MAFG-DT, is located at q25.3 of chromosome 17 with a transcript size of 1,895 bp. Subsequently, we used an online tool called iRNAm5C-PseDNC to identify 5-methylcytosine modified sites of MAFG-AS1, and the predicted sites included positions 4, 6, 10, 11, 13, and 17, as shown in Figure 6A.

Firstly, we detected the basal expression level of MAFG-AS1 in PCa cell lines by RT-qPCR, and it was found that MAFG-AS1 was highly expressed in all PCa cell lines with the benign prostate cell line RWPE1 to be the reference (Figure 6B). The expressions of MAFG-AS1 in PC-3 and DU145 cell lines are much higher than that of others, so we chose these two cell lines to construct the MAFG-AS1 knock-down phenotypes with two siRNAs, respectively. The results of RT-qPCR confirmed that both siRNA fragments significantly silenced the expression of MAFG-AS1 in PC-3 cell line and DU145 cell line (Figure 6C). After examining the transfection efficiency in PC-3 and DU145 cell lines, we performed a set of assays to verify this specific gene's function on cancer cells. CCK-8 assay and plate colony formation assay suggested that down-regulated MAFG-AS1 significantly inhibited PCa proliferation ability (Figures 6D, E). Furthermore, knockdown of MAFG-AS1 attenuated PCa cell's invasion ability via transwell assay (Figure 6F). These results indicated that MAFG-AS1 acts as a high-risk predictive factor in PCa, and somehow its high expression promotes cancer progression.

After determining MAFG-AS1's promotion in PCa's proliferation and invasion, we attempted to inspect its molecular function further. First of all, we checked out MAFG-AS1's expression profile with a pan-cancer analysis on GEPIA 2 (Tang et al., 2017) as shown in Figure 7A, and it turned out that MAFG-AS1 is highly expressed in most cancers. Additionally, we detected a significant difference in the expression of MAFG-AS1 in tumor tissue and adjacent normal tissue in another two datasets, CIT dataset and Cambridge Dataset (GSE70768), as shown in Figures 7B, C; the expression of MAFG-AS1 is higher in tumor tissue than in the adjacent normal tissue. Meanwhile, we also explored the relation between MAFG-AS1's expression and clinical characteristics in patients with PCa. It was found that patients with high GS tended to have high MAFG-AS1 expression level (GS = 6: 1.29 ± 0.70, GS = 7: 1.58 ± 0.84, and GS ≥ 8: 1.84 ± 0.92, p < 0.001), and the significant difference was also observed in the pathological T stage category, lymph node metastasis category, and BCR status category (≤T2c: 1.52 ± 0.81, T3/T4: 1.75 ± 0.92, p = 0.005; lymph node metastasis-positive: 1.89 ± 0.99, lymph node metastasis-negative: 1.62 ± 0.87, p = 0.018; BCR-positive: 1.94 ± 0.91, BCR-negative: 1.58 ± 0.86, p < 0.001), as shown in Supplementary Table S8. Although it was not included in the ceRNA network (Figure 4A), it interacts with 15 proteins in the lncRNA–protein network, ranking first, which indicates it might play an essential role in cancer progression. Then we implemented Spearman's correlation analysis between MAFG-AS1 and all other genes in TCGA-PRAD dataset to screen out the correlated genes of MAFG-AS1 for further functional enrichment analysis. The criteria were |Spearman R| > 0.4 as well as p-value < 0.05, and 833 correlated genes were found. Next, we utilized these genes to perform functional enrichment analysis via the website tool OmicSharetools. Consequently, most genes were enriched in these terms: RNA processing, mRNA metabolic process, ribonucleoprotein complex biogenesis, protein targeting, and mRNA catabolic process (GO Biological Process) and ribosome, spliceosome, and RNA transport (KEGG Pathway) (Supplementary Figures S5A–E). Complete functional enrichment documents can be found in Supplementary Table S7. These results suggested that MAFG-AS1 might influence cancer development by interacting with proteins to interrupt the activities of RNAs, especially mRNAs. Lastly, we used LNCipedia to predict MAFG-AS1's protein-coding potential. The website would carry out its prediction with various applications such as CPC, HMMER, PRIDE, PhyloCSF, CPAT, and Ribosome-profiling. As a result, only CPAT provided a coding probability of 76.79% (Figure 7D). It seems that MAFG-AS1 encodes protein with little potential.