Development an Immune-Associated LncRNA Prognostic Index for Papillary Thyroid Cancer


 BACKGROUND: Papillary thyroid cancer (PTC) is the most common pathological type of thyroid cancer and its incidence is still on the rise worldwide. Our study aimed to provide an understanding of papillary thyroid cancer in immunological aspects and prompt its clinical diagnosis and treatment.METHODS: We obtained RNA-FPKM data and clinical data of papillary thyroid cancer (PTC) samples from the TCGA data portal. Immune-related genes (IRGs) were selected from Molecular Signatures Database v7.0. The immune-related lncRNAs were identified by constructing DELs-IRGs co-expression networks. Cox regression and Kaplan-Meier methods were used to construct immune-related lncRNAs prognosis index (IRLPI) and analyze the impact of IRLPI on overall survival. Gene Set Enrichment Analysis (GSEA) was performed by clusterProfiler R package to explore the potential molecular mechanism in TCGA dataset. The relationship between immune cell types and IRLPI were explored by using Pearson correlation analysis.RESULTS: A total of 478 primary PTC samples were collected for the analysis. 162 immune-related lncRNAs were identified in immune-lncRNAs co-expression network. And 7 lncRNAs were extracted to develop the IRLPI. Patients with high-IRLPI showed poor survival probability (P =0.00052) and IRLPI was identified as an independent prognostic factor (HR: 2.329, 95%CI: 1.481-2.664, p =2.54e-4). ROC curve declared a promising prognosis ability of IRLPI. PCA proved that IRLPI reflected different immune environment and biological processes. GSEA indicated that immune system development is the main different biological process between high- and low-IRLPI group. Furthermore, IRLPI was proved that it related with multiple immune cells and immune responses by Pearson correlation analysis.CONCLUSION: Our results showed that immune-related lncRNAs prognosis index (IRLPI) could serve as a potential prognostic model in papillary thyroid cancer and play an important role in the immunotherapies and surveillance of PTC.


Introduction
Papillary thyroid cancer (PTC) was the most type of thyroid cancer and still increased over the past decades [1]. Statistics showed that the incidence of the thyroid cancer increased by 4.1% for annual and kept a high death in recent years [2]. The standardized protocols were used for the treatment of PTC [3].
However, some of PTC patients still had a poor prognosis because of various risk factors [4]. Therefore, novel therapeutic strategies and biomarkers were urgently established to ameliorate the prognosis of PTC.
Growing studies revealed the relationship of between human immune and tumor appeared complicated [5,6]. Tumor cells could evade immune elimination to suppress immune response when tumor microenvironment was disordered [7][8][9]. Bai Y et al. reported that PTC had a strong correlation with the immune molecule BRAF V600E, PD-L1, and PD-1 [10]. Further, Gunda V et al. performed immunotherapy on PTC mice by combining BRAF inhibitors and anti-PD1/PD-L1 antibodies, suggesting that immunotherapy could improve the survival rate of PTC mice [11]. These existed evidences indicated that immune process played a crucial role in the development of cancer.
LncRNAs had a wide range of biological function and played a crucial regulatory factor in immune processes [12]. Pei et al. reported that lncRNA SNHG1 could inhibit the differentiation of Treg cells by promoting miR-448 expression and reducing IDO level, thereby impeding the immune escape of breast cancer [13]. On the other hand, Lnc-EGFR was proved that it could stimulate Treg differentiation, suppress CTL activity and then promote hepatocellular carcinoma (HCC) growth in an EGFR-dependent manner [14]. In addition, lncRNA-AFAP1-AS1 was also found to be positively correlated with the expression of PD-1, while both the high expression of AFAP1-AS1 and PD-1 indicated poor prognosis of nasopharyngeal carcinoma [15]. These studies revealed lncRNAs were associated with immune processed and prognosis of various cancers. However, the relationship between lncRNAs with the immunity and prognosis of PTC was still poorly studied and needed further elucidation.
In the present study, comprehensive analysis was performed to explore the relationship between lncRNAs and the immune of PTC and identify the immune-related lncRNA associated with the development of PTC. An immune-related lncRNAs prognosis index (IRLPI) was further constructed for developing a novel individualized prognostic biomarker of PTC by using the immune-related lncRNAs. Moreover, we used bioinformatics analyses to explore underlying mechanisms of IRLPI. Results from our study could offer foundation for subsequent personalized diagnosis and treatment of PTC.

