Effect of the m6ARNA gene on the prognosis of thyroid cancer, immune infiltration, and promising immunotherapy

Article type: Research Article

Keywords: thyroid cancer, m6A, prognosis, risk model, immune infiltration, immunotherapy

Authors: Minqi Xia, Shuo Wang, Yingchun Ye, Yi Tu, Tiantian Huang, Ling Gao

Affiliations: Department of Endocrinology and Metabolism, Renmin Hospital of Wuhan University, Wuhan, China; Department of Breast and Thyroid Surgery, Renmin Hospital of Wuhan University, Wuhan, China

License: Copyright © 2022 Xia, Wang, Ye, Tu, Huang and Gao CC BY 4.0 This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.

Article links: DOI: 10.3389/fimmu.2022.995645  |  PubMed: 36389678  |  PMC: PMC9664221

Relevance: Relevant: mentioned in keywords or abstract

Full text: PDF (10.1 MB)

Introduction

Thyroid cancer is one of the most common tumors in human subjects (ref. 1, ref. 2). The main types of thyroid cancer include papillary carcinoma (PTC), follicular thyroid cancer (FTC), undifferentiated thyroid cancer (ATC) and medullary thyroid cancer (MTC) (ref. 3). The incidence of thyroid cancer has been increasing (ref. 4). At present, the incidence rates of thyroid cancer are 7.4 and 22.0 per 100,000 for males and females in the United States (US), respectively (ref. 5). In addition, thyroid cancer is usually asymptomatic (ref. 6), and at autopsy, 35.6% of the Finnish population had occult thyroid cancer (ref. 7). Currently, most thyroid cancers are curable with conventional treatments such as surgery, radioactive iodide (RAI) therapy,thyroid stimulating hormone (TSH) suppression therapy for local or localized disease (ref. 8) and bioinformatics is also a useful tool to analyze sequencing results and clinical data, construct a prognostic model and seek involved pathways or therapeutic targets, particularly in tumor research (ref. 9–ref. 11). Therefore, using a panel of susceptible genes with RNA/DNA sequencing techniques to identify high-risk patients is advocated for the early diagnosis of high-risk patients and precision medicine (ref. 12, ref. 13).

M6A is a ubiquitous RNA modification in eukaryotes (ref. 14) and is involved in a variety of biological processes, such as embryonic development, apoptosis, spermatogenesis, and circadian rhythms (ref. 15–ref. 17). The modification consists of three processes: catalysis, recognition, and removal, which are completed by m6A methyltransferase “writers”, m6A binding protein “readers” and demethylase “erasers”, respectively (ref. 18–ref. 20). Methyltransferases include METTL3/14, RBM15, RBM15B, WTAP, KIAA1429, and ZC3H13, which are critical in regulating stem cell pluripotency, cell differentiation, and the circadian cycle. RNA-binding proteins, including YTHDC1/2, YTHDF1/2, HNRNPC, ELF3, IGF2BP2 and CBLL1, play an important role in recognizing RNA methylation information and participating in the translation and degradation of downstream RNA. Demethylase consists of FTO and ALKBH5, which mediate RNA demethylation and play a role in energy homeostasis, adipocyte differentiation, and fertility in mice (ref. 21–ref. 24). Dysregulation of m6A has been linked to cancer progression (ref. 25–ref. 30). The effect of m6A-related genes on cancers has been widely studied, such as breast cancer (ref. 31), gastrointestinal cancer (ref. 32), urothelial carcinoma (ref. 33), gastric cancer (ref. 34). In pancreatic cancer, ALKBH5-mediated upregulation of DDIT4-AS1 maintains pancreatic cancer stemness and inhibits chemosensitivity through activation of the mTOR pathway (ref. 35). METTL14 regulates the miR-30c-2-3p/AKT1S1 axis to inhibit gastric cancer progression through mediated mA modifications of circORC5 (ref. 36). Less research has been done on m6A in thyroid cancer, but recent studies have revealed that METTL3 inhibits the progression of papillary thyroid carcinoma (ref. 37).With the development of recent years, in addition to surgical or conventional treatment, cancer immunotherapy has also made amazing progress in thyroid cancer (ref. 38, ref. 39), but the relationship between thyroid cancer and immunity is also unclear.

To further understand the function of m6A-related genes in thyroid cancer, we selected 17 m6A-related genes (ref. 40–ref. 42) for clinical correlation analysis, immune infiltration analysis and immunotherapy analysis. In this study, thyroid cancer patients were divided into 2 clusters using consensus clustering based on the expression of 17 m6A-related genes. Differences in clinical traits between the two cohorts were analyzed. In addition, we developed a risk model by LASSO regression for predicting overall survival (OS) and further explored the relationship between the risk score and biological processes (BPs), cellular components (CCs), molecular function (MF), pathways, outcomes, immune infiltration and immunotherapy.

Materials and methods

Ethics statement

Data for all subjects were obtained from the internet and signed informed consent by default. Analysis of thyroid tissues from this study was carried out under the recommendations of the Regional Ethics Committee guidelines and institutional policies and the Ethics Committee of Renmin Hospital of Wuhan University. Thyroid nodule tissue was obtained from patients with malignant or benign thyroid nodules. All patients signed informed consent forms. Freshly harvested thyroid tissues were used for immunohistochemical staining.

Data processing

