Prognostic Risk Model and Tumor Immune Environment Modulation of m5C-Related LncRNAs in Pancreatic Ductal Adenocarcinoma

RNA methylation modification is a key process in epigenetics that regulates posttranscriptional gene expression. With advances in next-generation sequencing technology, 5-methylcytosine (m5C) modification has also been found in multiple RNAs. Long non-coding RNAs (lncRNAs) were proved to have a key role in cancer progression and closely related to the tumor immune microenvironment. Thus, based on the PDAC patients’ clinical information and genetic transcriptome data from the TCGA database, we performed a detailed bioinformatic analysis to establish a m5C-related lncRNA prognostic risk model for PDAC patients and discovered the relationship between the risk model and PDAC immune microenvironment. Pearson correlation coefficient analysis was applied to conduct a m5C regulatory gene and m5C-related lncRNA co-expression network. Expression of m5C-related lncRNAs screened by univariate regression analysis with prognostic value showed a significant difference between pancreatic cancer and normal tissues. The least absolute shrinkage and selection operator (LASSO) Cox regression method was applied to determine an 8-m5C-related lncRNA prognostic risk model. We used principal component analysis to indicate that the risk model could distinguish all the samples clearly. The clinical nomogram also accurately predicted 1-, 1.5-, 2-, and 3-year survival time among PDAC patients. Additionally, this risk model was validated in the entire group and sub-test groups using KM analysis and ROC analysis. Combined with the clinical characteristics, the risk score was found to be an independent factor for predicting the survival of PDAC patients. Furthermore, the association between the risk model and tumor immune microenvironment was evaluated via the ESTIMATE R package and CIBERSORT method. Consequently, the results indicated that immune cells were associated with m5C-related lncRNA risk model scores and had different distribution in the high- and low-risk groups. Based on all these analyses, the m5C-related lncRNA risk model could be a reliable prognostic tool and therapeutic target for PDAC patients.


