Comprehensive Analysis of the PROSER2-AS1-Related ceRNA Network and Immune Cell Infiltration in Papillary Thyroid Carcinoma

Introduction

Recently, the morbidity of thyroid cancer has increased dramatically, accounting for 1–2%^1,2 of all cancers, and papillary thyroid carcinoma (PTC) is the most common histological type.³ Currently, the treatment of thyroid cancer includes surgery, hormone therapy and radioiodine therapy. Most PTC patients have a good prognosis; however, PTC patients with distant metastasis have poor prognosis, and available treatments cannot dramatically improve prognosis.⁴ It has been reported that the mortality of thyroid cancer in China is increasing gradually at a rate of 1.4% per year.⁵ Fine needle biopsy is a reliable preoperative qualitative method for thyroid cancer, but it is easy to misdiagnose and to miss a diagnosis when the volume of nodules is small. Because of early lymph node metastasis and postoperative recurrence, patients with PTC face a poor prognosis. Therefore, early diagnosis and comprehensive treatment of PTC are particularly important.

Long noncoding RNAs (lncRNAs), with a length of more than 200 bases, are unable to be translated into proteins. Growing data suggest that lncRNAs are involved in tumor development processes, including pretranscription, posttranscription, cell proliferation, differentiation, apoptosis and migration.^6,7 Thus, lncRNAs can act as new diagnostic biomarkers and therapeutic targets in cancer.⁸ lncRNAs associated with thyroid cancer are mainly divided into two categories: oncogenes and tumor suppressor genes. Among these lncRNAs, the main oncogenic lncRNAs are BANCR, HOTAIR, H19 and MALAT1,^9–12 and the main tumor-suppressing lncRNAs are PTCSC3, PTCSC2 and MEG.^13–16 Competing endogenous RNA (ceRNA) was proposed by Salmena et al¹⁷ who described a hypothesis that lncRNAs could induce gene silencing by binding to miRNAs. The ceRNA network plays crucial roles in tumorigenesis.^17,18 As a result, it is imperative to investigate the role of ceRNA networks in the poor prognosis of PTC.

The tumor microenvironment (TME) is composed of extracellular matrix (ECM) and nonmalignant normal cells, including fibroblasts, pericytes, immune cells and endothelial cells.^18–20 Tumor-infiltrating immune cells (TIICs) have emerged as a key factor that drives tumor initiation and metastasis.²¹ A recent study revealed that the occurrence and development of thyroid cancer are highly correlated with inflammation and immune cell infiltration, and exploring the immune cell infiltration pattern is helpful for judging the prognosis and treatment response of patients.²²

In this study, we analyzed the regulatory ceRNA network of lncRNAs to further explore the molecular mechanism of the occurrence and development of PTC to identify more biological diagnostic and prognostic monitoring markers of PTC. We revealed a potential role for PROSER2-AS1 in tumor immunology and its prognostic value, which will help to further explore the profound effects of lncRNAs in the process of cancer.

Materials and Methods

Data Acquisition and Screening of DELs