Data preprocessing
The RNA-FPKM data and clinical data of PTC samples were obtained from the TCGA data portal (https://cancergenome.nih.gov/). In the study, the samples with no follow-up information or follow-up time less than 30 days were excluded. Next, 478 primary PTC samples and 58 non-tumor tissues were included for further analysis. Additionally, the expression pro les of GSE64912 and GSE83520 were downloaded for validation set through GEOquery R package. The GSE64912 included 18 primary PTC samples and 4 non-tumor tissues, and GSE83520 included 12 primary PTC samples and 12 non-tumor tissues.

Obtainment of lncRNA and mRNA symbol
The Ensembl IDs were converted to gene symbols based the GRCh38.98 genome le which was downloaded from the Ensembl database(http://asia.ensembl.org/index.html). The Ensembl IDs with maximum mean were reversed when more than one Ensembl IDs had a matched gene symbol. lncRNA and mRNA symbols were extracted according to gene type. 11,877 IncRNAs with overlapping three datasets were retained for this study. Next, a total of 332 Immune-related genes (IRGS) were acquired from Molecular Signatures Database v7.0 (Immune system process M13664.Immune response M19817). Only 325 IRGs with overlapping three datasets were extracted for further analysis.

Identi cation of immune-related lncRNAs in TCGA dataset
To identify immune-related lncRNAs which were played important roles in the pathogenesis of PTC, we screened out the differentially expressed lncRNAs (DELs) in TCGA dataset through the edgeR R package. The cut-off was fold changes >1 and FDR (the adjusted p value) < 0.05. The immune-related lncRNAs were identi ed by constructing DELs-IRGs co-expression networks through the Pearson correlation analysis. The P < 0.05 and Pearson coe cient> 0.35 was were considered statistically signi cant. To make the correlation analysis more accurate, the low-expression DELs and IRGs with the average expression less than 0.5 were ltered before the Pearson correlation analysis.
Validation of Immune-related lncRNAs in GSE dataset GSE64912 and GSE83520 were combined and performed to batch correction through sva R package. To validate immune-related lncRNA co-expression network, the Pearson correlation analysis and the same cut-off criteria were conducted to the GSE dataset.

Development of the immune-related lncRNAs-based prognostic index
To develop the immune-related lncRNAs-based prognostic index (IRLPI), We must rstly calculate the risk score of each sample. The risk score was conducted for previously studies [16,17]. The risk score of each patient was calculated as follows: Risk score = βgene1 × exprgene1 + βgene2 × exprgene2 + ⋯ + βgenen × exprgenen. Next, the IRLPI was Identi ed as log2(Risk score) in our study.
Speci cally, a log2 (FPKM value+ 1) data format was used for the development of IRLPI. The immunerelated lncRNAs were conducted the univariate Cox analysis for screening out survival-associated lncRNAs. The immune-related lncRNAs with p < 0.01 were identi ed as the survival-associated lncRNAs and selected for constructing the IRLPI by using the multivariate Cox analysis. To streamline factors and improve accuracy, the step function was performed to analysis in the R platform and the lncRNAs with p > 0.05 in calculating risk score were removed in our result. Kaplan-Meier curve analysis was further performed to evaluate the impact of IRLPI on overall survival (OS). The cut-off value was determined by its median value. To identify IRLPI was an independent factor, univariate and multivariate Cox regression analysis were conducted on combination with clinical indicators.
Validation of the 7 Immune-related lncRNAs expression In order to verify the differential expression of the 7 immune-related lncRNA in PTC, we chose GSE83520 and GSE64912 to verify the expression differences of the 7 immune-related lncRNA in PTC. GSE83520 and GSE64912 was performed using the t-test or paired t-test.

IRLPI and tumor microenvironment
To explore the potential different immune microenvironment of between high and low-IRLPI, we quanti ed the enrichment levels of tumor microenvironment based on 29 immune gene set by using the single-sample gene-set enrichment analysis (ssGSEA) [18,19]. The Pearson correlation was conducted to further evaluate the relationship of between each component and IRLPI. P<0.05 and /R/ > 0.15 indicated that correlation existed between IRLPI and its component.

Statistical analysis
All statistical analyses were performed in R 3.6.1 (http://www.r-project.org/) and its corresponding R package. Pearson correlation analysis was completed through Hmisc R package. Immune-lncRNAs coexpression networks were constructed using Cytoscape 3.7.2 software. The AUC of ROC curve was calculated by timeROC package R. The principal components analysis (PCA) were completed in R platform. To identify biological processes and signaling pathways that were differentially activated between the high and low-IRLPI group in PTC, we identi ed an ordered list of genes through edgeR R package and conducted Gene Set Enrichment Analysis (GSEA) on the gene with adjusted p < 0.05. The c5.bp.v6.2.entrez.gmt le came from Molecular Signatures Database (MSigDB, http://software.broadinstitute.org/gsea/index.jsp) was downloaded for the GSEA analysis to identify biological processes. GSEA analysis was completed by using clusterpro ler R package.

Results
Constructed immune-lncRNAs co-expression network A total of 353 lncRNAs were identi ed as the differentially expressed lncRNAs(DELs) in PTC based the edgeR algorithm. (Fig 1a-1b). After ltering the low-expression DELs and IRGs, 176 DELs and 246 IRGs were retained to construct the immune-lncRNAs co-expression network through Pearson correlation analysis. 162 DELs were nally identi ed as the immune-related lncRNAs of PTC (Fig 1c) in immunerelated lncRNAs co-expression network. The result of correlation analysis was validated in the GSE dataset(S1).
As the Figure 2b showed that patients with high-IRLPI were correlated with worse overall survival in the Kaplan-Meier analysis (P=0.00052). And the time-dependent ROC curve analysis declared that the IRLPI had a promising prognosis ability for OS (Fig 2c). The IRLPI distribution and these lncRNAs expression patterns of PTC were shown in Figure 2d. At multivariate analysis, the IRLPI still maintained independently associated with OS, with a HR of 2.329 (95%CI:1.481-2.664, p =2.54e-4) along with other features ( Table 2).
Validation of the 7 Immune-related lncRNAs expression As the Figure 4 shown that all of the 7 immune-related lncRNAs were also highly expressed in PTC in GSE datasets. The result was also consistent with the result of TCGA dataset.
High and low-IRLPI groups displayed different tumor immune microenvironment We used the principal components analysis (PCA) to explore the ability of the discrimination between low-and high-IRLPI groups based on the 7 immune-related lncRNAs. As the Figure 4 showed that the 7 immune-related lncRNAs could be used to better divide PTC patients into two different directions comparing with the all immune-related lncRNAs of PTC, immune-related mRNA and whole gene expression pro les. This result indicated PTC patients with different IRLPI existed differential immune environment. We found that IRLPI was exactly associated with tumor immune microenvironment in further analysis (Fig 4a -4b). The APC co stimulation, pDCs, Type I IFN response and Type II IFN response were negative associated with the IRLPI (Fig 5c-5f).

GSEA analysis
Various biological processes were associated with IRLPI (S2). One of these biological processes was immune system development associated with immune, which indicated that the high and low-IRLPI groups had different immune process (Fig 6a). As Figure 6b-6i shown that IRLPI was also associated with multiple tumor-associated signal pathways, including Metabolic pathway, Focal adhesion, Calcium signaling pathway, cGMP-PKG signaling pathway, Adrenergic signaling in cardiomyocytes, Carbon metabolism, Pancreatic secretion, Glutamatergic synapse.

Discussion
With the in-depth study of tumor immunotherapy, tumor immune had been found that it played an increasingly important role in the recurrence and metastasis of PTC [20,21]. LncRNAs were con rmed to serve as an important role in immune process activities, such as antigen exposure, antigen recognition, immune activation, immune cell in ltration and tumor clearance by regulating the expression of immune genes5,12. In our study, we de ned 325 immune-related lncRNAs closely associated with PTC. Moreover, many of these lncRNAs were impacted on the prognosis of PTC.
In further analysis, the IRLPI established by 7 immune-related lncRNAs was a promising biomarker closely related to the prognosis of PTC through complex analysis. It was an independent factor for predicting the prognosis of PTC. High IRLPI group predicted a poor prognosis compared with low IRLPI.
Many scholars also proposed different prognosis models for PTC in previous studies [22,23], but our index was focused on immune-related lncRNAs of PTC. This had not explored at previous. As the deepening of the research on lncRNAs, it was con rmed that lncRNAs played an important control factor in tumor immune [7]. Therefore, we believed that IRLPI model based on immune-related lncRNA would show stronger predictive power and value in PTC.
CRNDE and LINC01614 were considered to be oncogenes in previous studies. CRNDE could promote cell proliferation, invasion and migration by competitively binding miR-384, and the CRNDE/miR-384/PTN played an important role in the regulation of PTC progression [24]. The results were consistent with our analysis. In Colorectal cancer, CRNDE had been demonstrated as a novel serum-based marker for the diagnosis and prognosis of colorectal cancer [25]. Furthermore, CRNDE also had been shown to be elevated in glioma [26], hepatocellular carcinoma [27], lung cancer [28] and other tumors [29,30]. For LINC01614, it could promote expression of FOXP1 to impact the prognosis of lung cancer by inhibiting the expression of mir-217 [31]. The same was true for breast cancer [32]. LINC01614 could affect the survival of breast cancer associated with TGF -and focal adhesion kinase (FAK) signaling. Choi et al [33] also found that LINC01614 could enhance the proliferation, colony formation and migration ability of gastric cancer cells. Unfortunately, only CRNDE and LINC01614 of 7 immune-related lncRNAs had been explored in cancers so far. Our study broadened the understanding of between lncRNAs and PTC, and also gave us some new insight into immune aspect.
The 7 immune-related lncRNAs could use to better divide PTC patients into two different directions in our analysis. The result indicated that various of signal pathways and biological processes were associated with the IRLPI in our result. Activation of cGMP-PKG signaling pathway could inhibit the proliferation and induce apoptosis of colon cancer [34]. Wang et.al indicated activation of Calcium signaling pathway could induce autophagy to promote apoptosis of tumor [35]. Glutamatergic synapse could promote Ca2+ entry into cells or regulate the PI3K/AKT pathway to impact tumor [36,37]. In cancer patients, the nutrient metabolism pathways existed changes. Ciavardelli et al. indicated that there was metabolic heterogeneity between tumor cells and stromal cells [38]. In addition, adhesion and invasion-related pathways such as focal adhesion, could affect the invasion and metastasis of tumor. Previously studies were declared that the adhesion kinase was difference in follicular thyroid carcinoma with distant metastases [16,39,40].
Furthermore, IRLPI was associated with immune system development, which was a biological process with immune-related biological progress. This result manifested that IRLPI was related with immune status and immune microenvironment of PTC.
To explore the immune mechanism of IRLPI, ssGSEA analysis and Pearson correlation analysis were performed. Our result showed that IRLPI was closely associated with the PTC immune microenvironment. The in ltration of pDCs, APC and regulation of two types of IFN response were considered as anti-tumor component of IRLPI with negative correlation. pDCs might induce the differentiation of naive CD4+ T cells into ICOS+Foxp3+Tregs and was one of the mechanisms underlying tumor escape in PTC plus multinodular non-toxic goiter patients [41]. However, other studies showed a signi cant reduction in pDCs in cancer tissues [42]. T cells required at least two signals for their full activation in immune response. The interaction of between T cell receptor with an MHC antigenic peptide complex was the rst factor. The interaction of between T cells receptors and co-stimulatory molecules of antigen presenting cells (e.g.LFA-3, CD40, ICOS, CD80 or CD86) was the second factor [43,44]. APC co-stimulus molecules served as the second factor for anti-tumor in PTC by our analysis. Li et al [45] declared that IFN-a (type I IFN) could induce Fas expression and apoptosis in hedgehog pathway through inhibiting Ras-Erk signaling in Basal cell carcinoma, which indicated IFN-a was Anti-cancer agent. This was consistent with our result. There was controversial for IFN-γ (Type II IFN) in PTC. IFN-γ (Type II IFN) could mediate the activation of T lymphocytes and played an anti-cancer role [46]. And in PTC, IFN-γ could induce the secretion of CXCL10 and inhibit the proliferation of PTC cells signi cantly. But in funnily, LV et al [47] found that TNF-a and IFN-γ could induce epithelia mesenchymal transition of the three PTC cell lines and malignant progression in human PTC cells. The topic of between immune and tumor was complex and still unclear at present. In our study, we identi ed the anti-tumor immune components associated with IRLPI in PTC. Our preliminary analysis result could broaden our understanding of the immunological aspects of PTC and provide a perspective to explore for further. This research comprehensively explored the impact on between IRLPI established by 7 immune-related and prognosis of PTC in R platform. IRLPI was novel index for the prognosis of PTC. Meanwhile, we investigated the latent mechanism of IRLPI by applying multiple analytical methods. However, some limitations still needed to be considered in our study. First, the lack of validation in another independent cohort was a limitation of our study. Second, the reliability of molecular mechanism and immune-related analysis results faced challenged because of the lack of in vivo and in vitro experiments.

Conclusion
IRLPI established by 7 immune-related lncRNAs was constructed and anticipated to serve as a potential prognostic biomarker. We further systematically analyzed the role of between IRLPI and the monitoring of prognosis in PTC. Our ndings would provide a novel insight for immunotherapies and surveillance in PTC. Availability of data and materials R 3.6.1(http://www.r-project.org/) was an open source software. The RNA-FPKM data and clinical data of PTC samples came from the TCGA data portal (https://cancergenome.nih.gov/). GSE64912 and GSE83520 gene expression pro les were downloaded through GEOquery R package. The GRCh38.98 genome le could acquire from the Ensembl database(http://asia.ensembl.org/index.html). Immunerelated genes (IRGS) were acquired from Molecular Signatures Database v7.0 (Immune system process M13664.Immune response M19817)

Competing interests
The authors declare that they have no competing interests.