INTRODUCTION
Pancreatic ductal adenocarcinoma (PDAC) is a common malignant tumor of the digestive system, which is characterized by high degree of malignancy, easy recurrence, and metastasis (1). Recently, the incidence rate of PDAC has increased year by year. Due to the lack of early diagnosis methods, most patients were diagnosed at an advanced stage and lost the chance of radical resection (2). Lately, personalized molecular therapies targeting specific tumor-associated genes and the application of immune checkpoint inhibitors have effectually meliorated the prognostic effect of patients with advanced cancer (3). However, pancreatic cancer is characterized by high rates of drug resistance and metastasis due to its specific tumor microenvironment, current therapies have not effectively contributed to the prognosis of patients with PDAC, and researchers still need more specific biomarkers for developing more effective diagnostic and treatment strategies (4). Hence, it is essential to find new molecular biomarkers for early diagnosis of PDAC.
Epigenetics study aims to research gene differential expression without changing nucleotide sequence (5), and epigenetic modification has been proven to play an important role in the progress of tumor biology (6), such as DNA and RNA methylation, genomic imprinting, gene silencing, and non-coding RNA modification. In recent years, N6-methyladenine (m6A) and 5methylcytosine (m5C) RNA methylation have become a vital direction in detecting the regulation of RNA epigenetic modification upon the growth of various malignant tumors (7,8). With the advance of high-throughput sequencing technology, the distribution and biological functions of RNA m5C modification at mRNA and non-coding RNA have been gradually discovered (9). RNA m5C methylation is catalyzed by a methyltransferase complex, which consists of three methylation-related enzymes, including methyltransferase ("Writer"), demethylase ("Erase"), and m5C binding protein ("Reader"). Similar to functions of m6A methylation, m5C methylation exerted biological effects mainly by regulating RNA stability, affecting transcription efficiency, and mediating RNA localization (10). In addition, abnormal expression of m5C methylation is associated with the occurrence and development of several malignant tumors (11,12). Recent reports have also demonstrated that m5C-related gene expression is associated with prognosis in patients with lung and pancreatic cancers, which further indicated that m5C methylation may have an essential role in malignant tumor growth (13,14).
Long non-coding RNAs are transcribed by corresponding genes and have a similar structure to mRNA including polyA tail and promoter structure (15). Due to different splicing methods, multiple lncRNAs were formed during the differentiation process. Several research unveiled that lncRNAs could modulate downstream genes in epigenetic, transcriptional, and post-transcriptional levels including gene silencing, histone processing, transcriptional regulation, transcriptional interference, and nuclear transport, which were closely related to the development of various human diseases (16). In malignant tumors, methylation-related genes could affect tumor progression by regulating the methylation level of lncRNA. For instance, LINC00942-specific sequence recruited methylated transferase METTL14 and stabilized the downstream targets of LNC942, which could promote tumorigenesis and development in breast cancer (17). Several reports have demonstrated that m5C modification existed widely in non-coding RNAs; nevertheless, the reports about m5C methylation regulation in lncRNAs were still few. Thus, it is essential to further research the relationship between m5C methylation and lncRNAs, especially in malignant tumor.
Tumor microenvironment (TME) is a mixture of fluid, immune cells, and blood vessels surrounding the tumor (18). The interaction between tumor cells and TME could determine the progression and fate of tumor. Therefore, understanding the composition and function of TME is crucial to research cancer progression (19). Studies have indicated that multiple genetic mutations could increase the incidence of cancer. However, the impact of TME on cancer development is still unclear. Accumulating evidence revealed that lncRNAs were important regulators in the immune system, which could control the distribution and activation of immune cells and mediate the tumor cells' immune evasion. For instance, LncRNA-MALAT1 could promote thyroid cancer progression via modulating tumor-associated macrophages secreting FGF2 protein (20). In colorectal cancer, lncRNA SATB2-AS1 could regulate tumor tissue immune cell density via combining with downstream proteins WDR5 and GADD45A and regulating TH1-type chemokines (21). Furthermore, the relationship between lncRNAs and tumor immune microenvironment has also been reported in pancreatic cancer (22). However, the relationship between m5C-related lncRNAs and tumor-infiltrating lymph cell distribution has never been reported in PDAC.
In our study, we summarized 177 PDAC patients' clinical information together with transcriptome expression data, and constructed an mRNA-lncRNA co-expression network between m5C-related mRNAs and lncRNAs. A total of 17 prognostic lncRNAs were screened via univariate Cox analysis. Consequently, we generated a PDAC prognostic risk model consisting of 8 m5C-related lncRNAs via LASSO regression analysis. We further verified the prognostic risk model in subgroups via Kaplan-Meier (KM) analysis and receiver operating characteristic curve (ROC) analysis. Additionally, the association between the risk model and tumor immune environment was measured via the ESTIMATE R package and CIBERSORT tool. Finally, 6 types of immune cells were identified in high-or low-risk groups. Moreover, m5C-related lncRNAs had statistically significant association with immune related genes and 5 types of tumor-infiltrating lymph cells were negatively correlated with the risk scores. In summary, our research indicated that the m5C-related lncRNA risk model might offer a promising prognostic tool and play an essential role in regulating PDAC immune cell distribution.

Data Collection and Processing
We downloaded a total of 182 patient data from the TCGA data website including clinical and transcriptome expression raw data

LASSO Analysis
The LASSO regression method was introduced to construct a m5C regulatory gene-related lncRNA predicted risk model via R package. It is a widely used high-dimensional index regression method that screened the m5C-related lnRNAs with prognostic value and constructed a predicted risk model by applying a penalty proportional to the contraction of the regression coefficient. We established a m5C-related lncRNA risk model consisting of 8 lncRNAs depending on this method.

GSEA
We divided all samples into two groups according to the median risk score of the m5C-related lncRNAs. GSEA software (GSEA_4.1.0) was performed to analyze the data. Enrichments of gene sets with a normalized p-value less than 0.05 were considered to be significant.

Identification Tumor-Infiltrating Immune Cells
CIBERSORT is a typical tool to deconvolute the expression matrix of immune cells based on linear support vector regression, which quantified the infiltrating immune cells proportion via detecting marker gene expression. In our study, we combined the expression of marker genes from 22 types of immune cells and all the PDAC patients' transcriptome data to acquire the tumor-infiltrating lymph cell distribution scores via CIBERSORT. ESTIMATE R package could use gene expression profile to predict stromal cell and immune cell scores, and then calculate the amount of these two types of cells, which was performed to analyze the pancreatic cancer purity in this research.