RNA expression quantification data (510 cancers vs. 58 normal) and the corresponding clinical data of the THCA cohort were obtained from the TCGA database (https://portal.gdc.cancer.gov/) (ref. 43). The fragments per kilobase of exon model per million mapped fragments (FPKM) values were normalized with the transcripts per million (TPM) method and then converted (log₂+1) in TCGA. The clinical information of all thyroid patients is shown in Supplementary Table S1. The GSE33630 dataset (49 cancers vs. 45 normal) and GSE60542 (33 cancers vs. 35 normal) cohort in our study were downloaded from the GEO database (http://www.ncbi.nlm.nih.gov/geo) (ref. 44).The protein-protein interaction (PPI) diagram of 17 m6A-related genes was constructed in the Search Tool for the Retrieval of Interacting Genes/Proteins (STRING) network (https://string-db.org/) (ref. 45). PPI genes were analyzed by Cytoscape software. The genetic database used for GSEA data analysis is “MSigdb” (ref. 46).

Differential expression analysis

Differentially expressed genes (DEGs) in normal and tumor thyroid samples were screened by the “limma” package. We extracted the expression matrix of the 17 m6A-related genes, including METTL3/14, RBM15, RBM15B, WTAP, KIAA1429, ZC3H13, YTHDC1/2, YTHDF1/2, HNRNPC, ELF3, IGF2BP2, CBLL1, FTO, and ALKBH5. The R packages “heat-map” and “vioplot” were used to show the differential expression of 17 m6A-related genes in normal and tumor samples. The correlation between 17 m6A-related genes was constructed by the R package “correlation”.

Protein-protein interaction

We uploaded 17 m6A-related genes to the Search Tool for the Retrieval of Interacting Genes/Proteins (STRING) network to construct a protein-protein interaction (PPI) network to show the association between m6A-related genes.

The relationship between clusters and clinical traits after clustering analysis

The “ConsensusClusterPlus” package was used to cluster data according to the expression of 17 m6A-related genes in 568 thyroid cancer samples. The consensus matrix and cumulative distribution function (CDF) were used to calculate the optimal number of clusters. The difference in clinical traits between the two clusters was illustrated by the R packages “survival” and “heat-map”.

Risk model construction and clinical correlation analysis

Univariate Cox analysis was used to evaluate the prognostic value of 17 m6A-related genes with which HR<1 or >1 was regarded as a protective gene or risk gene, respectively. Whether there was an association between genes and prognosis was determined by the P value (p<0.05). To avoid overfitting, least absolute shrinkage and selection operator (LASSO) regression analysis was performed using the “glmnet” package. The risk score of each patient was calculated by the following formula:

Coef(i) and x(i) are regression coefficients and gene expression levels, respectively.

Based on their risk score being below or above the median risk score, 449 THCA samples were divided into two subgroups: the high-risk group and the low-risk group. The distribution of the risk score, survival curve, heatmap of the 3 genes in the model, and receiver operating characteristic (ROC) curve were analyzed. The relationship between the risk score and clinical traits, such as fustat, age, gender, stage, and tumor, nodes and metastases (TNM) stage, was also assessed. A nomogram was created using the package “rms”, and the ROC curve was used to prove the validity of the model. We further combined clinical traits with univariate and multivariate Cox regression analyses using the “survival” package. Analysis variables included age, gender, stage, TNM stage, and risk score. Finally, Cox regression analysis was used to analyze the relationship between each clinical trait and survival probability. The Wilcoxon test was used to calculate the relationship between each clinical trait and risk scores.

Differential gene enrichment analysis between the high- and low-risk groups

Patients with THCA were divided into high- and low-risk groups based on the median risk score. |log₂FC|≥1和false discovery rate (FDR)<0.05 as a screening criterion for differential genes between the two groups (high- and low-risk groups). Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) enrichment analyses were performed using the “cluster profile” package based on the selected risk DEGs. GSEA was performed using GSEA version 4.2.3 based on data from the high- and low-risk groups. Set the minimum gene set as 15 and the maximum gene set as 500, and resampling a thousand times, a P value of< 0.05 and an FDR of< 0.25 were considered statistically significant.

Immune infiltration analysis and immunotherapy

