Investigation of risk signatures associated with anoikis in thyroid cancer through integrated transcriptome and Mendelian randomization analysis

Article type: Research Article

Keywords: thyroid cancer, anoikis, risk signature, prognosis, Mendelian randomization analysis

Authors: Xiang-Yi Chen, Jia-Ying Lai, Wen-Jun Shen, Dawei Wang, Zhi-Xiao Wei

Affiliations: Department of Nuclear Medicine, The First Affiliated Hospital of Guangxi Medical University, Nanning, China; Department of Medical Engineering, Medical Supplies Center of PLA General Hospital, Beijing, China; Department of Nuclear Medicine, The Sixth Medical Center of PLA General Hospital, Beijing, China

License: Copyright © 2024 Chen, Lai, Shen, Wang and Wei 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/fendo.2024.1458956  |  PubMed: 39568815  |  PMC: PMC11576184

Relevance: Relevant: mentioned in keywords or abstract

Full text: PDF (725 KB)

Introduction

Thyroid cancer (THCA) is a prevalent malignancy of the endocrine system and ranks as the ninth most common cancer worldwide (ref. 1). THCA, primarily originating from follicular cells, constitutes approximately 95% of THCA (ref. 2). The prognosis of THCA patients depends on factors such as tumor size, age, degree of metastasis, among which tumor size closely correlates with poor prognosis (ref. 3). In the United States alone, around 44,000 new cases of THCA are diagnosed annually. Among patients with highly differentiated THCA, about 54% have a low risk of recurrence and can be effectively treated through surgery (ref. 1). Despite having a favorable prognosis overall, once it progresses to locally advanced stages or metastatic radioiodine refractory (RAIR) which exhibits relatively poorer outcomes that significantly impact patient survival rates. Therefore, we believe that identifying novel genes associated with THCA not only plays a crucial role in pathological diagnosis but also provides valuable insights for treatment selection.

Anoikis is a form of cellular apoptosis triggered by the loss or improper attachment (ref. 4). Normal epithelial cells rely on extracellular matrix (ECM) contact for survival and proliferation, whereas cancer cells can bypass this requirement (ref. 4). Consequently, anti-anoikis mechanisms are essential for invasive metastasis and dissemination of cancer cells (ref. 4). Simultaneously, numerous studies have unveiled the association between anoikis and THCA. For instance, Lee JJ et al. discovered that galactin-3 inhibitors suppressed resistance to anoikis and invasion in THCA cells (ref. 5). Jensen K et al. demonstrated that inhibiting the gap junction pathway could enhance sensitivity to anoikis stimuli in THCA cells (ref. 6). Tang M et al. also identified anoikis-related gene CDKN2A as a predictor of prognosis and immune response while mediating proliferation and migration in THCA (ref. 7). Other investigations have revealed that MicroRNA-363-3p inhibits anoikis in human THCA by targeting integrin α6 (ref. 8), suggesting a close relationship between anoikis and malignant progression of THCA. Therefore, it is crucial to authenticate predictive genes for effective therapy against THCA.

In order to further investigate the role of anoikis-related genes (ARGs) in the occurrence and development of THCA, we utilized sample information from the TCGA dataset to conduct a comprehensive analysis on survival-associated genes. Through bioinformatics analysis, we successfully identified six prognostic genes for THCA, namely LRRC75A, METTL7B, ADRA1B, TPD52L1, TNFRSF10C, and CXCL8. Additionally, the risk signature was constructed to accurately predict patients’ survival time. To elucidate the causal relationship between these prognostic genes and THCA risk, a two-sample Mendelian randomization (MR) analysis was performed using genome-wide association studies (GWAS) data from the OpenGWAS database. The results revealed a significant causal association between CXCL8 and THCA. Collectively, our findings provide a solid foundation for future research on THCA and offer new insights into personalized treatment strategies.

Materials and methods

Data acquisition