The gene expression profile of GSE66783 (Platform: GPL19850) was selected from the Gene Expression Omnibus (GEO) (http://www.ncbi.nlm.nih.gov/geo/).²³ The GSE66783 dataset included 5 pairs of PTC and normal thyroid tissues, and normalized microarray data were reanalyzed using the GEO2R online tools (https://www.ncbi.nlm.nih.gov/geo/geo2r/).²⁴ The differentially expressed genes (DEGs) were considered significant when |logFC|>1, with an adjusted P value<0.05. A volcano plot was generated with the “ggplot2” package of R software. Differentially expressed lncRNAs (DELs) were identified by setting the cutoff value at |logFC| > 2 and adjusted P<0.05. The heatmap of the DELs was drawn by the complex Heatmap package of R software. We also gathered RNA-sequencing data of PTC from the TCGA database (https://genome-cancer.ucsc.edu/) and clinical information from 510 samples of tumors and 58 samples of paracancerous tissue. Our study procedure flowchart is shown in Figure 1. As all data in our study were obtained from public databases (GEO and TCGA database) following publication guidelines, the Ethics Committee of the First Affiliated Hospital of Shandong First Medical University & Shandong Provincial Qianfoshan Hospital required no further ethical approval for the study.

Figure 1 Flow diagram of data collection and analysis.

Univariate and Multivariate Cox Analyses

Univariate Cox regression analysis was used to determine the effect of DELs on progression-free survival (PFS) using the survival package of R software. To identify independent prognostic factors, all significant variables on univariate Cox regression analysis (P≤0.05) were subjected to multivariate Cox regression analysis.

Survival Analysis

The prognostic analysis was performed using the R survival package (Version 3.6.3). The survival result was displayed as a Kaplan–Meier (K-M) curve with the cutoff value determined by survminer and visualized by survival (Version 3.6.3). A value of p <0.05 was considered statistically significant.

Gene Set Enrichment Analysis (GSEA)

To identify the potential function of PROSER2-AS1, gene set enrichment analysis (GSEA) was conducted between the low-and high-expression groups using the R package cluster Profiler (version 3.6.3). Gene sets with false discovery rate (FDR)< 0.25 and P.adjust <0.05 were considered to be significantly enriched.^25,26

Immune Cell Infiltration Analysis

We conducted an immune cell infiltration analysis to evaluate the correlation between PROSER2-AS1 expression and the infiltrating levels of 24 immune cell types²⁷ on the basis of TCGA database and performed using the single-sample GSEA (ssGSEA) method and the R package GSVA (version 3.6.3).²⁸ Moreover, the relationships between PROSER2-AS1 expression and several immune cell biomarkers were estimated using the GEPIA database (http://gepia.cancer-pku.cn/index.html).²⁹ The immune cell composition in PTC was obtained from the TCIA database (https://www.tcia.at/home).³⁰

Construction of the ceRNA Network

Online tools of lncATLAS (http://lncatlas.crg.eu/)³¹ and LncSLdb (http://bioinformatics.xidian.edu.cn/lncSLdb/index.jsp)³² were used for cellular localization prediction of PROSER2-AS1. To confirm the ceRNA network, we screened the target miRNAs of PROSER2-AS1 based on lncRNASNP2 (http://bioinfo.life.hust.edu.cn/lncRNASNP#!/).³³ The mRNAs targeted by miRNAs were predicted using miRWalk3.0 (http://mirwalk.umm.uni-heidelberg.de/).³⁴ Predicted genes that were fitted to the miRanda, miRDB, and TargetScan databases were considered as target genes. We constructed a ceRNA network using Cytoscape software (Version 3.7.1).

Results

Identification of DELs in PTC Patients

We first obtained the gene expression profiles from the GSE66783 dataset, including 5 tumor samples and 5 paired normal thyroid samples. Volcano plots of the DERNAs (including DElncRNAs, DEmiRNAs, and DEmRNAs) were visualized using the ggplot2 package of R software (Figure 2A). Heatmap clustering was drawn by the Complex Heatmap package of R software (Figure 2B). By setting the cutoff value at |logFC| > 2 and adjusted P<0.05, we identified 105 DELs. We identified PROSER2-AS1 by conducting univariate and multivariate Cox regression analyses based on TCGA database (P≤0.05) (Table 1). In the univariate Cox regression, different levels of PROSER2-AS1 showed statistical significance with the variables (P=0.007). Multivariate analysis further revealed that PROSER2-AS1 (P<0.05) was an independent prognostic factor in the PFS of PTC patients. Figure 2D showed that multiple cancer types exhibited different PROSER2-AS1 expression, including breast, kidney, stomach, liver, uterus, prostate, brain, colon, rectum and thyroid cancer. However, the changing trends of PROSER2-AS1 expression in the above cancers were different. In thyroid cancer, the expression of PROSER2-AS1 was significantly reduced, as validated by TCGA database analysis (Figure 2C). The immune-related gene list was downloaded from ImmPort (https://www.immport.org/) (Table S2).

Table 1 Univariate and Multivariate Cox Regression Analyses of Clinical Characteristics Associated with PFS in PTC

Figure 2 Differentially expressed genes in thyroid papillary carcinoma and normal samples from GSE66783 datasets. (A) Volcano plot of DEGs in GSE66783. (B) Heatmap of DELs in GSE66783. (C) The expression of PROSER2-AS1 in thyroid papillary carcinoma and paired normal thyroid tissues. (D) The expression levels of PROSER2-AS1 in different types of tumors based on TCGA database (*P<0.05, **P<0.01, ***P<0.001). Abbreviations: DEGs, differentially expressed genes; DELs, differentially expressed lncRNAs.

Association Between PROSER2-AS1 Expression and Prognosis in PTC

Correlations between PROSER2-AS1 and clinical parameters were examined using the ggplot2 package (Table 2). There were significant differences in T stage (P=0.015), N stage (P<0.001) and pathologic stage (P=0.023) between the PROSER2-AS1 low expression group and high expression group. We further used Kaplan–Meier curves to estimate the correlation of PROSER2-AS1 expression with the prognosis of PTC. The results showed that the lower expression group of PROSER2-AS1 had a poorer prognosis of PTC (PFS: HR=0.48, 95% CI=0.28–0.82, P=0.01) (Figure 3A). Additionally, subgroup analysis of KM curves by different clinical features using R software further revealed the association between PROSER2-AS1 expression and prognosis in PTC. We found that its decline was significantly linked with poor prognosis in the subgroups of T stage (PFS: HR=0.49, 95% CI=0.28–0.86, P=0.016), N stage (PFS:HR=0.41, 95% CI=0.20–0.87, P=0.011), M stage (PFS:HR=0.46, 95% CI=0.20–1.04, P=0.03), pathological stage (PFS:HR=0.48, 95% CI=0.23–1.00, P=0.049), age (>45 years) (PFS:HR =0.46, 95% CI=0.23–0.93, P=0.024) and sex (female)(PFS: HR=0.44, 95% CI=0.22–0.86, P=0.028), as shown in Figure 3B–G. The results demonstrated that PROSER2-AS1 expression significantly correlated with poorer outcome in PTC. ROC analysis of PROSER2-AS1 expression revealed that the area under the ROC curve was 0.820, which suggests that PROSER2-AS1 may be a potential diagnostic biomarker of PTC (Figure 3H).

Table 2 Correlation Between PROSER2-AS1 Expression and Clinicopathological Characteristics in PTC

Figure 3 PROSER2-AS1 expression and prognosis in patients with PTC. (A) Kaplan–Meier (K-M) curves were drawn using the R package survminer to evaluate the prognostic value of PROSER2-AS1 in the PFS of PTC patients. (B–G) K-M curve subgroup analysis of T stage, N stage, M stage, pathological stage, sex and age. (H) ROC curve showing the efficiency of the PROSER2-AS1 expression level in distinguishing PTC tissue from nontumor tissue.

Functional Enrichment Analysis

To predict the function of PROSER2-AS1, gene set enrichment analysis (GSEA) was used to search for KEGG and Reactome pathways using TCGA database. The results revealed that the immune regulatory interactions, antigen activating B cell receptor, FCGR activation, signaling by the B cell receptor and the role of phospholipids in phagocytosis were significantly enriched (Figure 4A). In addition, cytokine receptor interactions, antigen processing and presentation, the JAK-STAT signaling pathway, the cell cycle and the p53 signaling pathway were significantly enriched by KEGG pathway analysis (Figure 4B). These results suggested that PROSER2-AS1 is associated with many malignancy-related pathways in thyroid cancer, especially immune-related pathways.

Figure 4 Enrichment plot from the GSEA. Gene set enrichment plots of the (A) REACTOME pathway and (B) KEGG pathway with low PROSER2-AS1 expression in PTC. Abbreviations: NES, normalized NS; Padj, adjusted P value; FDR, false discovery rate.

Correlation Between PROSER2-AS1 and Immune Cell Infiltration

We explored the association between PROSER2-AS1 expression and infiltrating immune cells in PTC using the GSVA package of R software. As shown in Figure 5A, the results revealed that the levels of many infiltrating immune cell types, such as dendritic cells, macrophages, neutrophils, Th1 cells, Th2 cells and Treg cells, decreased significantly in the high-expression group of PROSER2-AS1. Next, we further analyzed the relationship between PROSER2-AS1 and the infiltrating immune cell types. PROSER2-AS1 expression was negatively correlated with the infiltration levels of B cells (r=−0.130, P=0.004), neutrophils (r=−0.240, P<0.001), dendritic cells (r=−0.260, P<0.001), macrophages (r=−0.220, P<0.001), Th1 cells (r=−0.220, P<0.001), Th2 cells (r=−0.320, P<0.001) and Treg cells (r=−0.300, P<0.001) and positively correlated with NK cells (r=0.170, P<0.001) (Figure 5B–I). Similar results were verified in Table 3. PROSER2-AS1 was significantly negatively correlated with M1 macrophage biomarkers (NOS2, IRF5, and PTGS2), neutrophil biomarkers (CEACAM8, ITGAM and CCR7), dendritic cell biomarkers (HLA-DPB1, HLA-DQB1, HLA-DRA, HLA-DPA1, CD1C, NRP1 and ITGAX), and Th1 biomarkers (STAT4 and STAT1) in PTC. B cells (0.98%), dendritic cells (0.34%), M1 macrophages (2.40%), M2 macrophages (4.12%), natural killer cells (0.46%), neutrophils (11.90%), CD4⁺ T cells (1.50%), CD8⁺ T cells (1.06%), and regulatory T cells (1.15%) comprised the known immune cell components in PTC (Figure 5J). These results suggested that PROSER2-AS1 is involved in tumor immune cell infiltration.

Table 3 Correlation Between PROSER2-AS1 and Biomarkers of Immune Cells in PTC

Figure 5 Correlation of PROSER2-AS1 expression with immune cell infiltration in PTC. (A) Twenty-four immune cell infiltrates in different PROSER2-AS1 expression groups of PTC patients. (B–H) PROSER2-AS1 had significant negative correlations with the infiltrating levels of B cells, neutrophils, dendritic cells (DCs), macrophages, Th1 cells, Th2 cells and Treg cells. (I) PROSER2-AS1 had a positive correlation with infiltrating levels of NK cells. (J) Immune cell infiltration in papillary thyroid carcinoma. *P<0.05, **P<0.01, ***P<0.001.

Construction of the ceRNA Network

We examined the subcellular localization of PROSER2-AS1 and found that it was mainly localized in the cytoplasm (Figure S1). To determine the underlying mechanisms, we constructed a lncRNA-miRNA-mRNA network. The prediction of lncRNA-miRNA and miRNA-mRNA functional pairs was carried out based on LncRNASNP2 and miRWalk3.0, respectively. We identified 28 pairs of lncRNA-miRNA (Figure S2) and 270 pairs of miRNA-mRNA (Table S1). A ceRNA network centered on PROSER2-AS1 was visualized using Cytoscape, as shown in Figure 6.

Figure 6 The PROSER2-AS1 related ceRNA network in PTC. The yellow triangle represents lncRNA PROSER2-AS1. The green diamond represents miRNA. The blue triangle represents mRNA.

Predicted PROSER2-AS1 Downstream Target Genes Associated with Immune Infiltration

The heatmap of the highest ranked co-expressed genes of PROSER2-AS1 was presented in Figure 7A. We determined 10 genes (FURIN, FGF2, RORA, NR3C1, TGFBR3, VEGFA, FGFR1, TGFBR1, ACVR2A and MAP2K1) related to immunity by analyzing the intersection of PROSER2-AS1-targeting mRNAs, an immune-related gene list and co-expressed genes using a Venn diagram (Figure 7B). Finally, only 2 genes (TGFBR3 and VEGFA) remained according to the principle of differential expression (Figure 7C) and survival curve analysis (Figure 7D and E) based on TCGA database. The expression levels of TGFBR3 and VEGFA were positively correlated with the expression of PROSER2-AS1 (Figure 7F and G). We next focused on the relationship between TGFBR3 and immune infiltration in PTC. As shown in Figure 8A, aDC (P<0.001), DC (P<0.001), macrophages (P<0.05), Th2 (P<0.01) and Treg (P<0.001) cells exhibited lower levels in the high TGFBR3 expression group compared to the low group. However, NK (P<0.001), T helper cells (P<0.01), Tcm (P<0.001) and Tgd (P<0.001) cells showed higher proportions in the high expression group. TGFBR3 had a negative relationship with infiltrating levels of aDC cells (r = −0.214, p < 0.001), DC cells (r = −0.224, p < 0.001), macrophages (r = −0.135, p =0.002), Th2 (r = −0.186, p < 0.001) and Treg cells (r = −0.284, p < 0.001) (Figure 8B–F), and positively correlated with T helper cells (r = 0.141, p = 0.001), Tcm (r = 0.218, p < 0.001), Tgd (r = 0.226, p < 0.001) and NK (r = 0.209, p < 0.001) in PTC (Figure 8G–J).

Figure 7 Expression of PROSER2-AS1 target mRNA related to immune cell infiltration in PTC. (A) Heatmap of the top 50 co-expressed genes of PROSER2-AS1. (B) Overlap between the immune-related gene list, the PROSER2-AS1 target mRNA and co-expressed genes of PROSER2-AS1 shown in Venn diagrams. (C) mRNA gene expression of selected genes from (B) in PTC based on TCGA database. (D and E) KM curves of VEGFA and TGFBR3 in PTC patients. (F and G) PROSER2-AS1 has significant positive correlations with VEGFA and TGFBR3. ***P<0.001.

Figure 8 Correlation of TGFBR3 expression with immune cell infiltration in PTC. (A) Twenty-four immune cell infiltrates in different TGFBR3 expression groups of PTC patients. (B–F) TGFBR3 had significant negative correlations with the infiltrating levels of aDC, dendritic cells (DCs), macrophages, Th2 cells and Treg cells. (G–J) TGFBR3 had a positive correlation with infiltrating levels of T helper cells, Tgd cells, Tcm cells and NK cells. *P<0.05, **P<0.01, ***P<0.001.

Discussion

With the application of high-resolution ultrasound and increasing attention to health, the incidence rate of thyroid carcinoma is increasing yearly. Papillary thyroid carcinoma (PTC) accounts for the majority of thyroid tumors. The prognosis of PTC is satisfactory; however, because of its biological characteristics, PTC is prone to spread via early lymph node metastasis or even distant metastasis, resulting in poor prognosis. Therefore, it is important to identify PTC accurately and early. In recent years, the biological function of lncRNAs has attracted much attention. Previous studies have shown that lncRNAs play a role in the occurrence and development of tumors, including papillary thyroid carcinoma.^35–38 In this study, differentially expressed lncRNAs were identified through the GEO database. Among them, PROSER2-AS1 was identified by univariate and multivariate Cox analyses. The expression of PROSER2-AS1 was not consistent in different tumors. Through Pan Cancer analysis in TCGA database, PROSER2-AS1 expression was up-regulated in bladder urothelial carcinoma (BLCA), cholangio carcinoma (CHOL), colon adenocarcinoma (COAD), rectum adenocarcinoma (READ), and stomach adenocarcinoma (STAD) and down-regulated in breast invasive carcinoma (BRCA), glioblastoma multiforme (GBM), head and neck squamous cell carcinoma (HNSC), kidney chromophobe (KICH), kidney clear cell carcinoma (KIRC), prostate adenocarcinoma (PRAD) and uterine corpus endometrial carcinoma (UCEC). In thyroid carcinoma, its expression was down regulated, which suggested that it is a suppressor gene in thyroid cancer. Kaplan–Meier survival analysis implied that PTC patients with low expression of PROSER2-AS1 had shorter PFS (p <0.05). Multivariate Cox analysis confirmed that low expression of PROSER2-AS1 was an independent risk factor for FPS in PTC patients. In addition, K-M subgroup analysis showed that clinicopathologic features were also associated with poor prognosis in PTC, including advanced stage and lymph node metastases. ROC analysis verified the diagnostic value of PROSER2-AS1 in PTC. To further investigate the functions of PROSER2-AS1 in PTC, we performed GSEA based on TCGA database. The GSEA results showed that cytokine receptor interaction, antigen processing and presentation, JAK-STAT signaling pathway, cell cycle, p53 signaling pathway and immune regulatory interactions were differentially enriched in PROSER2-AS1-low cells. These results suggest that PROSER2-AS1 may be considered a PTC prognostic biomarker.

lncRNAs have unique advantages compared with ordinary genes, including differential expression in different tissues, diseases, and various stages, which suggests that lncRNAs may be biomarkers for the specificity of a certain tissue or a definitive stage of a disease. Increasing evidence has shown that lncRNAs play roles in the development of many diseases, including cancers.^39–41 There are two main mechanisms of lncRNAs: transcriptional regulation and posttranscriptional regulation. The former plays a regulatory role in the nucleus by interacting with nuclear factors, while the latter plays a regulatory role in the cytoplasm through a ceRNA mechanism. In this study, PROSER2-AS1 was predicted to be located in the cytoplasm. To further study the mechanism, we established a PROSER2-AS1-related ceRNA network.

To date, accumulating evidence suggests that lncRNAs play important roles in the regulation of immune-related processes.^42,43 Therefore, we further analyzed the relationship between PROSER2-AS1 and immune cell infiltration. PROSER2-AS1 expression was negatively correlated with the degrees of macrophage, B cell and neutrophil infiltration in PTC. These results implied that low-PROSER2-AS1 patients might have more active immune responses than high-PROSER2-AS1 patients. In the following study, we selected target mRNAs from the ceRNA network by taking the intersection of immune-related genes. Finally, we identified two mRNAs (TGFBR3 and VEGFA) by screening the differential expression and the effect on prognosis in TCGA database. Correlation analysis showed that there was a positive correlation between PROSER2-AS1 and the two target mRNAs. Vascular endothelial growth factor A (VEGFA) has been reported to promote angiogenesis previously. In recent years, the role of VEGFA in tumor immunity received much attention, suppression of T-cell responses⁴⁴ and induction of Treg proliferation,⁴⁵ resulting in tumor-induced immunosuppression. However, these results were not consistent with our bioinformatics analysis.

We next focused our attention on TGFBR3. Transforming growth factor beta receptor III (TGFBR3) is a coreceptor of the TGF-β superfamily. As a suppressor, TGFBR3 plays a role in tumor metastasis.^46,47 Brent et al provided novel findings of decreased-expression of TGFBR3 suppressing tumor antigen-specific immune responses in both breast cancer and melanoma.⁴⁸ In our study, we observed low TGFBR3 expression in PTC and low expression of TGFBR3 was associated with a poor prognosis, suggesting anti-oncogenic role of TGFBR3 in PTC. Further analysis revealed that there was positive correlation between TGFBR3 expression and NK cell infiltration and negative correlation between TGFBR3 expression and Treg cell infiltration in PTC. However, few studies have focused on the role of TGFBR3 in PTC. Additional studies are needed to clarify the impact of TGFBR3 on the metastasis of PTC.

Our study has the following limitations. First, the study was performed primarily based on bioinformatic analyses, which need further experimental verification. Second, the study was a retrospective study, which needs to be supported by the results of further prospective studies. Finally, there was no evaluation of the potential diagnostic role, and it needs to be further explored whether PROSER2-AS1 can be used as a diagnostic marker. As a result, follow-up studies are needed to further validate our findings.

Conclusion

The results of this study indicate that reduced expression of PROSER2-AS1 is correlated with clinical outcomes in PTC and provide evidence for the role of PROSER2-AS1 in immune cell infiltration. TGFBR3 maybe the downstream gene of PROSER2-AS1. In summary, PROSER2-AS1 may be an important regulator of immune cell infiltration and a potential prognostic biomarker in PTC patients.