We used the R package “limma” of the “normalizeBetweenArrays” function to reduce batch effects that may exist between or within the two cohorts to merge the gene information in GSE33630 and GSE60542.The “ssGSEA” package was used to calculate the immune infiltrating cell score for the TCGA and GEO cohorts, and the score was used to compare immune function and immunological pathways between the high- and low-risk groups, and the relationship between risk score and immune function and immunological pathway was also discussed. To assess the potential impact of the risk score on the immune checkpoint response, we used the “limma” package to calculate the difference between the high- and low-risk groups, |log₂FC|≥1 and false discovery rate (FDR)<0.05 as significance. We used the data from Genomics of Drug Sensitivity in Cancer (GDSC) (https://www.cancerrxgene.org/) database (ref. 47, ref. 48)and The Cancer Therapeutics Response Portal (CTRP)(https://portals.broadinstitute.org/ctrp/) database (ref. 49) to obtain the relationship between the expression levels of three genes (RBM15, KIAA1429, FTO) and immunotherapeutic drugs. In addition, drug sensitivity and CCLE expression data were obtained from the PRISM Repurposing dataset (https://depmap.org/repurposing) (ref. 50) to assess the relationship between drugs and risk scores.

Immunohistochemistry

Immunohistochemistry (IHC) was performed as previously described (ref. 51). Rabbit antihuman antibodies (FTO antibody (No. DF8421, Affinity Biosciences, Cincinnati, OH, USA), RBM15 antibody (No. DF12061, Affinity Biosciences USA), KIAA1429 antibody (No. 25712-1-AP, Proteintech, Rosemont, IL, USA), were used at 1:500 dilutions. The secondary antibody was peroxidase-labelled antibody (rabbit IgG (H+L) KPL, Baltimore, MD, USA).

Results

Research technology road map

The flow chart of this study is shown in f1. By using expression data from the TCGA-THCA cohort, 16 m6A-related genes were selected that differed between normal and tumour samples (|Log₂FC|>1, FDR<0.05). Risk model based on the three m6A-related genes were then constructed using and LASSO regression analysis, and risk scores were also calculated. Patients were classified into high- and low-risk groups based on the median value of their risk scores in the TCGA-THCA cohort and were further used for subsequent clinical, immune infiltration, mutation, enrichment analysis and immunotherapy analysis. GEO cohort was also used to validation the immune function with risk score. Finally, in vitro validation was performed using thyroid cancer and normal tissues.

Expression of 17 m6A-related genes in the TCGA-THCA cohort

First, the TCGA-THCA cohort in the TCGA database was analyzed to compare the expression of 17 m6A-related genes in normal and tumor tissues (|Log₂FC|>1, FDR<0.05). As shown in the heatmap and violin diagram, there were significant differences in 16 genes between normal and tumor tissues, including significantly upregulated expression of HNRNPC, IGFBP2, ELF3, and RBM15B and downregulated expression of ZC3H13, FTO, KIAA1429, WTAP, YTHDC1, YTHDC2, ALKBH5, RBM15, METTL3, METTL14, YTHDC1, and CBLL1 in the TCGA-THCA cohort (f2). Spearman correlation analysis further supported the relationship between m6A-related genes (P<0.05), such as KIAA1429 and FTO, KIAA1429 and METTL14, ZC3H13 and METTL14, which were significantly positively correlated; ALKBH5 was negatively correlated with IGFBP1; and CBLL1 was negatively correlated with ELF3. Taking KIAA1429 and METTL14 as an example, the correlation coefficient=0.78, P<0.0001, indicating that there is a strong correlation between them (f2). A protein-protein interaction (PPI) network of the 17 m6A-related genes was constructed (f2). Next, we used the “Cytohubba” software to calculate the 17 genes and obtained 11 hub genes with the highest scores after calculation: ALKBH5, FTO, METTL14, KIAA1429, RBM15, RBM15B, METTL3, ZC3H13, YTHDC1, YTHDF1, and YTHDC (f2).

Consensus clustering of 17 m6A-related genes yielded two clusters

To further understand the overall role of m6A-related genes in the TCGA-THCA cohort, consensus clustering analysis was performed on 568 TCGA-THCA samples based on the expression profiles. Duplicate samples were removed, and similar samples were grouped into one class. To determine the optimal number of clusters, we varied the total number of clusters from 2 to 9 and examined the cumulative density function (CDF) curve of the consensus matrix (f3). Finally, K=2 is taken as the optimal cluster number according to the consensus matrix (f3). These subgroups can also be distinguished in principal component analysis (PCA) (f3). The Kaplan-Meier curves showed a significant difference in survival probability between the two clusters, with Cluster 1 having a longer OS (P<0.05) (f3). To determine whether there was a significant difference in clinical traits between the two clusters, a clinically relevant heatmap was produced. Based on the heatmap, we observed significant differences in stage (P<0.01), T stage (P<0.01), and N stage (P<0.001) between the two subgroups (f3).

Construction and value of the risk model

Univariate Cox regression analysis produced four genes (IGFBP2, KIAA1429, RBM15, and FTO). The forest map showed that KIAA1429, RBM15, and FTO were risk factors (HR>1), and IGFBP2 was a protective factor (HR<1). IGFBP2 was not included in the LASSO regression analysis because it was not a hub gene in 17 m6A-related genes (f4). Three risk factors (KIAA1429, RBM15, and FTO) were selected, and LASSO regression analysis was performed. These 3 genes (KIAA1429, RBM15, and FTO) were retained according to the minimum partial likelihood bias (f4). For the TCGA-THCA cohort, the risk score was calculated by the following formula:

The 449 samples were scored according to a risk scoring formula and were classified as high-risk and low-risk by the median of the risk score (there were 224 high-risk samples and 225 low-risk samples, with the middle 2 values 2.209 and 2.211, respectively). We also described the sample distribution for the three m6A-related genes, risk assessment subgroups, and survival conditions (f4). The accuracy of the risk score was evaluated by calculating the area under the ROC curve (AUC=0.744) (f4). There were significant differences in survival probability between the high- and low-risk groups (P=1.784e-02), and the overall survival of patients with thyroid cancer was poor in the high-risk group (f4). We constructed a prognostic nomogram using risk score and clinical traits such as survival, age, gender, stage, T stage, N stage and M stage. (f5). The c-index of the model population was 0.959 95% confidence interval (CI) (0.939-0.987), P value=1.091e-223. The nomogram can effectively predict the overall survival time of patients with THCA. ROC curve analysis showed that the AUC was 0.538 at 1 year, 0.791 at 5 years, and 0.749 at 10 years (f5). Then, univariate Cox regression analysis showed that age, stage, T stage and risk score were correlated with the prognosis of patients (P<0.05). Multivariate Cox regression analysis showed that age and risk score could be independent prognostic indicators (P<0.05) (f5).

Clinical subgroup analysis

We extracted the clinical data of age, gender, stage, and TNM stage to further explore the relationship between clinical traits and survival probability. Then, we eliminated the samples for which information was missing and carried out Cox regression analysis on the remaining samples to obtain the clinical traits and survival probability as well as the relationship between the clinical traits and risk score. According to our findings, age>65 had a significantly lower chance of survival than age ≤ 65 (P<0.001) (f6). Stage I-II, T1-T2 and M0 patients had a significantly higher chance of surviving than stage III-IV, T3-T4 and M1 patients (P<0.05). In contrast, neither gender nor N stage seemed to significantly affect survival probability (f6). Regarding the differences in clinical traits between the high- and low-risk groups, the high- and low-risk groups’ survival probabilities were significantly different (P<0.05) for patients who were over 65 and in Stages III-IV (Supplementary Table S2). Additionally, the risk scores between the T1-T2 and T3-T4 groups differed significantly (P<0.05). However, other clinical traits (age, gender, stage, N, and M) did not appear to be associated with a difference in risk scores (Supplementary Table S3).

Functional analysis based on the risk model

To further explore the functional difference between the high- and low-risk groups, we used the R package “limma” to screen DEGs, and the screening standard was set to FDR<0.05 and |log₂FC|≥1. A total of 76 DEGs were screened between the high-risk and low-risk groups in the TCGA cohort. Among these genes, 49 genes in risk DEGs were upregulated, and the remaining 27 genes were downregulated. Next, we performed GO, KEGG, and GSEA analyses based on 76 risk DEGs. GO and KEGG analyses showed that risk DEGs were related mainly to proteolysis, lipid metabolism, and immune response (f7). GSEA showed that risk DEGs were related mainly to the cancer pathway (f7).

Immune infiltration and immunotherapy

We further compared the enrichment fractions of immune cells and the activity of immune-related pathways between the high- and low-risk groups in the TCGA-THCA cohort. A comparison of the various infiltrative immune cells in tumors showed that the degree of DC, iDC and NK-cell infiltration was higher in the low-risk group than in the high-risk group, while other immune cells did not differ between the two groups (f8). In the TCGA cohort, the activity of the Type_II_IFN pathway was significantly higher than the activity of the low-risk group, while the activity of the APC_costimulation pathway was higher in the low-risk group. The other 11 immune pathways showed little difference between the high-risk and low-risk groups (f8). Next, we further investigated the relationship between risk scores and individual immune cells and validated it using GEO data. TCGA data showed that risk score was significantly positively correlated with Type-II-IFN-Response pathway, and this conclusion was also confirmed in GEO database (f8). We used volcano plots to show the differences in immune checkpoints between high- and low-risk groups, P value of < 0.05 and |log₂FC|>1 were considered statistically significant (f8). Furthermore, we mapped the relationship between the risk groups and immune checkpoints (f8). Using the GDSC and CTRP databases, we summarized the correlation between gene expression and drug sensitivity across cancers, and the 30 drugs with the strongest correlation are listed (f9). Furthermore, we analyzed the association between risk scores and responsiveness to nontumor drugs through the PRISM database (f9). The results showed that 20 drugs, including PD-168393, ibrutinib, lapatinib, AS703026, trametinib, PD-0325901, cobimetinib, binimetinib, RO4987655, TAK733, SCH900776, inosine, NSC23766, sulforaphane, BMS3445541, diphenhydramine, 2,3-DCPE, niraparib, mercaptopurine and tacedinaline, had correlations with risk scores.

Expression of RBM15, KIAA1429, and FTO in benign thyroid nodule tissue vs thyroid cancer tissues

Immunohistochemistry result revealed that the staining intensity of 3 genes (RBM15, FTO, and KIAA1429) were higher in thyroid cancer tissue than in benign thyroid nodule tissue in our specimens which is consistent with bioinformatic analysis (f10).

Discussion

In this study, based on data from the TCGA-THCA cohort, we obtained DEGs between normal and tumor tissues and identified 3 m6A-related genes (KIAA1429, RBM15, and FTO) to construct a risk model and verified its credibility. Based on the risk model, the differences in enrichment pathways, immune function, and immunotherapy between the high- and low-risk groups were compared. Small molecule drugs can be designed for intervention based on the risk score and immune checkpoints, which has the potential to improve the treatment of chemotherapy, radiotherapy, and even immunotherapy. The risk model also provides a new possible avenue for exploring the treatment of thyroid cancer based on the risk score.

We analyzed the expression of 17 m6A-related genes in thyroid cancer and adjacent tissues and found that except for YTHDF2, the other 16 genes were differentially expressed in thyroid cancer. In particular, HNRNPC, RBM15B, IGFBP2, and ELF3 were upregulated, while the other 12 genes were downregulated in thyroid cancer tissues. In addition, there were multiple correlations between m6A-related genes. For example, KIAA1429 and FTO, KIAA1429 and METTL14, and ZC3H13 and METTL14 were significantly positively correlated. ALKBH5 was negatively correlated with IGFBP1 and CBLL1 with ELF3. The correlation between some of these genes has been reported, which is consistent with the conclusions of this paper (ref. 52–ref. 55).

Based on the expression of 17 m6A-related genes in TCGA-THCA cohort, we identified two distinct molecular clusters. Compared with Cluster 1, Cluster 2 had worse survival probability. The results showed that there were significant differences between the two groups in the clustering analysis. Therefore, we further built the risk model. To further select genes related to prognosis and establish a risk model, we used univariate Cox analysis to select 3 genes with unfavorable prognosis: RBM15, KIAA1429, and FTO. LASSO regression analysis was further used to calculate the risk scores of each clinical dataset, and the data were divided into low- and high-risk groups based on the median risk scores of all samples. Three genes for constructing the risk model have been shown to be associated with cancer. A study on breast cancer showed that KIAA1429 can be a carcinogenic factor of breast cancer by regulating CDK1 independently of m6A (ref. 56). RBM15 is an important factor in X chromosome silencing and is expressed in breast tissue (ref. 40, ref. 57–ref. 59). Studies have also shown that estrogen can independently increase the expression of the RBM15 gene and is associated with the occurrence of breast diseases (ref. 59). The FTO gene is highly expressed in breast cancer (ref. 60), acute myeloid leukaemia (AML) (ref. 27), and glioblastoma (ref. 61) and promotes the progression of these cancers, and related drugs have been developed. However, FTO showed low expression in bladder cancer (ref. 62). The carcinogenic role of FTO seems to be controversial. In addition to being associated with cancer, FTO is also associated with immunity. KIAA1429 has been shown to affect a variety of cancers through different pathways, including promoting the progression of hepatocellular cancer through the posttranscriptional modification of m6A relying on GATA3 and regulating cell proliferation in gastric cancer by directly targeting c-jun mRNA (ref. 62, ref. 63). Immunosuppressants have been reported to be widely used in the clinical treatment of malignant tumors, including thyroid cancer (ref. 64). RBM15 expression was positively correlated with immune-infiltrating cells in renal clear cell carcinoma (KIRC), brain low-grade glioma (LGG), and pancreatic adenocarcinoma (PAAD). In addition, RBM15 expression was strongly correlated with immune checkpoint markers in PAAD (ref. 65). In addition, some studies have found that downregulation of FTO in the chorionic villi destroys immune tolerance and angiogenesis at the maternal-fetal interface, leading to abnormal methylation and oxidative stress and ultimately to spontaneous abortion (ref. 66).

We found that the survival rate of the high-risk group was significantly lower than that of the low-risk group. Next, we constructed a nomogram to further improve the predictive power, which showed that the predictive nomogram could be applied to the TCGA-THCA cohort. As the total clinical score increased, the survival probability of the patient gradually decreased. The C-index of this nomogram reached 0.959, indicating the reliability of the nomogram. We also analyzed the differences in survival probability for each clinical subgroup, and there were significant differences in survival among age, stage, and T-stage subgroups. In order to explore the impact of multiple clinical traits and risk score on the prognosis of patients, we further performed a multivariate Cox analysis, which showed that age and risk score were indeed independent risk factors for the prognosis of thyroid cancer patients. The risk score, as an independent risk factor for prognosis, has also been studied in other diseases (ref. 40).In conclusion, this risk model has higher diagnostic power than previous cluster analyses and is more helpful in assessing patient prognosis.

In addition, according to the GO, KEGG and GSEA analyses of the DEGs in the high- and low-risk groups, the analysis results showed that the functions of the risk DEGs were related mainly to proteolysis, immune response, and cancer pathway. In the process of tumor migration, proteolytic enzymes can degrade the basement membrane, make cancer cells through the basement membrane, and migrate elsewhere, in which proteolytic enzymes play an important role in the process of tumor invasion and metastasis. For example, some proteases in thyroid cancer are elevated in tumors (ref. 67):transmembrane protease serine 4 promotes thyroid cancer proliferation through cAMP response element-binding protein (CREB) phosphorylation (ref. 68), and the HIV protease inhibitor nefinavir induces apoptosis of thyroid medulla cancer cells by downregulating RET signalling (ref. 69). These results suggest that the formation mechanism of thyroid cancer in the high-risk group was related to the action of proteolytic enzymes. In our study, the GO results showed that differences in immune function between the high-risk and low-risk groups were manifested mainly in monocyte differentiation, lymphocyte differentiation and B-cell activation. GSEA showed that there were significant differences between the high-risk and low-risk groups in the ERBB pathway and small cell lung cancer pathway, which overlapped with our results of immune function analysis and further demonstrated that our prognostic model was inseparable from cancer and immunity.

Currently, the treatment of thyroid cancer is limited mainly to surgical resection (ref. 70, ref. 71), which can improve the survival probability of patients. However, the treatment of thyroid cancer has expanded, and great progress has been made with anticancer drugs. Immunotherapy has become an important part of thyroid cancer treatment (ref. 72–ref. 77). Therefore, we performed immune analysis for the high- and low-risk groups. In this study, compared with the low-risk group, the immune infiltration levels of dendritic cells, iDCs and NK cells in the high-risk group were lower, but Type II IFN pathway activity was higher. Moreover, we can see that Type II IFN pathway activity are significantly positively correlated with risk score in both TCGA and GEO cohorts. The decrease in dendritic cells, iDCs and NK cells, which are key cells in the immune response and suppress the development of primary tumors and metastases (ref. 78–ref. 82), may be related to the lower survival rate in the high-risk group. Type II IFN has a rejection effect on highly immunogenic tumors (ref. 83). IFN-γ, the only type II IFN, binds to the type II IFN receptor and activates Jak1 and Jak2, which then further activates STAT1. The activated STAT1 homodimer binds to the IFN-γ activation site in the promoter of some ISGs to regulate immune function (ref. 84). However, it has also been shown that long-term exposure to IFN-γ can lead to immune escape due to cell desensitization and immune editing (ref. 85, ref. 86). Its expression was elevated in the high-risk group, suggesting that Type-II IFN is involved in the antithyroid cancer effect in the high-risk group of thyroid cancer. The increased expression of Type-II IFN in high-risk group indicates that it may be involved in the anti-thyroid cancer effect, but it cannot exclude that it is involved in the immune escape of thyroid cancer.

Compared with the low-risk group, CD44 and NPR1 were significantly increased in the high-risk group. CD44 is encoded by 20 exons that are alternatively spliced to generate the CD44 standards (CD44s) and the CD44 variants (CD44v) (ref. 87). Recent studies have shown that abnormal levels and variant forms of CD44 are expressed in a variety of tumor types (ref. 88–ref. 95). High expression of CD44 is thought to be associated with poor prognosis in cancer, and in addition to this, CD44 v6 and v3 have been reported to be associated with cancer metastasis. Inhibition of CD44 has a cancer-inhibiting effect. It has been shown that inhibition of CD44 inhibits the development of colon tumors (ref. 96) in mice and suppresses the proliferation and metastasis of liver and ovarian cancer cells (ref. 97, ref. 98). NRP1 is also involved in the cell proliferation, migration, and invasion associated with cancer progression (ref. 99). In colorectal cancer, inhibition of NRP1 was able to inhibit metastasis of colorectal cancer cells (ref. 100). Moreover, NRP1/mdm2-targeted d-peptide supramolecular nanodrugs have been shown to be highly effective and less toxic for the treatment of hepatocellular carcinoma, with strong anti-cancer activity against SK-Hep-1 cells in vitro and in vivo, without significant host toxicity, making them a promising treatment for hepatocellular carcinoma (ref. 101). In summary, the CD44 and NRP1 inhibitors to treat cancer patients with high-risk scores are promising. From the perspective of the relationship between immunotherapy drugs and risk score, the drugs for patients with high risk score mainly focus on the class of teninib drugs, and these drugs may provide new insights into the treatment of thyroid cancer.

Finally, we verified the expression of 3 genes in thyroid cancer and thyroid nodule tissues by immunohistochemistry. Compared with the thyroid nodule tissues, the expression of 3 genes was increased in thyroid cancer tissues, indicating that the 3 genes constructed as models were also associated with the development of thyroid cancer.

Conclusion

Based on bioinformatics and in vitro experiments, we determined that there were differences in the expression of the m6A-related gene in thyroid cancer. The risk model was able to predict the survival probability of patients with thyroid cancer and confirmed that the risk score was an independent risk factor for thyroid cancer. In addition, enrichment analysis of the high-low risk group confirmed that there were significant differences in proteolysis, immune response, and cancer formation in the high/low-risk group. Further analysis of immune function showed that there were differences in DCs, iDCs, NK cells, and Type II IFN response between the two groups. In addition, immune checkpoint analysis yielded two immune (CD44 and NRP1) checkpoints with therapeutic potential. Finally, the immunodrug analysis identified 20 drugs that were associated with risk scores.

There is no doubt that there are some limitations to the study. First, all analyses were based on retrospective data from a public database. In addition, we have only performed a simple study on the relationship between BP, CC, MF, pathway, outcome, immune function and immunological pathways, immune checkpoint, and THCA and the actual effects of proteolytic, iDCs, NK cells, and Type-II IFN. However, this study systematically analyzed the expression differences of m6A-related genes in thyroid cancer and constructed a risk model that could evaluate the survival, prognosis, and immunotherapy of patients.

Data availability statement

The raw data supporting the conclusions of this article will be made available by the authors, without undue reservation.

Ethics statement

The studies involving human participants were reviewed and approved by Ethics Committee of Renmin Hospital of Wuhan University. The patients/participants provided their written informed consent to participate in this study.

Author contributions

MX and TH performed the data analysis and wrote the manuscript. SW ,YY,LG participated in the study design, figures,writing and revised the manuscript. YT provided tissue samples.All authors read and approved the final manuscript.

Funding

LG supported from the National Natural Science Foundation of China (No. 81571376), the Fundamental Research Funds for the Central Universities, #2042020kf1079 and China Young Scientific Talent Research Fund for diabetes 2017.

Conflict of interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

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.

References

1. BR Haugen, EK Alexander, KC Bible, GM Doherty, SJ Mandel, YE Nikiforov. 2015 American Thyroid association management guidelines for adult patients with thyroid nodules and differentiated thyroid cancer: The American thyroid association guidelines task force on thyroid nodules and differentiated thyroid cancer.. Thyroid (, 2016. [DOI | PubMed]
2. Y Zheng, P Nie, D Peng, Z He, M Liu, Y Xie. m6AVar: a database of functional variants involved in m6A modification.. Nucleic Acids Res (, 2018. [DOI]
3. M Xing. Molecular pathogenesis and mechanisms of thyroid cancer.. Nat Rev Cancer (, 2013. [DOI]
4. S Lohia, M Hanson, RM Tuttle, LGT Morris. Active surveillance for patients with very low-risk thyroid cancer.. Laryngoscope Investig Otolaryngol (, 2020. [DOI]
5. J Gudmundsson, G Thorleifsson, JK Sigurdsson, L Stefansdottir, JG Jonasson, SA Gudjonsson. A genome-wide association study yields five novel thyroid cancer risk loci.. Nat Commun (, 2017. [DOI | PubMed]
6. DS Cooper, GM Doherty, BR Haugen, RT Kloos, SJ Mandel, EL Mazzaferri. Revised American thyroid association management guidelines for patients with thyroid nodules and differentiated thyroid cancer.. Thyroid (, 2009. [DOI]
7. HR Harach, KO Franssila, VM Wasenius. Occult papillary carcinoma of the thyroid. a "normal" finding in finland. a systematic autopsy study.. Cancer (, 1985. [DOI]
8. IM Min, E Shevlin, Y Vedvyas, M Zaman, B Wyrwas, T Scognamiglio. CAR T therapy targeting ICAM-1 eliminates advanced human thyroid tumors.. Clin Cancer Res (, 2017. [DOI]
9. A Chang, NH Chakiryan, D Du, PA Stewart, Y Zhang, Y Tian. Proteogenomic, epigenetic, and clinical implications of recurrent aberrant splice variants in clear cell renal cell carcinoma.. Eur Urol (, 2022. [DOI]
10. H Schinke, E Shi, Z Lin, T Quadt, G Kranz, J Zhou. A transcriptomic map of EGFR-induced epithelial-to-mesenchymal transition identifies prognostic and therapeutic targets for head and neck cancer.. Mol Cancer (, 2022. [DOI | PubMed]
11. FA Buttner, S Winter, V Stuhler, S Rausch, J Hennenlotter, S Fussel. A novel molecular signature identifies mixed subtypes in renal cell carcinoma with poor prognosis and independent response to immunotherapy.. Genome Med (, 2022. [DOI | PubMed]
12. I De Vlaminck, L Martin, M Kertesz, K Patel, M Kowarsky, C Strehl. Noninvasive monitoring of infection and rejection after lung transplantation.. Proc Natl Acad Sci U.S.A. (, 2015. [DOI]
13. I De Vlaminck, HA Valantine, TM Snyder, C Strehl, G Cohen, H Luikart. Circulating cell-free DNA enables noninvasive diagnosis of heart transplant rejection.. Sci Transl Med (, 2014. [DOI]
14. R Desrosiers, K Friderici, F Rottman. Identification of methylated nucleosides in messenger RNA from novikoff hepatoma cells.. Proc Natl Acad Sci U.S.A. (, 1974. [DOI]
15. M Frye, BT Harada, M Behm, C He. RNA Modifications modulate gene expression during development.. Science (, 2018. [DOI]
16. IA Roundtree, ME Evans, T Pan, C He. Dynamic RNA modifications in gene expression regulation.. Cell (, 2017. [DOI]
17. BS Zhao, X Wang, AV Beadell, Z Lu, H Shi, A Kuuspalu. m(6)A-dependent maternal mRNA clearance facilitates zebrafish maternal-to-zygotic transition.. Nature (, 2017. [DOI]
18. PJ Batista, B Molinie, J Wang, K Qu, J Zhang, L Li. m(6)A RNA modification controls cell fate transition in mammalian embryonic stem cells.. Cell Stem Cell (, 2014. [DOI]
19. JM Fustin, M Doi, Y Yamaguchi, H Hida, S Nishimura, M Yoshida. RNA-Methylation-dependent RNA processing controls the speed of the circadian clock.. Cell (, 2013. [DOI | PubMed]
20. S Geula, S Moshitch-Moshkovitz, D Dominissini, AA Mansour, N Kol, M Salmon-Divon. Stem cells. m6A mRNA methylation facilitates resolution of naive pluripotency toward differentiation.. Science (, 2015. [DOI]
21. C Dina, D Meyre, S Gallina, E Durand, A Korner, P Jacobson. Variation in FTO contributes to childhood obesity and severe adult obesity.. Nat Genet (, 2007. [DOI]
22. R Do, SD Bailey, K Desbiens, A Belisle, A Montpetit, C Bouchard. Genetic variants of FTO influence adiposity, insulin sensitivity, leptin levels, and resting metabolic rate in the Quebec family study.. Diabetes (, 2008. [DOI]
23. TM Frayling, NJ Timpson, MN Weedon, E Zeggini, RM Freathy, CM Lindgren. A common variant in the FTO gene is associated with body mass index and predisposes to childhood and adult obesity.. Science (, 2007. [DOI]
24. G Zheng, JA Dahl, Y Niu, P Fedorcsak, CM Huang, CJ Li. ALKBH5 is a mammalian RNA demethylase that impacts RNA metabolism and mouse fertility.. Mol Cell (, 2013. [DOI | PubMed]
25. I Barbieri, K Tzelepis, L Pandolfini, J Shi, G Millan-Zambrano, SC Robson. Promoter-bound METTL3 maintains myeloid leukaemia by m(6)A-dependent translation control.. Nature (, 2017. [DOI]
26. X Deng, R Su, H Weng, H Huang, Z Li, J Chen. RNA N(6)-methyladenosine modification in cancers: current status and perspectives.. Cell Res (, 2018. [DOI]
27. Z Li, H Weng, R Su, X Weng, Z Zuo, C Li. FTO plays an oncogenic role in acute myeloid leukemia as a N(6)-methyladenosine RNA demethylase.. Cancer Cell (, 2017. [DOI]
28. J Liu, MA Eckert, BT Harada, SM Liu, Z Lu, K Yu. m(6)A mRNA methylation regulates AKT activity to promote the proliferation and tumorigenicity of endometrial cancer.. Nat Cell Biol (, 2018. [DOI]
29. R Su, L Dong, C Li, S Nachtergaele, M Wunderlich, Y Qing. R-2HG exhibits anti-tumor activity by targeting FTO/m(6)A/MYC/CEBPA signaling.. Cell (, 2018. [DOI | PubMed]
30. LP Vu, BF Pickering, Y Cheng, S Zaccara, D Nguyen, G Minuesa. The N(6)-methyladenosine (m(6)A)-forming enzyme METTL3 controls myeloid differentiation of normal hematopoietic and leukemia cells.. Nat Med (, 2017. [DOI]
31. H Chen, Y Yu, M Yang, H Huang, S Ma, J Hu. YTHDF1 promotes breast cancer progression by facilitating FOXM1 translation in an m6A-dependent manner.. Cell Biosci (, 2022. [DOI | PubMed]
32. Q Zhao, Y Zhao, W Hu, Y Zhang, X Wu, J Lu. m(6)A RNA modification modulates PI3K/Akt/mTOR signal pathway in gastrointestinal cancer.. Theranostics (, 2020. [DOI]
33. J Kong, S Lu, L Zhang, Y Yao, J Zhang, Z Shen. m6A methylation regulators as predictors for treatment of advanced urothelial carcinoma with anti-PDL1 agent.. Front Immunol (, 2022. [DOI | PubMed]
34. C Yuan, J Zhang, C Deng, Y Xia, B Li, S Meng. Crosstalk of histone and RNA modifications identified a stromal-activated subtype with poor survival and resistance to immunotherapy in gastric cancer.. Front Pharmacol (, 2022. [DOI | PubMed]
35. Y Zhang, X Liu, Y Wang, S Lai, Z Wang, Y Yang. The m(6)A demethylase ALKBH5-mediated upregulation of DDIT4-AS1 maintains pancreatic cancer stemness and suppresses chemosensitivity by activating the mTOR pathway.. Mol Cancer (, 2022. [DOI | PubMed]
36. HN Fan, ZY Chen, XY Chen, M Chen, YC Yi, JS Zhu. METTL14-mediated m(6)A modification of circORC5 suppresses gastric cancer progression by regulating miR-30c-2-3p/AKT1S1 axis.. Mol Cancer (, 2022. [DOI | PubMed]
37. Y Zhu, X Peng, Q Zhou, L Tan, C Zhang, S Lin. METTL3-mediated m6A modification of STEAP2 mRNA inhibits papillary thyroid cancer progress by blocking the hedgehog signaling pathway and epithelial-to-mesenchymal transition.. Cell Death Dis (, 2022. [DOI | PubMed]
38. Z Lu, S Bai, Y Jiang, S Wu, D Xu, J Zhang. Amplifying dendritic cell activation by bioinspired nanometal organic frameworks for synergistic sonoimmunotherapy.. Small (, 2022. [DOI]
39. W Xie, Y Zeng, L Hu, J Hao, Y Chen, X Yun. Based on different immune responses under the glucose metabolizing type of papillary thyroid cancer and the response to anti-PD-1 therapy.. Front Immunol (, 2022. [DOI | PubMed]
40. ZH Yu, ST Feng, D Zhang, XC Cao, Y Yu, X Wang. The functions and prognostic values of m6A RNA methylation regulators in thyroid carcinoma.. Cancer Cell Int (, 2021. [DOI | PubMed]
41. ZY Feng, T Wang, X Su, S Guo. Identification of the m(6)A RNA methylation regulators WTAP as a novel prognostic biomarker and genomic alterations in cutaneous melanoma.. Front Mol Biosci (, 2021. [DOI | PubMed]
42. W Tan, S Liu, Z Deng, F Dai, M Yuan, W Hu. Gene signature of m6A-related targets to predict prognosis and immunotherapy response in ovarian cancer.. J Cancer Res Clin Oncol (, 2022. [DOI]
43. Comprehensive genomic characterization defines human glioblastoma genes and core pathways.. Nature (, 2008. [DOI]
44. D Sia, Y Hoshida, A Villanueva, S Roayaie, J Ferrer, B Tabak. Integrative molecular analysis of intrahepatic cholangiocarcinoma reveals 2 classes that have different outcomes.. Gastroenterology (, 2013. [DOI]
45. D Szklarczyk, AL Gable, D Lyon, A Junge, S Wyder, J Huerta-Cepas. STRING v11: protein-protein association networks with increased coverage, supporting functional discovery in genome-wide experimental datasets.. Nucleic Acids Res (, 2019. [DOI]
46. A Subramanian, P Tamayo, VK Mootha, S Mukherjee, BL Ebert, MA Gillette. Gene set enrichment analysis: a knowledge-based approach for interpreting genome-wide expression profiles.. Proc Natl Acad Sci U.S.A. (, 2005. [DOI]
47. F Iorio, TA Knijnenburg, DJ Vis, GR Bignell, MP Menden, M Schubert. A landscape of pharmacogenomic interactions in cancer.. Cell (, 2016. [DOI]
48. W Yang, J Soares, P Greninger, EJ Edelman, H Lightfoot, S Forbes. Genomics of drug sensitivity in cancer (GDSC): a resource for therapeutic biomarker discovery in cancer cells.. Nucleic Acids Res (, 2013. [DOI]
49. Y Zou, MJ Palte, AA Deik, H Li, JK Eaton, W Wang. A GPX4-dependent cancer cell state underlies the clear-cell morphology and confers sensitivity to ferroptosis.. Nat Commun (, 2019. [DOI | PubMed]
50. SM Corsello, RT Nagari, RD Spangler, J Rossen, M Kocak, JG Bryan. Discovering the anti-cancer potential of non-oncology drugs by systematic viability profiling.. Nat Cancer (, 2020. [DOI]
51. H Yang, S Wang, Y Ye, M Xie, Y Li, H Jin. GLP-1 preserves β cell function via improvement on islet insulin signaling in high fat diet feeding mice.. Neuropeptides (, 2021. [DOI | PubMed]
52. H Su, G Wang, L Wu, X Ma, K Ying, R Zhang. Transcriptome-wide map of m(6)A circRNAs identified in a rat model of hypoxia mediated pulmonary hypertension.. BMC Genomics (, 2020. [DOI | PubMed]
53. MB Uddin, Z Wang, C Yang. Dysregulations of functional RNA modifications in cancer, cancer stemness and cancer therapeutics.. Theranostics (, 2020. [DOI]
54. JY Youn, WH Dunham, SJ Hong, JDR Knight, M Bashkurov, GI Chen. High-density proximity mapping reveals the subcellular organization of mRNA-associated granules and bodies.. Mol Cell (, 2018. [DOI | PubMed]
55. J Zhu, M Wang, D Hu. Deciphering N(6)-Methyladenosine-Related genes signature to predict survival in lung adenocarcinoma.. BioMed Res Int (, 2020. [DOI | PubMed]
56. JY Qian, J Gao, X Sun, MD Cao, L Shi, TS Xia. KIAA1429 acts as an oncogenic factor in breast cancer by regulating CDK1 in an N6-methyladenosine-independent manner.. Oncogene (, 2019. [DOI]
57. J Chen, K Yu, G Zhong, W Shen. Identification of a m(6)A RNA methylation regulators-based signature for predicting the prognosis of clear cell renal carcinoma.. Cancer Cell Int (, 2020. [DOI | PubMed]
58. F Lanigan, G Gremel, R Hughes, DJ Brennan, F Martin, K Jirstrom. Homeobox transcription factor muscle segment homeobox 2 (Msx2) correlates with good prognosis in breast cancer patients and induces apoptosis in vitro .. Breast Cancer Res (, 2010. [DOI | PubMed]
59. TA Stevens, R Meech. BARX2 and estrogen receptor-alpha (ESR1) coordinately regulate the production of alternatively spliced ESR1 isoforms and control breast cancer cell growth and invasion.. Oncogene (, 2006. [DOI]
60. Y Niu, Z Lin, A Wan, H Chen, H Liang, L Sun. RNA N6-methyladenosine demethylase FTO promotes breast tumor progression through inhibiting BNIP3.. Mol Cancer (, 2019. [DOI | PubMed]
61. D Dixit, Q Xie, JN Rich, JC Zhao. Messenger RNA methylation regulates glioblastoma tumorigenesis.. Cancer Cell (, 2017. [DOI]
62. T Lan, H Li, D Zhang, L Xu, H Liu, X Hao. KIAA1429 contributes to liver cancer progression through N6-methyladenosine-dependent post-transcriptional modification of GATA3.. Mol Cancer (, 2019. [DOI | PubMed]
63. R Miao, CC Dai, L Mei, J Xu, SW Sun, YL Xing. KIAA1429 regulates cell proliferation by targeting c-jun messenger RNA directly in gastric cancer.. J Cell Physiol (, 2020. [DOI]
64. S Fazeli, E Paal, JH Maxwell, KD Burman, ES Nylen, SG Khosla. Salutary response to targeted therapy in anaplastic thyroid cancer.. J Invest Med High impact Case Rep (, 2019. [DOI]
65. Z Zhao, Q Ju, J Ji, Y Li, Y Zhao. N6-methyladenosine methylation regulator RBM15 is a potential prognostic biomarker and promotes cell proliferation in pancreatic adenocarcinoma.. Front Mol Biosci (, 2022. [DOI | PubMed]
66. W Qiu, Y Zhou, H Wu, X Lv, L Yang, Z Ren. RNA Demethylase FTO mediated RNA m(6)A modification is involved in maintaining maternal-fetal interface in spontaneous abortion.. Front Cell Dev Biol (, 2021. [DOI | PubMed]
67. T Stepien, M Brozyna, K Kuzdak, E Motylewska, J Komorowski, H Stepien. Elevated concentrations of SERPINE2/Protease nexin-1 and secretory leukocyte protease inhibitor in the serum of patients with papillary thyroid cancer.. Dis Markers (, 2017. [DOI | PubMed]
68. H Guan, W Liang, J Liu, G Wei, H Li, L Xiu. Transmembrane protease serine 4 promotes thyroid cancer proliferation via CREB phosphorylation.. Thyroid (, 2015. [DOI | PubMed]
69. Y Kushchayeva, K Jensen, A Recupero, J Costello, A Patel, J Klubo-Gwiezdzinska. The HIV protease inhibitor nelfinavir down-regulates RET signaling and induces apoptosis in medullary thyroid cancer cells.. J Clin Endocrinol Metab (, 2014. [DOI]
70. W Chen, J Li, S Peng, S Hong, H Xu, B Lin. Association of total thyroidectomy or thyroid lobectomy with the quality of life in patients with differentiated thyroid cancer with low to intermediate risk of recurrence.. JAMA Surg (, 2022. [DOI]
71. RQ Yang, KL Lou, PY Wang, YY Gao, YQ Zhang, M Chen. Surgical navigation for malignancies guided by near-Infrared-II fluorescence imaging.. Small Methods (, 2021. [DOI | PubMed]
72. M Abhishek, A Renuka, A Ujjwal, C Amit, P Vijay, N Vanita. Atypical posterior reversible encephalopathy syndrome associated with lenvatinib therapy in a patient with metastatic thyroid cancer-a case report.. Cancer Rep (Hoboken) (, 2022. [DOI | PubMed]
73. L Agate, L Puleo, C Giani, L Valerio, E Molinaro, R Elisei. Re: "Symptomatic biliary disorders during lenvatinib treatment for thyroid cancer: An underestimated problem".. Thyroid (, 2021. [DOI]
74. F Ahmaddy, C Burgard, L Beyer, VF Koehler, P Bartenstein, MP Fabritius. (18)F-FDG-PET/CT in patients with advanced, radioiodine refractory thyroid cancer treated with lenvatinib.. Cancers (Basel) (, 2021. [DOI | PubMed]
75. K Alshehri, Y Alqurashi, M Merdad, S Samargandy, R Daghistani, H Marzouki. Neoadjuvant lenvatinib for inoperable thyroid cancer: A case report and literature review.. Cancer Rep (Hoboken) (, 2022. [DOI | PubMed]
76. M Aoyama, H Takizawa, T Otani, S Inoue, N Kawakita, M Tsuboi. Noninvasive monitoring of paclitaxel and lenvatinib efficacy against anaplastic thyroid cancer in orthotopic SCID mouse models using smallanimal FDGPET/CT.. Oncol Rep (, 2020. [DOI]
77. N Arai, H Fujiwara, H Sasaki. Letter to the Editor regarding "Unusual magnetic resonance imaging findings of a glioblastoma arising during treatment with lenvatinib for thyroid cancer".. World Neurosurg (, 2018. [DOI]
78. G Multhoff, K Pfister, C Botzler, A Jordan, R Scholz, H Schmetzer. Adoptive transfer of human natural killer cells in mice with severe combined immunodeficiency inhibits growth of Hsp70-expressing tumors.. Int J Cancer (, 2000. [DOI]
79. JP Veluchamy, N Kok, HJ van der Vliet, HMW Verheul, TD de Gruijl, J Spanholtz. The rise of allogeneic natural killer cells as a platform for cancer immunotherapy: Recent innovations and future developments.. Front Immunol (, 2017. [DOI | PubMed]
80. A Collignon, F Silvy, S Robert, M Trad, S Germain, J Nigri. Dendritic cell-based vaccination: powerful resources of immature dendritic cells against pancreatic adenocarcinoma.. Oncoimmunology (, 2018. [DOI | PubMed]
81. BA Hanks, A Holtzhausen, KS Evans, R Jamieson, P Gimpel, OM Campbell. Type III TGF-beta receptor downregulation generates an immunotolerant tumor microenvironment.. J Clin Invest (, 2013. [DOI]
82. UK Scarlett, MR Rutkowski, AM Rauwerdink, J Fields, X Escovar-Fadul, J Baird. Ovarian cancer progression is controlled by phenotypic changes in dendritic cells.. J Exp Med (, 2012. [DOI | PubMed]
83. D Mittal, MM Gubin, RD Schreiber, MJ Smyth. New insights into cancer immunoediting and its three component phases–elimination, equilibrium and escape.. Curr Opin Immunol (, 2014. [DOI | PubMed]
84. LC Platanias. Mechanisms of type-i- and type-II-interferon-mediated signalling.. Nat Rev Immunol (, 2005. [DOI]
85. JL Benci, B Xu, Y Qiu, TJ Wu, H Dada, C Twyman-Saint Victor. Tumor interferon signaling regulates a multigenic resistance program to immune checkpoint blockade.. Cell (, 2016. [DOI | PubMed]
86. V Shankaran, H Ikeda, AT Bruce, JM White, PE Swanson, LJ Old. IFNgamma and lymphocytes prevent primary tumour development and shape tumour immunogenicity.. Nature (, 2001. [DOI]
87. S Zhao, C Chen, K Chang, A Karnad, J Jagirdar, AP Kumar. CD44 expression level and isoform contributes to pancreatic cancer cell plasticity, invasiveness, and response to therapy.. Clin Cancer Res (, 2016. [DOI]
88. F Gansauge, S Gansauge, A Zobywalski, C Scharnweber, KH Link, AK Nussler. Differential expression of CD44 splice variants in human pancreatic adenocarcinoma and in normal pancreas.. Cancer Res (, 1995
89. M Kaufmann, KH Heider, HP Sinn, G von Minckwitz, H Ponta, P Herrlich. CD44 variant exon epitopes in primary breast cancer and length of survival.. Lancet (, 1995. [DOI]
90. H Kuniyasu, Y Chihara, T Kubozoe, T Takahashi. Co-Expression of CD44v3 and heparanase is correlated with metastasis of human colon cancer.. Int J Mol Med (, 2002. [DOI]
91. Z Li, K Chen, P Jiang, X Zhang, X Li, Z Li. CD44v/CD44s expression patterns are associated with the survival of pancreatic carcinoma patients.. Diagn Pathol (, 2014. [DOI | PubMed]
92. M Muramaki, H Miyake, S Kamidono, I Hara. Over expression of CD44V8-10 in human bladder cancer cells decreases their interaction with hyaluronic acid and potentiates their malignant progression.. J Urol (, 2004. [DOI]
93. M Todaro, M Gaggianesi, V Catalano, A Benfante, F Iovino, M Biffoni. CD44v6 is a marker of constitutive and reprogrammed cancer stem cells driving colon cancer metastasis.. Cell Stem Cell (, 2014. [DOI]
94. SJ Wang, VB Wreesmann, LY Bourguignon. Association of CD44 V3-containing isoforms with tumor cell growth, migration, matrix metalloproteinase expression, and lymph node metastasis in head and neck cancer.. Head Neck (, 2007. [DOI]
95. M Zoller. CD44: can a cancer-initiating cell profit from an abundantly expressed molecule?. Nat Rev Cancer (, 2011. [DOI]
96. V Subramaniam, IR Vincent, M Gilakjan, S Jothy. Suppression of human colon cancer tumors in nude mice by siRNA CD44 gene therapy.. Exp Mol Pathol (, 2007. [DOI]
97. X Huang, Y Sheng, M Guan. Co-Expression of stem cell genes CD133 and CD44 in colorectal cancers with early liver metastasis.. Surg Oncol (, 2012. [DOI]
98. ME Ritchie, B Phipson, D Wu, Y Hu, CW Law, W Shi. Limma powers differential expression analyses for RNA-sequencing and microarray studies.. Nucleic Acids Res (, 2015. [DOI | PubMed]
99. GJ Prud’homme, Y Glinka. Neuropilins are multifunctional coreceptors involved in tumor initiation, growth, metastasis and immunity.. Oncotarget (, 2012. [DOI]
100. X Liu, X Meng, X Peng, Q Yao, F Zhu, Z Ding. Impaired AGO2/miR-185-3p/NRP1 axis promotes colorectal cancer metastasis.. Cell Death Dis (, 2021. [DOI | PubMed]
101. Y Zhou, Y Chen, Y Tan, R Hu, MM Niu. An NRP1/MDM2-targeted d-peptide supramolecular nanomedicine for high-efficacy and low-toxic liver cancer therapy.. Adv Healthc Mater (, 2021. [DOI | PubMed]