A total of 497 cases of THCA and 56 control cases were gathered from the UCSC Xena database (https://xenabrowser.net/datapages/). Among these, 496 THCA samples were accompanied by survival information and comprehensive clinical characteristics, encompassing age, M/N/T stage, and stage, etc. A set of 338 ARGs were acquired from the published literature (ref. 9). Additionally, the GWAS data for THCA and eQTL GWAS data for ARGs in MR analysis were all collected from the OpenGWAS database (https://gwas.mrcieu.ac.uk/). The THCA’s GWAS dataset, namely ebi-a-GCST90018929, consisted of 1,054 cases with THCA and 490,920 control cases (European), encompassing a total of 24,198,226 single-nucleotide polymorphisms (SNPs).

Differential expression analysis and construction of the WGCNA network

“DESeq2” package (v1.36.0) (ref. 10) was utilized to mine the differentially expressed genes (DEGs) between THCA and control cases in the TCGA-THCA dataset with log₂fold change (FC) > 1 and adj.p< 0.05. To identify key module genes associated with ARGs in THCA, the ssGSEA algorithm in the “GSVA” package (v1.44.5) (ref. 11) was utilized to calculate ARG scores for THCA and control specimens in the TCGA-THCA dataset. Subsequently, depending on the gene expression matrix of all samples within the TCGA-THCA dataset and ARG score (clinical trait), a co-expression network was established using “WGCNA” package (v1.70-3) (ref. 12). The samples in the TCGA-THCA dataset were subjected to clustering analysis in order to identify and exclude any outlier samples, thereby ensuring the accuracy of subsequent analyses. Afterthat, in order to maximize the adherence of gene interactions to a scale-free distribution, we employed the expression matrix encompassing all genes as input data and conducted computations to determine the optimal soft threshold. By calculating gene adjacency and deriving a coefficient of dissimilarity based on gene similarity, a systematic clustering tree of genes was acquired. Finally, all samples were treated as individual traits and subjected to Pearson correlation analysis with modules obtained from WGCNA analysis in order to identify the module gene exhibiting the highest correlation. The intersection of DEGs and key module genes was performed using the “ggvenn” package (v0.1.9) for visualization, resulting in a set of differentially expressed ARGs (DE-ARGs).

Functional enrichment analysis for DE-ARGs

The protein-protein interaction (PPI) network of DE-ARGs was established through the STRING database (https://www.string-db.org/). Subsequently, “clusterProfiler” package (v4.7.1.001) (ref. 13) was employed to perform GO and KEGG enrichment analyses on the DE-ARGs, aiming to mine shared functions and associated pathways among genes.

Screening of prognostic genes and construction of the risk signature

Subsequent analysis involved performing univariate Cox regression using “survival” package (v3.4-0) (ref. 14) on DE-ARGs within THCA samples from TCGA-THCA dataset. Genes that met the criteria of hazard ratio (HR) ≠ 1 and p< 0.05 were categorized as survival-associated genes. These genes then underwent proportional hazards (PH) assumption testing. The genes that met the PH assumption (p > 0.05) were further integrated into a LASSO regression analysis using the “glmnet” package (v4.1-6) (ref. 15). Optimal model fit was achieved when the lambda value was minimized, identifying the prognostic genes for this study. The THCA samples in TCGA-THCA dataset were randomly allocated to the training set and the internal validation set at a ratio of 7:3 for constructing the risk signature. Based on the prognostic genes and coefficients obtained from the LASSO analysis, the risk score for each THCA sample was calculated according to the formula below:

Following the calculation of risk scores for each THCA patient, they were categorized into high-risk and low-risk groups via the median risk score. Subsequently, the Kaplan-Meier (K-M) curves for survival analysis between these two risk groups were plotted using “survminer” (v0.4.9) (ref. 14) package, followed by a comparison of survival differences through log-rank test. Besides, the receiver operating characteristic (ROC) curves for the risk signature were drawn by “survivalROC” package (v1.0.3) (ref. 16), with survival times of 1, 3, and 5 years as the time points. To further assess the generality of the risk signature, we validated it using the same method as described above in the internal validation set.

Independent prognostic analysis

The nomogram was generated using the “rms” package (v6.3-0) (ref. 17), and then calibration curves were built to estimate the accuracy of nomograms in predicting 1/3/5-year survival in THCA patients. In order to assess the relevant between prognostic genes and clinical characteristics, we integrated the clinical information (including age, T/M/N stage, and disease stage) of THCA patients obtained from the TCGA database to examine the differences in prognostic gene expression among different clinical features.

Immune infiltration analysis

The infiltration of immune cell types in patients with THCA was estimated with CIBERSORT (https://cibersort.stanford.edu/), and the Wilcoxon test was employed to make comparisons between the high- and low-risk groups. To investigate the association between tumor immune evasion and risk score, we calculated and compared the TIDE scores of THCA samples through TIDE’s official website (http://tide.dfci.harvard.edu/) and Wilcoxon test, respectively. In addition, Spearman correlation analysis was employed to examine the relationship between TIDE and risk scores.

Functional enrichment analysis for prognostic genes

Using “h.all.v7.4.symbols.gmt” obtained from the MSigDB database (http://www.broadinstitute.org/msigdb) as the reference, we employed “GSVA” package (v1.44.5) (ref. 11) to calculate scores for various pathways enriched in high- and low-risk groups (p< 0.05). Subsequently, we utilized the “limma” package to identify significantly different pathways between THCA samples in two risk groups. Moreover, the “psych” package (v2.2.9) (ref. 18) was utilized for the analysis of Spearman correlation between the risk score and all genes. Depending on risk score and coefficient serving as ranking criteria, gene set enrichment analysis (GSEA) was conducted by the “clusterProfiler” package (|NES| > 1 & adj.p< 0.05). The “c2.cp.kegg.v2023.1.Hs.symbols” gene set obtained from the GSEA website (http://www.gsea-MSigdb.org/gsea/MSigdb) was served as the background gene set.

Drug sensitivity analysis

The oncoPredict package (v0.2) (ref. 19) was utilized to calculate the IC50 values of 198 chemotherapeutic/targeted therapy drugs provided by GDSC database (https://www.cancerrxgene.org/) on THCA samples. Subsequently, we conducted a comparative analysis of the IC50 values for 198 drugs between two distinct groups.

Study design of MR analysis

Using a two-sample MR analysis, we examined the causal relationship between prognostic genes and THCA. Therefore, the selection of effective instrumental variables (IVs) should adhere to three crucial assumptions: (1) IVs must exhibit strong associations with exposure; (2) IVs must be independent of confounding factors; and (3) The association between IVs and outcome should be solely mediated by exposure factors.

MR analysis of THCA and prognostic genes

The extract_instruments function in the TwoSampleMR package (v0.5.7) (ref. 20) was utilized to identify IVs (SNPs) that exhibited significant correlation with prognostic genes (exposure factors), employing a screening criterion of p< 5*10^-8. Subsequently, the clump parameter was set to TRUE, and IVs in linkage disequilibrium (LD) were removed by adjusting the parameters as follows: r² = 0.001; kb = 100. To address potential issues arising from IVs exhibiting significant correlation with THCA (outcome), we employed the harmonise_data function in the TwoSampleMR package to standardize effect equiele and effect size. MR analysis was performed using 5 algorithms (MR egger, weighted median, inverse variance weighted [IVW], simple mode, and weighted mode). Notably, the primary outcome considered was based on IVW results, where a P< 0.05 indicated a significant causal relationship. Finally, ensuring accuracy in directionality was crucial and achieved through utilizing the TwoSampleMR package to conduct Steiger directivity analysis.

Sensitivity analysis

A sensitivity analysis was undertaken to assure the robustness of MR analysis. Initially, heterogeneity was evaluated using Cochran’s Q statistic and its corresponding p value. Then, the MR egger intercept was examined to investigate potential horizontal pleiotropy. It is worth noting that if p< 0.05 in egger intercept, it necessitates reconsideration of the study design. Furthermore, to demonstrate reliability in our analysis results, a stepwise elimination process was employed where each SNP was gradually removed and its impact on outcome variables assessed. This leave-one-out sensitivity analysis approach served as an additional validation for our findings.

THCA patients and tissue specimens

A total of 15 pairs of THCA tumor and corresponding adjacent normal tissues were collected from First Affiliated Hospital of Guangxi Medical University between August 2024 and September 2024. Samples were snap-frozen and stored at −80°C until used in real-time qPCR (RT-qPCR) experiments. The Research Ethics Committee of The First Affiliated Hospital of Guangxi Medical University approved this study, which was consistent with the Declaration of Helsinki.

RNA Isolation, cDNA Synthesis, and RT-qPCR

Total RNAs were isolated by Triquick Reagent (Trizol Substitute, Solarbio, China). RNA (1 μg), quantified by NanoDrop2000 (Thermo Fisher Scientific), was reversely transcribed to cDNA using the first-strand cDNA synthesis kit (Monad, China). Quantitative PCR was applied using the SYBR Green dye (Yeasen, China) on Real-Time PCR System (Agilen, United States). All primers were synthesized by Sangon Biotech and their sequences were listed in Supplementary Table S1. ACTB served as an internal control. The 2^-ΔΔCt method was used to determine the relative expression between cancer and normal tissue for each selected prognostic gene.

Statistical analysis

This study was statistically analyzed using the R software package (v4.2.1). Various variables, including expression quantity and infiltration ratio, were evaluated for differences between groups using either the Wilcoxon test or t-test. The significance of discrepancies between groups was typically defined by a p value or adj.p value below 0.05, unless otherwise specified.

Results

Identification of DE-ARGs in the TCGA-THCA dataset

Through differential expression analysis, a set of 3,772 DEGs were detected in the TCGA-THCA cohort, comprising 1,611 down-regulated and 2,161 up-regulated genes. The volcano plot displayed all DEGs, while the heat map depicted the top 80 genes with the most significant up-regulation and down-regulation between control and THCA samples (f1, Supplementary Table S2). The ARG scores in THCA samples and healthy control samples were compared using the Wilcoxon test, revealing statistically significant differences (p< 0.05, f2). To screen out genes with high ARG correlation, we constructed a WGCNA network based on ARG scores. Depending on the results of the cluster analysis, no outlier samples were identified in the TCGA-THCA dataset, indicating its suitability for conducting WGCNA analysis (f2). When the R2 value approached the threshold of 0.85 (red line) and the average connectivity in the graph on the right tends towards 0, it indicated that the network was approaching a scale-free distribution (f2). Consequently, we determined the optimal soft threshold (β) to be 8. Furthermore, the minimum number of genes was set to 100 according to the standard of the hybrid dynamic tree cutting algorithm, resulting in a total of 14 modules obtained (f2). Among the 14 modules, MEbrown exhibited the strongest association with ARGs score (cor = 0.89 & p< 0.001) (f2). Consequently, MEbrown module was identified as the pivotal module in this study, encompassing a total of 1,443 genes and designated as the key module gene. Finally, 570 DE-ARGs were acquired by overlapping 3,772 DEGs and 1,443 ARGs-related genes (f2).

Investigation of potential action mechanism for DE-ARGs

Based on the functional enrichment analysis, wound healing, developmental growth involved in morphogenesis and pathways associated with nervous system function, such as axon extension, synapse organization, axonogenesis, were enriched by DE-ARGs (f3). Moreover, they were also associated with ECM-receptor interaction, transcriptional misregulation in cancer, cell adhesion molecules, p53 signaling pathway according to KEGG pathways (f3). In order to investigate potential interactions among DE-ARGs, we constructed a PPI network. The network comprised a total of 405 nodes and 1,371 interaction pairs, revealing a highly intricate interplay between 570 ARGs (f3).

Risk signature based on six prognostic genes for THCA

After performing univariate Cox regression analysis, a total of 40 survival-associated genes were identified (HR ≠ 1 and p< 0.05) (Supplementary Table S3). These genes were subsequently subjected to tests for the PH assumption, which they all satisfied (p > 0.05) (Supplementary Table S4). According to LASSO regression analysis, six prognostic genes (LRRC75A, METTL7B, ADRA1B, TPD52L1, TNFRSF10C, and CXCL8) were used to construct the risk signature with the best performance (f4). The THCA samples in the training set were ranked based on their risk scores, and a risk curve was plotted (f4). This curve demonstrated a higher number of deceased samples in the high-risk group compared to the low-risk group, indicating that the mortality rate was lower among low-risk THCA patients. Additionally, f4 displayed expression levels of prognostic genes in both high and low-risk groups, in which LRRC75A, METTL7B, ADRA1B, TPD52L1, and TNFRSF10C were low expressed in high-risk group, and CXCL8 showed an opposite trend. The K-M curves revealed a statistically significant disparity in the survival rates of patients with THCA between the high- and low-risk groups (p< 0.05), with the high-risk group exhibiting a higher mortality rate (f4). The area under the curve (AUC) values of the ROC curve were 0.938, 0.890, and 0.863 at 1, 3, and 5 years, respectively, indicating a strong ability for predicting survival time (f4). Furthermore, the risk signature constructed by the prognostic gene in the internal validation set exhibited consistent findings with those observed in the training set, thereby providing further evidence of its validity and universality (Supplementary Figure S1A–D).

Risk score was an independent predictor of THCA

Based on the clinical correlation analysis, the expression levels of CXCL8 exhibited significant variations among different subtypes of N stage (0/1) and T stage (1/2/3/4), while METTL7B demonstrated differential expression between the TNM stages (I/II/III/IV) and T stage (1/2/3/4), and TPD52L1 displayed alterations at M stage (0/1) (Supplementary Figure S2A–E, Supplementary Table S5). To demonstrate the clinical value of the risk score more directly, we constructed a nomogram model that mapped the risk score onto the total number axis to determine patients’ overall scores and subsequently derive their corresponding predicted survival probabilities (f5). Meanwhile, the calibration curve of the nomogram demonstrated its excellent predictive capability (Supplementary Figure S2F).

Relevance of prognostic genes and immune environment

To investigate the correlation between risk signature and the immune microenvironment in THCA patients, an analysis of immunoinfiltration was conducted. f6 exhibited the proportion of immune infiltration cells in the high- and low-risk groups. The high-risk group exhibited significantly elevated expression of activated dendritic cells, eosinophils, monocytes, neutrophils, and activated memory CD4 T cells among the total of 9 kind of differential immune cells. Conversely, the low-risk group demonstrated highly expressed M2 macrophages, resting mast cells, resting NK cells, and gamma delta T cells (f6). Correlation analysis revealed that CXCL8 displayed the strongest positive correlation with resting dendritic cells (cor = 0.5) (f6). What’s more, the TIDE scores presented significant disparities between the low- and high-risk groups, with the latter demonstrating a heightened likelihood of immune evasion (f6). Additionally, there was a positive correlation observed between the risk score and TIDE score (f6). Therefore, the risk signatures related to anoikis can be utilized for evaluating immune infiltration and immune evasion in patients with THCA.

Exploration of potential functions of prognostic genes

To determine the pathways for which prognostic gene enrichment differed between high- and low-risk groups, we conducted GSEA and GSVA enrichment analyses. f7 displayed the pathways exhibiting significant disparities between two risk groups, highlighting the most prominent distinction where TNFA, IL2 STAT5, KRAS, P53, IL6 JAK STAT3 signaling pathways, inflammatory, and apoptosis were significantly inhibited while oxidative phosphorylation, fatty acid metabolism, and adipogenesis is significantly activated (f7). According to the results of GSEA, risk score was associated with chemokine, T cell receptor, and nod like receptor signaling pathways, cytokine receptor interaction, and hematopoietic cell lineage (f7).

Comparison of drug sensitivity between groups

Based on the data from GDSC database, the IC50 values of 107 drugs exhibited significant variations between two risk groups (high and low), implying potential divergent therapeutic effects in patients at different stages of tumor progression. For the purpose of demonstration, we presented the top10 drugs exhibiting significant differences in f8, including ABT737_1910, AZD3759_1915, Erlotinib_1168, etc. In conclusion, our findings suggested that the risk score played a beneficial role in predicting drug effectiveness.

Causal effect of prognostic genes on THCA

To investigate the causal relationship between prognostic genes and THCA, we conducted a two-sample MR analysis, utilizing prognostic genes as exposure factors and THCA as the outcome. Unfortunately, no GWAS IDs were found for ADRA1B, TPD52L1, and TNFRSF10C, hence they were excluded from MR analysis. According to the results from IVW method, only CXCL8 demonstrated a causal relationship with THCA and was identified as a risk factor for THCA [odds ratio (OR) = 1.1520, 95% confidence interval (CI): 1.0333-1.2843, p = 0.0108]. There was no significant causal association between LRRC75A or METTL7B and THCA (p > 0.05) (T1, f9). The Steiger direction check was confirmed as TRUE, indicating the presence of correct causation and the absence of reverse causation (Supplementary Table S6). The MR analysis included a total of 27 SNPs, all of which exhibited F value exceeding 10, indicating their suitability as strong IVs (T2). Meanwhile, the risk effects of the estimated exposure factors at each SNP site were combined to predict the diagnostic efficacy on THCA within the forest plot. As depicted in f9, IVW analysis revealed that the impact values of CXCL8 on THCA were all greater than 0. In order to assess the randomness of the analysis, a funnel plot was constructed by combining the β coefficient and standard error (SE) for each IVs (f9). The distribution of IVs in the plot appeared slightly uneven, which could be attributed to heterogeneity and potential bias resulting from a limited number of IVs. Consequently, a subsequent test for heterogeneity was conducted, yielding the p value greater than 0.05, indicating it had no impact on the overall findings (Supplementary Table S7). Moreover, the horizontal pleiotropy test conducted on CXCL8 indicated a p value exceeding 0.05, suggesting the absence of any significant horizontal pleiotropic effects (Supplementary Table S8). The absence of any confounding factors was found in this study, thereby further enhancing the robustness of our analysis. Additionally, there were no abnormal SNP results observed in relation to the leave-one-out analysis (f9). In summary, these findings suggested a significant causal relationship between CXCL8 and THCA, indicating that upregulation of CXCL8 expression was associated with an increased risk of THCA.

Verification of prognostic genes expression by qPCR

The expression of the six prognostic genes was then verified in tumor and normal samples by qPCR. Consistent with the prediction, the results showed that the expression levels of the selected genes were remarkable elevated in THCA patients (p< 0.01) (f10).

Discussion

The incidence of THCA is steadily increasing recent decades, making it the most prevalent endocrine cancer with relatively low mortality rates, partly due to overdiagnosis over world (ref. 21, ref. 22). However, almost 10% of THCA exhibit recurrent or metastatic disease, the early appropriate treatment and accurate diagnosis of metastatic TCHA is still challenging (ref. 23, ref. 24). When cells fail to attach to the ECM or adhere in inappropriate locations, they undergo a specific form of apoptosis known as anoikis. If the anoikis process is not properly executed, it may result in the rapid proliferation of cells in ectopic sites. This dysregulation of anoikis has become a hallmark of cancer cells, contributing to their metastasis to distant organs (ref. 25, ref. 26). Metastasis, particularly to the lungs and bones, is the primary cause for worsen survival (ref. 27). Therefore, constructing risk signature associated with metastasis can offer valuable insights for disease management and early intervention for THCA.

In this study, we identified a total of 3,772 DEGs between normal and tumor samples from the TCGA-THCA dataset. After intersecting with ARGs, we obtained 570 DE-ARGs. Functional analysis revealed that these DE-ARGs primarily enriched in P53 signaling pathway, focal adhesion, transcriptional dysregulation, and proteoglycan pathways in cancer, highlighting the pivotal role of ARGs in tumorigenesis. Subsequently, utilizing univariate Cox regression, PH assumption test and LASSO analysis techniques, we further identified six prognostic genes: LRRC75A, METTL7B, ADRA1B, TPD52L1, TNFRSF10C, and CXCL8. The PCR results showed that all six prognostic genes were significantly overexpressed in THCA tissues. Notably, these findings were consistent with the differential expression analysis, with each gene exhibiting a log₂FC greater than 1. These findings confirm our initial conclusions and underscore the importance of these genes in the context of THCA.

Leucine-rich repeat-containing protein 75A (LRRC75A) is a LRRC-rich protein belonging to the LRRC superfamily, yet its functionality remains largely unknown (ref. 28). However, based on scRNA-seq profiling, Miura T et al. have suggested that LRRC75A may positively regulate VEGF expression and serve as a potential biomarker for predicting BM-MSC angiogenesis under ischemic conditions (ref. 28). Notably, some effective therapeutic drugs for THCA act through the angiogenesis pathway involving VEGFR2, the receptor of VEGF (ref. 29). Several studies have demonstrated the high expression of Methyltransferase-like 7B (METTL7B) in THCA, which is consistent with our findings. It is suggested that METTL7B may promote TGF-β1-induced epithelial-mesenchymal transition (EMT), thereby enhancing the migration and invasion of THCA cells (ref. 30, ref. 31). Adrenoceptor alpha 1B (ADRA1B) belongs to the adrenergic receptor alpha1 (ADRA1) subfamily, and adrenergic receptor antagonists have been reported to play a role in treating various cancers such as prostate and breast cancer, besides its significant upregulation in THCA samples (ref. 32–ref. 34). Therefore, we conclude that ADRA1B plays a crucial role in the occurrence and progression of THCA; however, further investigation is required to elucidate its underlying mechanism. Tumor protein D52-like 1 (TPD52L1) belongs to the TPD52 family; previous research has shown that TPD52 can act as a novel regulator of LKB1-AMPK pathway by negatively regulating AMPK activity and thus affecting cancer cell metabolism (ref. 35). As early as 2014, Andrade BM’s team proposed that AMPK regulated iodide and glucose uptake in normal thyroid cells and was highly activated in THCA (ref. 36, ref. 37). However, immunohistochemical analysis revealed no activation of the AMPK pathway in adenoma tissue compared to non-tumor tissue from the same patient (ref. 36). Nevertheless, we still hypothesize that TPD52L1 may be involved in the development of THCA through its regulation of the AMPK pathway, and further investigation into its mechanism is urgently needed. Additionally, TNF receptor superfamily member 10C (TNFRSF10C) inhibits intracellular signaling pathways for apoptosis, and its high methylation is significantly associated with an increased risk in CRC patients (ref. 38). Although the precise role of TNFRSF10C in THAC remains unclear, further investigation is needed to explore the potential link between the significantly inhibited TNF-α/NF-κB signaling pathway identified in GSVA analysis and TNFRSF10C as one of the decoy receptors for TNF-associated apoptosis-inducing ligand (TRAIL). C-X-C motif chemokine ligand 8 (CXCL8 or IL-8) is a chemokine that exerts autocrine or paracrine regulation on tumor proliferation, invasion, and migration. It has been implicated in various cancers such as breast cancer, prostate cancer, lung cancer, and others (ref. 39). The relationship between CXCL8 and thyroid gland was confirmed as early as 1992. Experimental evidence has demonstrated that the presence of CXCL8 is essential for inducing EMT in THCA cells (ref. 40). Recent studies have further revealed higher levels of IL-6 and CXCL8 in THCA tissues compared to normal tissues, which aligns with our findings (ref. 41). Furthermore, GSEV result spot that tumor promotive signal such KRAS, EMT and IL-6/JAK/STAT3 are restricted under anoikis-related condition. Neoplastic diseases involve constant interactions, adaptations, and evolutions between tumor cells at different stages, involving surrounding normal cells, immune cells, integrity with sequential activation tumorigenesis signa-pathways (ref. 42, ref. 43). During the early stages of tumor detachment from extracellular matrix environment, types are linked to promote progression (growth/infiltration into surrounding tissues) that could be summarized among T stage category. Subsequently, when overwhelming metastasis was aroused, the EMT-related and other pathways is initially activated by detachment from the extracellular matrix into circulation (lymph or blood circulation), triggering activation of anoikis-related signal-pathways. Our findings demonstrate that EMT and anoikis are two sequential and independent procedures promote tumorigenesis. This is also the main reason for further MR analysis after the construction of the prognostic signature contained anoikis-related genes.

The risk signature was constructed based on prognostic genes, and the THCA samples were stratified into high- and low-risk groups. The prognosis of patients in the low-risk group exhibited significantly better outcomes compared to those in the high-risk group. Utilizing the ROC curve analysis, we observed that the risk signature demonstrated a robust predictive ability for 1, 3, and 5-year survival probabilities in THCA patients. Furthermore, we developed a nomogram model incorporating risk scores to further elucidate the clinical significance of these risk signatures. In our clinical correlation analysis, we found that risk scores varied to different extents across distinct clinical characteristics. Consequently, can risk scores be employed to differentiate between subtypes of THCA patients? If feasible, this study would provide valuable insights for staging diagnosis and personalized treatment strategies for THCA patients; thus, warranting further investigation.

It is widely acknowledged that immunotherapy holds great promise as a treatment modality for various tumors. To identify personalized treatments for patients with different risk levels of THCA, we conducted an analysis on the disparities in the immune microenvironment between these two groups. Our immunoinfiltration analysis revealed significant discrepancies in the proportions of infiltration among nine types of immune cells, including M2 macrophages, NK cells, and diverse T cells. Shi ZY et al. (ref. 44), discovered a correlation between CD4 T cell count and clinical treatment efficacy, as well as adaptability and recovery potential in THCA patients; this correlation exhibited a negative trend. In our findings, we observed that the proportion of activated memory CD4 T cells were significantly higher in the high-risk group compared to the low-risk group. On the one hand, M2 macrophages have obvious immunosuppressive function to mediate immune escape and promote the proliferation and metastasis of THCA; on the other hand, it can also secrete VEGF and EGF to promote tumor microangiogenesis (ref. 45). Through the calculation and comparison of TIDE scores, we observed significant disparities between high- and low-risk groups, with samples from the former exhibiting a higher likelihood of immune escape. Moreover, there was a positive correlation between risk score and TIDE score. Consequently, we posit that aberrant infiltration of immune cells may foster the progression of THCA, thereby offering a theoretical foundation for exploring THCA immunotherapy.

However, our research has certain limitations. Firstly, the samples we analyzed were obtained from public databases with a small sample size and limited clinical information, which may restrict the extent to which our risk model can be researched and validated in terms of its clinical value. Secondly, although the expression of these six hub genes had been verified by clinical samples, further studies are needed to confirm the role of these prognostic genes in THCA. Moreover, additional investigation is urgent to uncover mechanisms of anoikis and its association with THCA.

Conclusion

In conclusion, we employed bioinformatics methods to screen six prognostic genes (LRRC75A, METTL7B, ADRA1B, TPD52L1, TNFRSF10C, and CXCL8) associated with anoikis for THCA. Utilizing these prognostic genes as a basis, we developed a risk signature to accurately predict the prognosis of THCA patients and made predictions regarding potential functions and targeted drugs. Significantly advancing the field of research on THCA prognostication is our novel integration of MR with transcriptome analysis in order to identify prognostic genes exhibiting significant causal relationships. Of particular interest is our discovery through MR analysis that there exists a substantial causal relationship between CXCL8 and THCA; indeed, CXCL8 emerges as a major risk factor for this disease. This finding not only provides a theoretical foundation for further exploration into the association between ARGs and THCA but also offers an innovative target for enhancing prognosis and facilitating precise and effective treatment strategies for individuals afflicted by this condition.

References

1. L Boucai, M Zafereo, ME Cabanillas. Thyroid cancer: A review.. Jama. (, 2024. [DOI]
2. F Nabhan, PH Dedhia, MD Ringel. Thyroid cancer, recent advances in diagnosis and therapy.. Int J Cancer. (, 2021. [DOI]
3. J Li, Z Li, P Zhao. Diagnosis and prognosis of thyroid cancer by immune-related genes.. Am J Clin Oncol. (, 2024. [DOI | PubMed]
4. ML Taddei, E Giannoni, T Fiaschi, P Chiarugi. Anoikis: an emerging hallmark in health and diseases.. J Pathol. (, 2012. [DOI]
5. JJ Lee, YC Hsu, YS Li, SP Cheng. Galectin-3 inhibitors suppress anoikis resistance and invasive capacity in thyroid cancer cells.. Int J Endocrinol. (, 2021. [DOI | PubMed]
6. K Jensen, A Patel, J Klubo-Gwiezdzinska, A Bauer, V Vasko. Inhibition of gap junction transfer sensitizes thyroid cancer cells to anoikis.. Endocr Relat Cancer. (, 2011. [DOI]
7. M Tang, W Luo, Y Zhou, Z Zhang, Z Jiang. Anoikis-related gene CDKN2A predicts prognosis and immune response and mediates proliferation and migration in thyroid carcinoma.. Transl Oncol. (, 2024. [DOI | PubMed]
8. Y Pan, X Zhu, K Wang, Y Chen. MicroRNA-363-3p suppresses anoikis resistance in human papillary thyroid carcinoma via targeting integrin alpha 6.. Acta Biochim Biophys Sin (Shanghai). (, 2019. [DOI]
9. X Diao, C Guo, S Li. Identification of a novel anoikis-related gene signature to predict prognosis and tumor microenvironment in lung adenocarcinoma.. Thorac Cancer. (, 2023. [DOI]
10. MI Love, W Huber, S Anders. Moderated estimation of fold change and dispersion for RNA-seq data with DESeq2.. Genome Biol. (, 2014. [DOI | PubMed]
11. S Hänzelmann, R Castelo, J Guinney. GSVA: gene set variation analysis for microarray and RNA-seq data.. BMC Bioinf. (, 2013. [DOI]
12. P Langfelder, S Horvath. WGCNA: an R package for weighted correlation network analysis.. BMC Bioinf. (, 2008. [DOI]
13. G Yu, LG Wang, Y Han, QY He. clusterProfiler: an R package for comparing biological themes among gene clusters.. Omics. (, 2012. [DOI]
14. Y Shi, Y Wang, H Dong, K Niu, W Zhang, K Feng. Crosstalk of ferroptosis regulators and tumor immunity in pancreatic adenocarcinoma: novel perspective to mRNA vaccines and personalized immunotherapy.. Apoptosis. (, 2023. [DOI]
15. S Engebretsen, J Bohlin. Statistical predictions with glmnet.. Clin Epigenetics. (, 2019. [DOI | PubMed]
16. S Zhang, L Sun, D Cai, G Liu, D Jiang, J Yin. Development and validation of PET/CT-based nomogram for preoperative prediction of lymph node status in esophageal squamous cell carcinoma.. Ann Surg Oncol. (, 2023. [DOI]
17. M Liu, Y Gao, Y Yuan, X Liu, Y Wang, L Li. A comprehensive study of clinicopathological and genetic features of neuronal intranuclear inclusion disease.. Neurol Sci. (, 2023. [DOI]
18. S Grieder, MD Steiner. Algorithmic jingle jungle: A comparison of implementations of principal axis factoring and promax rotation in R and SPSS.. Behav Res Methods. (, 2022. [DOI | PubMed]
19. D Maeser, RF Gruener, RS Huang. oncoPredict: an R package for predicting in vivo or cancer patient drug response and biomarkers from cell line screening data.. Brief Bioinform. (, 2021. [DOI | PubMed]
20. Y Huang, W Lin, X Zheng. Causal association between 637 human metabolites and ovarian cancer: a mendelian randomization study.. BMC Genomics. (, 2024. [DOI | PubMed]
21. Y Liu, J Wang, X Hu, Z Pan, T Xu, J Xu. Radioiodine therapy in advanced differentiated thyroid cancer: Resistance and overcoming strategy.. Drug Resist Updat. (, 2023. [DOI | PubMed]
22. M Li, LD Maso, S Vaccarella. Global trends in thyroid cancer incidence and the impact of overdiagnosis.. Lancet Diabetes Endocrinol. (, 2020. [DOI]
23. L Fugazzola, R Elisei, D Fuhrer, B Jarzab, S Leboulleux, K Newbold. 2019 European thyroid association guidelines for the treatment and follow-Up of advanced radioiodine-Refractory thyroid cancer.. Eur Thyroid J. (, 2019. [DOI]
24. L Shobab, C Gomes-Lima, A Zeymo, R Feldman, J Jonklaas, L Wartofsky. Clinical, pathological, and molecular profiling of radioactive iodine refractory differentiated thyroid cancer.. Thyroid. (, 2019. [DOI]
25. J Wang, Z Luo, L Lin, X Sui, L Yu, C Xu. Anoikis-associated lung cancer metastasis: mechanisms and therapies.. Cancers (Basel). (, 2022. [DOI | PubMed]
26. K Zhao, Z Wang, T Hackert, C Pitzer, M Zöller. Tspan8 and Tspan8/CD151 knockout mice unravel the contribution of tumor and host exosomes to tumor progression.. J Exp Clin Cancer Res. (, 2018. [DOI | PubMed]
27. EA Toraih, MH Hussein, M Zerfaoui, AS Attia, A Marzouk Ellythy, A Mostafa. Site-specific metastasis and survival in papillary thyroid cancer: the importance of brain and multi-organ disease.. Cancers (Basel). (, 2021. [DOI | PubMed]
28. T Miura, T Kouno, M Takano, T Kuroda, Y Yamamoto, S Kusakawa. Single-cell RNA-seq reveals LRRC75A-expressing cell population involved in VEGF secretion of multipotent mesenchymal stromal/stem cells under ischemia.. Stem Cells Transl Med. (, 2023. [DOI]
29. T Enokida, M Tahara. Management of VEGFR-targeted TKI for thyroid cancer.. Cancers (Basel). (, 2021. [DOI | PubMed]
30. D Ye, Y Jiang, Y Sun, Y Li, Y Cai, Q Wang. METTL7B promotes migration and invasion in thyroid cancer through epithelial-mesenchymal transition.. J Mol Endocrinol. (, 2019. [DOI | PubMed]
31. J Ali, W Liu, W Duan, C Liu, J Song, S Ali. METTL7B (methyltransferase-like 7B) identification as a novel biomarker for lung adenocarcinoma.. Ann Transl Med. (, 2020. [DOI | PubMed]
32. RS Freudenberger, J Kim, I Tawfik, FA Sonnenberg. Optimal medical therapy is superior to transplantation for the treatment of class I, II, and III heart failure: a decision analytic approach.. Circulation. (, 2006. [DOI]
33. AM Harris, BW Warner, JM Wilson, A Becker, RG Rowland, W Conner. Effect of alpha1-adrenoceptor antagonist exposure on prostate cancer incidence: an observational cohort study.. J Urol. (, 2007. [DOI]
34. LK Zhong, XY Deng, F Shen, WS Cai, JH Feng, XX Gan. Identification of a 3-gene prognostic index for papillary thyroid carcinoma.. Front Mol Biosci. (, 2022. [DOI | PubMed]
35. Y Chen, C Peng, W Tan, J Yu, J Zayas, Y Peng. Tumor protein D52 (TPD52) affects cancer cell metabolism by negatively regulating AMPK.. Cancer Med. (, 2023. [DOI]
36. BM Andrade, DP de Carvalho. Perspectives of the AMP-activated kinase (AMPK) signalling pathway in thyroid cancer.. Biosci Rep. (, 2014. [DOI | PubMed]
37. AP Vidal, BM Andrade, F Vaisman, J Cazarin, LF Pinto, MM Breitenbach. AMP-activated protein kinase signaling is upregulated in papillary thyroid cancer.. Eur J Endocrinol. (, 2013. [DOI]
38. C Zhou, R Pan, H Hu, B Li, J Dai, X Ying. TNFRSF10C methylation is a new epigenetic biomarker for colorectal cancer.. PeerJ. (, 2018. [DOI | PubMed]
39. Q Liu, A Li, Y Tian, JD Wu, Y Liu, T Li. The CXCL8-CXCR1/2 pathways in cancer.. Cytokine Growth Factor Rev. (, 2016. [DOI | PubMed]
40. M Rotondi, F Coperchini, F Latrofa, L Chiovato. Role of chemokines in thyroid cancer microenvironment: is CXCL8 the main player?. Front Endocrinol (Lausanne). (, 2018. [DOI | PubMed]
41. L Zhang, S Xu, X Cheng, J Wu, Y Wang, W Gao. Inflammatory tumor microenvironment of thyroid cancer promotes cellular dedifferentiation and silencing of iodide-handling genes expression.. Pathol Res Pract. (, 2023. [DOI | PubMed]
42. R Tedja, AB Alvero, A Fox, C Cardenas, M Pitruzzello, H Chehade. Generation of stable epithelial-mesenchymal hybrid cancer cells with tumorigenic potential.. Cancers (Basel). (, 2023. [DOI | PubMed]
43. X Wang, PJA Eichhorn, JP Thiery. TGF-β, EMT, and resistance to anti-cancer treatment.. Semin Cancer Biol. (, 2023. [DOI | PubMed]
44. ZY Shi, SX Zhang, CH Li, D Fan, Y Xue, ZH Cheng. Differential distribution and prognostic value of CD4(+) T cell subsets before and after radioactive iodine therapy in differentiated thyroid cancer with varied curative outcomes.. Front Immunol. (, 2022. [DOI | PubMed]
45. Q Liu, W Sun, H Zhang. Roles and new insights of macrophages in the tumor microenvironment of thyroid cancer.. Front Pharmacol. (, 2022. [DOI | PubMed]