Analysis of the regulatory mechanisms of prognostic immune factors in thyroid cancer

Objective To explore the regulatory mechanism of immune prognostic factors in thyroid cancer. Methods Based on the TCGA database and GEO database, this study used bioinformatics methods to study the potential regulatory mechanism of thyroid cancer prognosis, analyzed the differentially expressed genes and differential miRNAs between thyroid cancer and normal paracancerous tissues by R software, and constructed lasso risk factors. The immune prognostic factors of thyroid cancer were obtained from the model, and the miRDB website was used to predict the possibility of differential miRNA target binding of the immune prognostic factors and correlation analysis was performed, and finally verified by cell experiments. Results There were 1413 differentially expressed genes between thyroid cancer and normal paracancerous tissues, among which 21 immune-related genes were prognostic factors with significant differences in expression; lasso risk model obtained AKAP12, APOC1, TIMP3, ADAMTS9, ANK2, HTRA3, SYNDIG1 , ADAMTS5 and DACT1 were nine prognostic factors. A total of 58 differential miRNAs were found in thyroid cancer tissues and non-cancerous tissues. The possibility of differential miRNA targeting and binding of immune prognostic factors on the miRDB website and cell experiments was analyzed. Conclusions The potential miRNA regulatory mechanism of immune prognostic factors in thyroid cancer has been explored.


Introduction
Thyroid cancer is a common endocrine malignancy (1,2).Although THCA has a low mortality rate and a relatively good prognosis, the existing therapeutic modalities for locally advanced or recurrent metastatic THCA are still ineffective in improving the prognosis of patients.Therefore, new therapeutic modalities targeting the molecular mechanisms of thyroid carcinogenesis and progression are still being explored.
Recent studies have revealed that normal cells must possess the intrinsic characteristics of tumor cells to develop into cancer cells; most importantly, tumor cells can evade antitumor immune responses and form tumor foci (3).Antitumor immunity is an essential mechanism for the host's inhibition of tumorigenesis.However, tumor cells can evade the body's immune attack, survive, and proliferate to form tumor foci (4,5).Tumor cells can play a crucial role in disease development through different immune mechanisms.Immune-related genes have been shown to play a role in cancer development (6)(7)(8) and, in some studies, can also can predict the prognosis of cancer patients and as therapeutic targets (9)(10)(11)(12).In addition to being an effective prognostic biomarker by constructing the immune prognostic features of cancer using bioinformatics, it is also a promising and novel therapeutic target.Wang et al. (13) found that an immune-related prognostic signature consisting of SLC10A2, FGF2, CCL28, NDRG1, ESM1, UCN, UTS2, and TRDC could effectively predict 3-and 5-year overall survival rates.Heterogeneous expression of immune-related genes can be used as a prognostic factor, but the regulatory mechanism for the heterogeneous expression of immune factors is unknown.
Micro ribonucleic acids (miRNAs) are conserved endogenous short-stranded noncoding ribonucleic acid (RNA) molecules consisting of 21-23 nucleotides.miRNAs have been shown to target multiple messenger RNAs (mRNAs) and are involved in almost every biological process (14, 15).In recent years, studies targeting the role of miRNAs in cancer based on miRNA expression profiles have increased annually.Many miRNAs are down-or up-regulated in human cancers and act as oncogenic or intra-cancer suppressors (16).In ovarian cancer, many miRNAs regulate the epithelial-mesenchymal transition program (17).At the level of molecular mechanisms, miR-125a and miR-125b can regulate immune cell development and function, thus acting as tumor suppressors or promoters (18).
In this study, a total of nine prognostic immune factors, AKAP12, APOC1, TIMP3, ADAMTS9, ANK2, HTRA3, SYNDIG1, ADAMTS5, and DACT1, were obtained by screening differentially-expressed genes and differential miRNAs using LASSO analysis and analyzing the differential miRNA.The target-binding relationship between differential miRNAs and prognostic immune factors was also analyzed.We validated them using correlation analysis to find the regulatory mechanism of immune prognostic factors in THCA.