Statistical Analysis
The data flow chart is shown in Supplementary Table 2. The m5C-related lncRNAs with correlation index > 0.3 and p < 0.05 were regarded as statistically significant. Univariate and multivariate Cox regression analyses were performed to identify the m5C-related lncRNA risk score as an independent prognostic factor for PDAC. The Kaplan-Meier method was performed to compare OS time of PDAC patients. ROC analysis was expended to measure the prognostic competence between the prognostic risk scores and other clinical parameters. The qRT-PCR experiments were analyzed by PRISM 7. The relationship between the m5C-related lncRNA risk score and tumor-infiltrating lymph cell distribution was analyzed by CIBERSORT and ESTIMATE R package.

Identification of m5C Regulator-Related LncRNAs in PDAC Patients
Initially, we screened the m5C methylation modification-related genes based on the published articles. A total of 13 m5C regulators were selected, namely, YBX1, ALYREF, DNMT1, NSUN4, TRDMT1, TET2, NSUN7, NSUN6, NSUN5, NSUN3, NSUN2, DNMT3a, and DNMT3b. We summarized the m5C regulator genes and then compared the expression of m5C regulator genes in 179 normal pancreas tissues and 171 pancreatic cancer tissues using online data base GEPIA (http://gepia.cancer-pku.cn/). We found that among the m5C regulators, the expression of YBX1, ALYREF, DNMT1, and NSUN4 was obviously higher in PDAC patient's tissues (p < 0.05), whereas no notable variances were detected in TRDMT1, TET2, NSUN7, NSUN6, NSUN5, NSUN3, NSUN2, DNMT3a, and DNMT3b (Supplementary Figure 1). Moreover, the m5C-related lncRNAs were recognized by correlation analysis depending on the expression of m5C regulators and lncRNAs in a total of 177 PDAC samples. Use correlation coefficient > 0.3 and p-value > 0.05 as filter criterion; finally, 242 m5C-related lncRNAs were screened. We also constructed a network between m5C-related genes and their coexpressing lncRNAs to show m5C-related lncRNA co-expressing relationship ( Figure 1A). Based on the screened m5C-related lncRNAs, we use univariate Cox analysis (p < 0.001) to further filter the prognostic m5C-related lncRNA in PDAC patients combined with the patient survival data. The results indicated that a total of 17 PDAC prognostic lncRNAs were screened, and most of the m5C-related lncRNAs were protective factors (hazard ratio, HR < 1), only CASC8 was a risk factor (HR > 1) in PDAC. In addition, the Circos picture showed that the protective m5C-related lncRNAs had positive expression correlation and risk factor CASC8 was negatively correlated with others ( Figures 1B, C). Furthermore, we tried to clarify the expression of the screened prognostic m5Crelated lncRNAs in PDAC patients depending on the TCGA database. The heatmap revealed that all the lncRNAs had statistical differences between the normal and tumor pancreatic tissues, which indicated that the m5C-related lncRNAs might play a key role in PDAC progression ( Figures 1D, E).

Establishing an 8 m5C-Related LncRNA Risk Model for PDAC Patients
Depending on the primary screen of the m5C-related lncRNAs, the LASSO regression technique was further performed randomly to construct a prognostic risk model for PDAC, which showed that eight m5C-related lncRNAs were appropriate for building the prognostic risk model (Figures 2A, B and Table 1). The risk score was calculated as follows: risk score = (-0.780839063578865 * AC022098.1. expression) + (-0.220638925265728 * AL031775.  Figure 2C). Furthermore, all the PDAC patients were divided into low-or high-risk groups based on the median risk scores calculated by the above formula. The principal component analysis (PCA) and three-dimensional PCA showed that patients with different risks were well separated into two clusters ( Figures 2D, E). Patients in the high-risk group had significantly shorter survival time than the patients in the lowrisk group. A distinguished difference in OS time was detected between the low-and high-risk groups via Kaplan-Meier survival analysis (HR = 2.47, 95% CI: 1.64-3.73, p < 0.001, Figure 2F).

Relationships Between m5C-Related LncRNAs and Clinical Pathological Parameters
We further detected the impact of the selected eight lncRNAs on PDAC patients' overall survival time, respectively. The OS curve of the m5C-related lncRNAs showed that patients with risk lncRNA high expression (CASC8) had shorter survival time, whereas those with protective lncRNA high expression (AC022098.1, AL031775.1, AC005332.6, AC096733.3, AC025165.1, AC009974.1, and PAN3-AS1) had longer survival time ( Figure 3A). According to the heatmap, the pancreatic  tumor size was significantly different between the high-and lowrisk groups (p < 0.05). Nevertheless, the other clinical factors including patient age, gender, tumor stage, and grade had no statistical differences ( Figure 3B). Moreover, we further subdivide these clinical indicators severally and analyzed the risk scores in the subgroups. The KM survival curve showed that the patients with high-risk scores had a shorter OS in male patients' group, younger patients (age less than 65 years), Stage I-II groups, Grade I-II or Grade III-IV groups, T I-II or T III-IV groups, N0 or N I-III groups, and patients without any metastasis group ( Figure 3C and Supplementary Figure 2).
The m5C-Related LncRNA Risk Score Was an Independent Prognostic Factor for PDAC We ordered the patients based on the risk scores. Heatmap revealed that the expression of 8 m5C-related lncRNAs were considerably different in the PDAC patients with different risk scores. Furthermore, the scatter dot plot also unveiled that the mortality of the PDAC patients increased with the risk scores, and patients with lower risk scores exhibited longer survival time ( Figure 4A). Next, we performed univariate and multivariable Cox regression analyses to validate if the risk score calculated by the m5C-related lncRNA risk model could work as an independent prognostic factor for PDAC patients. As shown in Figure 4B, univariate Cox regression analysis demonstrated that only the m5C-related lncRNA risk scores were obviously positively related with OS (HR: 197.088, 95% CI: 10.283-3777.290, p < 0.001). Multivariate analysis also revealed that m5C-related lncRNAs prognostic risk score (HR: 116.786, 95% CI: 3.668-3718.415, p = 0.007) had a significant relationship with PDAC patients' OS and could act as an independent prognostic factor ( Figure 4C). Our results suggested that the m5C-related lncRNA risk model was an independent prognostic factor for PDAC and took advantage of the traditional clinicopathological indicators including patient age, gender, and tumor stage. Moreover, the 1-year ROC curve proved that the AUC value for the m5C-related lncRNA risk model was 0.716, which was    Figure 4D). In addition, the 3-year ROC curve analysis also suggested that the risk score AUC value (AUC = 0.639) was a superior prognostic factor for PDAC patients ( Figure 4E).

Exanimating the Prediction Value of the m5C-Related LncRNA Risk Model
In order to predict the PDAC patient overall survival accurately, a nomogram was constructed to reveal the 1-, 3-, and 5-year survival rates based on the m5C-related lncRNAs expression risk scores together with clinicopathologic features (including age, gender, grade, T stage, M stage, and N stage) ( Figure 5A). In addition, we used a calibration curve to compare the consistency of the actual and the predicted 1-, 1.5-, 2-, and 3-year patient survival. We found that the actual and the predicted line were almost in accordance within 3 years ( Figure 5B). The above results proved that the nomogram we generated via m5C-related lncRNA prognostic risk scores was dependable. In addition, we randomly divided all the PDAC patients from TCGA into two subgroups (group A and B) in a 1:1 ratio and performed an internal validation for the m5Crelated lncRNA risk model. The KM survival curve and 5-year ROC curve were examined in each subgroup. Our results showed that the patients with higher m5C-related lncRNA risk scores had shorter OS in group A (HR: 2.05, 95% CI: 1.13-3.71, p = 0.019) and the AUC value of the 5-year ROC curve was 0.814 ( Figures 5C, D). The patients in group B had a similar OS trend (HR: 2.56, 95% CI: 1.44-4.56, p = 0.001) and the AUC value was 0.903 ( Figures 5E, F). Moreover, the scatter dot plot also revealed an obvious difference expression of the screened lncRNAs in high-or low-risk groups, and the patient survival time was positively related with the risk scores in all the subgroups (Supplementary Figure 3). The above results demonstrated that the m5C-related lncRNA risk model was a reliable predictive factor for the PDAC patient.

Detecting LnRNA Expression In Vitro and Functional Enrichment Analysis
The m5C-related lncRNA expression profiles from the TCGA database indicated that AC022098.1, AL031775.1, AC005332.6, AC096733.3, AC025165.1, AC009974.1, and PAN3-AS1 were downregulated, and the CASC8 was overexpressed in PDAC. In order to examine the lncRNAs' expression level in vitro, three pancreatic cancer cell lines and a normal pancreatic duct cell line were used to perform the qRT-PCR experiments. The in vitro results were not completely consistent with the data from TCGA ( Figure 6A). We found that CASC8 was upregulated in pancreatic  cancer cells (Mia-PaCa-2, CFPAC-1, and Panc-1 cells) and AC096733.3 was in low expression compared with HPNE cells, which were with the same expression profile with TCGA data. AC096733.3, AC025165.1, PAN3-AS1, and AC009974.1 were downregulated in at least two cancer cell lines and may play as tumor suppressor genes in PDAC. The above results were partly consistent with the patients' OS analysis data from TCGA. However, AC022098.1 and AC005332.6 were suppressed in Mia-PaCa-2 and overexpressed in CFPAC-1, and their expression was inconsistent in different cells. In addition, AL031775.1 was only downregulated in Panc-1 cells. Hence, the underlying mechanisms that m5C-related lncRNAs regulating PDAC patient OS time still needs to be further explored. To seek the signal pathways that m5C-related lncRNAs might involve in the low-or high-risk groups, GSEA was performed. We found that the MTORC1 signaling pathway was activated in the highrisk group; nevertheless, MYOGENESIS and KRAS signaling pathways were activated in the low-risk group (Figures 6B, C). Furthermore, we summarized differential genes in the low-or high-risk groups to conduct GO and KEGG enrichment analysis. The top five GO enrichment biological processes were T-cell activation, Calcium ion homeostasis, Hormone transport, Modulation of chemical synaptic transmission, and Regulation of trans-synaptic signaling pathway (Figures 6D-F). Additionally, KEGG analysis revealed that several immune-related pathways were enriched such as Primary immunodeficiency, T-cell receptor signaling pathway, Leukocyte transendothelial migration, and Th1/Th2 cell differentiation ( Figures 6G, H).

The Relationship Between m5C-Related LncRNAs and Tumor-Infiltrating Lymph Cells
We analyzed the differences of tumor microenvironment in the high-or low-risk groups. By using CIBERSORT tool with values of p < 0.05, a total of 22 tumor-infiltrating immune cells were screened ( Figure 7A). The results demonstrated that naïve B cells, CD8 + T cells, regulatory T (Treg) cells, and resting NK cells exhibited a higher expression in the low-risk group (p < 0.05), whereas the M0 and M2 phenotype macrophages had a higher expression in the high-risk group. M2 phenotype macrophage has been proven as a co-carcinogenic factor in pancreatic cancer. Hence, m5C-related lncRNAs may promote M2 phenotype macrophage polarization or infiltration in PDAC. Additionally, we further detected the relationship between m5C-related lncRNAs and immune-related genes (LMTK3, LAG3, CD27, CD28, CD86, and BTLA), which indicated that most of the lncRNAs had statistical correlation with immune-related genes except AC025165.1 (p < 0.05) ( Figure 7B). Gene expression analysis also showed that expression of the above immunerelated genes are different in low-or high-risk groups ( Figure 7C). Moreover, tumor microenvironment scores were performed to analyze the stromal and immune cell proportion in  the tumor environment via ESTIMATE R package. We obviously found that the PDAC tissues with low m5C-related lncRNA risk scores had higher immune scores, stromal scores, and overall ESTIMATE scores, which indicated that the PDAC in the lowrisk group had lower pancreatic tumor proportion ( Figure 7D). Lastly, we also investigated the association between the m5Crelated lncRNA risk score and tumor lymph cells. Spearman's correlation analysis revealed that the risk score was positively related with 3 tumor-infiltrating lymph cells (NK cells, M0, and M2 phenotype macrophages), but was negatively correlated with regulatory T cells, CD8 + T cells, activated memory CD4 + T cells, naïve B cells, and plasma cells (p < 0.05) ( Figure 7E). These outcomes revealed that the m5C-related risk scores could discriminate the different features of tumor-infiltrating lymph cells in PDAC.

DISCUSSION
Pancreatic cancer is an extremely malignant cancer, and its incidence is close to mortality rate. In recent years, although comprehensive treatments based on surgery have made some progress in the treatment of PDAC patients, the diagnosis and prognosis of PDAC are still unsatisfactory, and the survival rate is still less than 10%. Early diagnosis is the key to the treatment of pancreatic cancer (23). However, the etiology of pancreatic cancer is still unclear and there is a lack of effective biomarkers for early diagnosis. Therefore, it is urgent to explore the genetic or epigenetic factors of the occurrence and development of pancreatic cancer to identify new therapeutic targets or biomarkers. RNA methylation is one of the most important epigenetic modification in RNA post-transcriptional regulation and has been a research frontier and hotspot recently. RNA methylation includes several types such as m6A, m1A, m5C, m7G, and Nm, among which RNA m5C methylation refers to the modification of the fifth cytosine of RNA by methylation (24)(25)(26). At present, relevant studies have shown that RNA m5C methylation widely existed in cells and played an important role in various physiological and pathological processes (27)(28)(29). In this study, we downloaded 177 PDAC patients' gene expression profiles from the TCGA database and constructed an 8-m5Crelated lncRNA risk model. As far as we know, this is the first prognostic analysis of the m5C regulator-related lncRNAs in PDAC.
Long non-coding RNAs are more than 200 nucleotides, which were initially considered as "junk sequence" without specific biological functions. However, recent reports have proved that lncRNAs were widely expressed in human cells and closely related to the occurrence and development of tumors, which suggested that lncRNAs may be potential prognostic markers and therapeutic targets for malignant tumors (30). In addition, lncRNAs could regulate tumor growth in several ways including in situ regulation and molecular sink, combining with protein and molecular scaffold. Among them, the most essential way was regulating related gene expression (31,32). Hence, we focused on the lncRNAs that had co-expression relationship with m5C regulators in PDAC and used bioinformatic and statistical methods to generate a prognostic risk model for PDAC. In our research, we compared the expression of 13 m5C regulator genes in pancreatic cancer and normal tissues via the GEPIA online tool and identified a total of 242 m5C-related lncRNAs depending on the co-expression relationship. Furthermore, we analyzed the PDAC patients' clinical information together with gene transcriptome data downloaded from TCGA. Finally, 17 prognostic m5C-related lncRNAs were detected and all the prognostic m5C-related lncRNA expression was significantly different between pancreas normal and tumor tissues. In addition, LASSO regression was performed, and we constructed an 8-m5Crelated lncRNA risk model and separated all the 177 PDAC patients into high-and low-risk groups based on the risk scores. The PCA and patients' overall survival analysis confirmed that the 8-m5C-related lncRNA risk model could predict PDAC patient prognosis effectively. Typically, clinicians would make a rough judgment on PDAC patients' long-term survival based on clinical indicators such as general patient profile, pancreatic tumor size, invasion of peripheral blood vessels and nerves, presence of metastasis, and tumor pathological type. However, with advances in cancer gene research technology, scientists began to pay more attention to the changes in gene expression levels during tumor progression (33,34). In our work, univariate and multivariate Cox regression analysis indicated that the risk scores depending on the m5C-related lncRNAs expression were the independent risk factor for PDAC and better than traditional indexes including age, tumor grade, and stage on predicting patient survival rate. The currently constructed PDAC prognostic risk model mainly focused on a single molecular biomarker. However, only one biomarker could not work as a precise tumor biomarker in clinical diagnosis such as the typical PDAC biomarker CA19-9 that was not aberrant in every PDAC patient. It has been reported that lncRNAs could be released in serum, saliva, and urine (35). Hence, it has the potential to be a stable biomarker in human blood or other body fluid. A lot of studies have demonstrated that lncRNAs could work as a new type of biomarker in several malignant tumors including hepatocellular carcinoma, breast cancer, lung cancer, and pancreatic cancer (36). Nevertheless, our research was the first   (37)(38)(39). We checked the above lncRNA expression in the human pancreatic ductal cell line and three pancreatic cancer cell lines. The results showed that CASC8 was overexpressed in all pancreatic cancer lines, which indicated that CACS8 may also play a key role in promoting pancreatic cancer progression. Recent research showed that m5C-related genes could regulate tumor immune microenvironment in breast cancer (40). LncRNAs are also known to be involved in regulating tumor-infiltrating lymph cells in malignant tumors (41).
In our study, we carefully evaluated the relationship between m5C-related lncRNAs and PDAC-infiltrating immune cells. Our results exhibited that the distribution of tumor-infiltrating immune cells in PDAC were different in high-or low-risk groups. The immune scores and stromal scores of the high-risk group were also lower than the low-risk group. We found that naïve B cells, CD8+T cells, Treg cells, and resting NK cells had a higher expression level in the low-risk group. CD8+T cells, Treg cells, and NK cells were considered as tumor suppressors in malignant tumors. In pancreatic cancer, the depletion of Treg cells made the abnormal distribution of fibroblast, which could recruit myeloid cells and restore the immune suppression environment via secreting multiple cell chemokine including CCL3, CCL6, and CCL8 (42). In addition, the M0 and M2 phenotype macrophages had a higher expression level in the high-risk group. Based on the previous study, M2 phenotype macrophage has been proven as a co-carcinogenic factor in pancreatic cancer (43). Therefore, m5C-related lncRNAs may regulate pancreatic cancer progression via promoting M2 phenotype macrophage polarization or infiltration in PDAC. As we know, there were no published reports about detecting the potential function of m5C-related lncRNAs in regulating the distribution of tumor-infiltrating lymph cells in PDAC at present.
However, there are still some limitations in our study. For example, all the PDAC patients' transcriptome expression data and clinical information were downloaded and analyzed from the TCGA database or GEPIA online database. We only evaluated the m5C-related lncRNA risk model via internal validation, lacking external database validation. Nevertheless, we tried to download the PDAC patients' information from ICGC (International Cancer Genome Consortium) or GEO database, but the corresponding risk scores could not be calculated due to the lack of expression data of the selected m5C-related lncRNAs in these databases. In addition, we only detected the screened m5C-related lncRNA expression in vitro, and the in vivo experiments and PDAC patients' tumor tissue validation should still be performed in the future to make the prediction results more reliable. Moreover, the nomogram we constructed only calibrated within 3 years, because of the small number of PDAC patients and the lack of T1 stage patients, which suggested that our findings need further validation in larger PDAC patient follow-up cohorts. Additionally, there was a lack of molecular mechanism research in our study. Identifying which m5C-related lncRNAs could regulate pancreatic cancer survival is just the beginning. We will explore the specific mechanism of the screened lncRNAs that affected PDAC progression in our next work.

CONCLUSION
We constructed an 8-m5C-related lncRNA prognostic risk model for PDAC patients based on the transcriptome expression and clinical data from TCGA. The m5C-related lncRNA risk model was proved to have an independent prognostic value and provided an accurate survival prediction for PDAC patients. In addition, our research also offered us a better understanding of the regulation of tumor-infiltrating lymph cells in PDAC. In brief, the m5C-related lncRNA risk model may provide us the potential biomarkers or treatment targets for PDAC.

DATA AVAILABILITY STATEMENT
The datasets presented in this study can be found in online repositories. The names of the repository/repositories and accession number(s) can be found in the article/Supplementary Material.

ETHICS STATEMENT
All the patient's raw data in our study were downloaded from TCGA, so there are no ethical issues and other conflicts of interest. The patients/participants provided their written informed consent to participate in this study.

AUTHOR CONTRIBUTIONS
HY, JL, and LZ summarized and analyzed the TCGA patients' data. HY and PW drafted the manuscript. KJ and BX designed and monitored the research. All authors contributed to the article and approved the submitted version.