Data sources
The data sources for this study were the TCGA and GEO databases.Those databases included RNA sequencing (RNAseq) data and clinicopathological information from 510 THCA samples; mRNA data from 58 normal tissues adjacent to the cancer were obtained.The tissue from 29 THCA patients and five non-THCA patients was obtained from the GEO database GSE103996 dataset miRNA expression profile.Immune-related genes was got from the Immunology Database and Analysis Portal database (https://www.immport.org).

Differential expression analysis
Differential expression genes(DEGs) and differentiallyexpressed miRNAs of THCA tissues and paraneoplastic tissues were analyzed differentially using the R package limma package, where differentially-expressed genes were screened at p < 0.05, | log2FC|>1, and differentially-expressed miRNAs were screened at p < 0.05, |log2FC|>1, and represented in volcano and heat maps using Venn diagrams to take the intersection of the screened DEGs with immune-related genes.

Prognostic analysis
The transcriptional expression profile data and information of thyroid cancer patients obtained from the TCGA database were plotted using Kaplan-Meier curves with R software.Pvalues and hazard ratios with 95% confidence intervals (CI) were derived using a log-rank test and univariate Cox regression; the data were presented as forest plots.

The prognostic signature model
The relationship between prognostic immune-related gene expression and overall survival (OS) was assessed using LASSO COX analysis.the prognostic risk prediction model for THCA was based on the LASSO risk score calculation formula.patients with THCA were divided into high-risk and low-risk groups.KM curves were plotted to compare OS in the high-risk and lowrisk fractions.receiver operating characteristic survival analysis was performed using the R package SURVIVAL; decision curve analysis was performed using the risk model decision analysis (rmda) package.The relationship between risk score models and tumor immune infiltrating cells was also observed.

Target gene miRNA prediction and correlation analysis
Target gene miRNA prediction was performed using the microRNA target prediction database (miRDB) website (http:// mirdb.org/).Transcriptional profiles and miRNA expression profile data were obtained from the TCGA database.Correlation plots between the two were implemented using the R ggstatsplot package.Spearman's correlation analysis can perform correlations between quantitative variables, which were statistically significant with a p-value less than 0.05.

Cell culture and RT-qPCR
RT-qPCR assays for human normal thyroid cell line Nthyori 3-1 and human THCA cell line ACT-1 were purchased from the American Type Culture Collection.All cell lines were cultured in RPMI-1640 containing 10% fetal bovine serum (FBS, Gibco), 1% penicillin/streptomycin, and maintained in a 37°C, 5% CO 2 incubator.Total RNA was extracted from each group of cells using TRIzol reagent (Invitrogen, Carlsbad, CA) following the manufacturer's guidelines.RNA was synthesized into cDNA using Primescript RT Reagent Kit (Takara, Dalian, China).cDNA amplification and quantification were performed on an Applied Biosystems 7500 instrument using SYBR Green mix (Takara).The relative expression of mRNA was calculated using the 2 -DDCt method, and the experiment was repeated three times independently for each sample.miRNA was used as an internal reference for U6, and glyceraldehyde-3-phosphate dehydrogenase (GAPDH) was used as an internal reference for genes.The PCR primers are shown in Table 1.

Statistical analysis
Cellular tests were performed using the IBM Statistical Package for the Social Sciences (SPSS) 24.0 and GraphPad Prism8.0.1 software to process the data.Normally distributed measures were expressed as (x ± s), and an independent samples t-test was used to compare the two parts.A one-way analysis of variance can comparison between multiple parts, and Tukey's multiple comparisons test was used for post hoc tests.P was a two-sided test, and differences were statistically significant at p < 0.05.

Immunoprognostic analysis of differentially significant genes in thyroid cancer
We obtained transcriptomic data from 510 thyroid cancer patients and 58 paraneoplastic tissues from the TCGA database for differential gene analysis; a total of 676 up-regulated genes and 737 downregulated genes were obtained by |log2FC|>1 with p < 0.05. Figure 1A shows the volcano plot, and Figure 1B shows the heat map.We then obtained 264 identical genes by taking the intersection of the DEGs and immune gene sets (Figure 1C).We performed a prognostic analysis of these 264 genes based on individual gene expression using both one-way Cox and

Immunogenic prognostic model construction
To further explore the clinical value of significant differences in immune-related expression, prognostic features were established for 21 immune-related genes based on LASSO Cox analysis (Figures 2A, B

Significant differences exist in miRNAs regulating the heterogeneous expression of prognostic factors
To understand the heterogeneous expression mechanism of immune prognostic factors in thyroid cancer and adjacent tissues, we introduced the miRNA transcriptome information of thyroid cancer and non-cancer tissues into the GSE103996 dataset of the GEO database, and obtained 58 differential miRNAs through analysis.Among them, 21 were up-regulated and 37 were down-regulated.The heat map is shown in Figure 3B.We then screened for miRNAs that bind to postimmune factors, and their targeting relationships are shown in the Sankey diagram in Figure 3C.

Discussion
Tumor immunity is associated with proto-oncogene activation and tumor immune escape; therefore, understanding tumor-associated immune genes is crucial for tumor immunity research.In recent years, THCA immune studies have focused primarily on developing immune prognostic models (19)(20)(21), and less research has been involved in the regulatory mechanisms of prognostic immune factors.This study mainly analyzed the relationship between immune prognostic factors and miRNAs in THCA.We provided a theoretical basis for further understanding of the immune gene regulation (A-V) is a significantly different miRNA in the GEO database that has a target binding relationship with a prognostic factor and the impact on the prognosis of thyroid cancer.
mechanism of THCA and a reference for assessing the prognosis of THCA by analyzing the correlation between the two.This study analyzed the relationship between immune genes and prognosis using THCAGA and immunoinformatics databases.Twenty-one genes with significant differences in the immune prognosis of THCA were screened, and nine prognostic immune factors were identified by constructing a LASSO risk model (AKAP12, APOC1, TIMP3, ADAMTS9, ANK2, HTRA3, SYNDIG1, ADAMTS5, and DACT1).This risk model predicted OS after five years with high accuracy (AUC = 0.832).
The nine immune genes in the model have biological characteristics and are related to the clinical prognosis of tumor patients.AKAP12 is an A-kinase scaffold protein with a characteristic binding structural domain of the protein kinase A regulatory subunit.AKAP12 deletion is associated with increased cancer susceptibility (22).AKAP1 expression is associated with APOC1 overexpression, promotes tumor progression, and has a poorer prognosis for patient survival (23).TIMP3 is a metalloproteinase that belongs to the TIMP family.It has a high affinity for proteoglycans in the ECM and its broadest range of substrates, including all MMPs, ADAMs (a disintegrin and metalloproteinase), and ADAMTSs (ADAMs with thrombospondin motifs) (24).
TIMP-3 has been shown to have anti-metastatic effects by inhibiting matrix metalloproteinases and members of the ADAM family and down-regulating angiogenesis (25).ADAMTS5 and ADAMTS9 belong to the ADAMTS (disintegrin-like and metalloproteinase with platelet abasicin motif) proteins.They are extracellular zinc metalloproteinases that play an essential role in extracellular matrix assembly and degradation, connective tissue structure, angiogenesis generation, and cell migration (26).
ADAMTSL5 plays a role in maintaining the function of critical oncogenic signaling pathways, suggesting that it may act as a significant regulator of tumorigenicity and drug resistance (27).adamts9 regulates cancer cell growth and metastasis (28) and is associated with patient survival (29).ANK2 has a high mutation frequency in some cancers (30) and silencing of ANK2 expression reduces the growth and invasion of cancer cell type.HTRA3 has been implicated as a tumor suppressor in cancer progression in several cancer types, HTRA3 expression is negatively correlated with adaptive immune cell abundance (T helper cell 17 cells) and positively correlated with innate immune cells (natural killer cells, macrophages, etc.); abundance is positively correlated (31).SYNDIG1 is a prognostic immune factor in diffuse large B-cell lymphoma and breast cancer (32,33).DACT1 belongs to the DACT (Disheveled-associated antagonist of b-catenin) family and is a methylation biomarker for DACT1 in esophageal squamous cell carcinoma (34).
The above data suggest that the nine immune genes in our constructed risk model can exist as prognostic factors in some tumors, reflecting the accuracy of our constructed risk model from the side.
In many studies, the regulatory mechanisms of prognostic factors remain unclear.Based on this starting point, this study analyzed the correlation between immune prognostic factors and miRNAs through TCGA and GEO databases to find the miRNA mechanism of immune prognostic factors in thyroid cancer species.
It has been shown that miRNAs can regulate the target genes involved in the prognosis of tumor patients.For example, in nasopharyngeal carcinoma, the negative correlation between miRNA-19a-3p and PDCD5 expression levels, miRNA-19a-3p targeting to suppress PDCD5 expression, miRNA-19a-3p levels correlated with N classification and clinical stage of nasopharyngeal carcinoma patients, and PDCD5 levels correlated with T classification, pathological grade, and clinical stage.Survival analysis showed that the expression of high levels of miRNA-19a-3p or low levels of PDCD5 had a poorer prognosis in patients with nasopharyngeal carcinoma (35).
Our study identified 16 groups of miRNAs that showed negative correlations with prognostic immune factors, with specific target-binding relationships confirmed in other cancers.For example, TIMP3 was downregulated in cervical cancer and was related to poor prognosis in cervical cancer patients, and miR-21-5p target binding to TIMP3 modulated the development of cervical cancer (36).In gastric cancer, the miR- 21-5p binding relationship with TIMP3 is related to drug resistance in gastric cancer (37).In prostate cancer, the miR-181b-5p-TIMP3 axis regulates prostate cancer proliferation, migration, and invasive ability (38).
This study also found that by constructing a miRNA regulatory network with certain selected candidate factors (39-42), we could provide a new approach to cancer diagnosis and treatment.The miRNA regulatory network of prognostic immune factors was screened.We also validated at the cellular level that miR-15a-5p showed a significant negative correlation with ADAMTS5, ANK2, and SYNDIG1.miR-30c-5p and miR-34b-5p showed a significant negative correlation with ADAMTS9, but no functional gene.However, no further validation of the functional role of the genes was performed, which will be the direction of our later research.
In summary, nine immunoprognostic factors (AKAP12, APOC1, TIMP3, ADAMTS9, ANK2, HTRA3, SYNDIG1, ADAMTS5, and DACT1) were screened as biomolecular markers for predicting THCA.A prognostic assessment model and risk score system were constructed to predict THCA 5-year OS rates with high and low-risk groups using risk scores.The 5year OS rate of THCA was predicted with high accuracy.By constructing a prognostic immune factor miRNA regulatory network, we can provide a theoretical basis for the prognostic immune regulatory mechanism of THCA.

Publisher's note
All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers.Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

1
FIGURE 1 Screening of immune-related prognostic genes in DEGs.(A) Volcano plot, (B) Heat map, (C) Wayne plot, and (D) 21 differential genes with prognostic value.

2
FIGURE 2Construction of prognostic features of immune-related genes in thyroid cancer (A, B) distribution of LASSO coefficients of 21 pivot genes, the adjustment parameter l.min= 0.0065 is obtained, and the vertical black dotted line in B defines the optimal l value; (C) Thyroid cancer Risk score (Riskscore) distribution of patients, survival status and duration of thyroid cancer patients (D) KM survival curves of high and low-risk groups; (E): ROC curves of 1 year, 3 years and 5 years.

3
FIGURE 3Targeted regulation of differentially significant miRNAs that regulate heterogeneous expression of prognostic factors.(A, B) volcano plot of differentially significantly expressed miRNAs in the GSE103996 dataset, heat map; (C) Sankey plot of miRNAs that target predicted binding to prognostic factors.
FIGURE 4 (A-Y) is the expression of significantly different miRNAs in the GEO database that have a targeted binding relationship with the prognostic factors in the TCGA database.ns means no statistical significance, ** means p<0.01, *** means p<0.00.1.

TABLE 1
Sequence of primers.