Front. Oncol.Frontiers in OncologyFront. Oncol.2234-943XFrontiers Media S.A.10.3389/fonc.2025.1641195OncologyOriginal ResearchMulti-disease transcriptomic analysis of sex hormone genes reveals a novel prognostic model for thyroid cancer with breast cancer correlationsQiaoLixue1†LiHao2†YinKeyu3MaRunsheng145ZhangYifei45GuoYue1YinDetao145*1Thyroid Surgery Department, The First Affiliated Hospital of Zhengzhou University, Zhengzhou, China2Department of Hepatobiliary and Pancreatic Surgery, The First Affiliated Hospital of Zhengzhou University, Zhengzhou, China3School of Basic Medical Sciences, Lanzhou University, Lanzhou, Gansu, China4Engineering Research Center of Multidisciplinary Diagnosis and Treatment of Thyroid Cancer of Henan Province, Zhengzhou,, China5Key Medicine Laboratory of Thyroid Cancer of Henan Province, Zhengzhou, China
Edited by: Mamunur Rahaman, University of New South Wales, Australia
Reviewed by: Tomohiro Chiba, Cancer Institute Hospital of Japanese Foundation for Cancer Research, Japan
Hao Wang, Shenzhen University General Hospital, China
*Correspondence: Detao Yin, detaoyin@zzu.edu.cn
†These authors have contributed equally to this work and share first authorship
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.
Background
There is a potential bidirectional pathogenicity between thyroid and breast cancers. The association between sex hormones and two types of malignant tumors has emerged as a topic of intense academic debate in recent years. However, the role of sex hormone metabolism-related genes in thyroid cancer still needs to be further explored.
Methods
We obtained thyroid and breast cancer transcriptome data from the TCGA database and sex hormone metabolism-related gene sets from the MSigDB database, thus screening for sex hormone metabolism-related genes linked to the two malignant tumors. Univariate cox regression analysis was used for the screening of disease-free survival (DFS)-associated genes. The TCGA-THCA patients were classified as two categories via a consistent clustering algorithm, and the differential genes between the two categories were subsequently screened. A sex hormone metabolism-related prognostic model (TBSMRPM) of thyroid cancer versus breast cancer consisting of 10 genes was developed by Cox regression analyses and least absolute shrinkage with selection operator (LASSO) cox regression analysis. Finally, we performed clinicopathological subgroup analyses to analyze the correlation between TBSMRPM and clinical characteristics, immune infiltration, tumor mutation burden (TMB), and chemosensitivity, and verified the expression of TBSMRPM signature genes by qRT-PCR.
Results
We identified 2 clusters correlated with sex hormone metabolism, and screened 10 prognostic differential genes related to thyroid cancer, breast cancer and sex hormone metabolism. After establishing the two risk groups for thyroid cancer originated from TBSMRPM, the results showed that the high-risk group exhibited the shorter DFS (P<0.05). In further clinical stratification analysis, immune infiltration analysis, TMB and drug sensitivity analysis, the two TBSMRPM groups showed significant differences. The qRT-PCR results showed that C2CD4A, CERS1, MMP9, SLC5A1, HORMAD2 were highly expressed in the IHH4, KTC-1, and TPC-1 cell lines, while SLITRK2, ARHGEF37, PLP1, RNF223, and F3 were lowly expressed.
Conclusion
The TBSMRPM established in this study has a certain value for the prognosis of thyroid cancer and contributes to refine clinicians’ treatment protocols.
Thyroid cancer (THCA) and breast cancer (BRCA), as common hormone-dependent malignancies, their bidirectional pathogenic association and the regulatory role of sex hormones therein have become a focus of oncology research. Epidemiological data clearly demonstrate a significant interaction between these two organ malignancies: Nielsen et al. confirmed that individuals with a history of thyroid cancer have a 1.32-fold higher risk of developing breast cancer compared to the general population, while breast cancer patients have a 1.55-fold increased risk of subsequent thyroid cancer (1). This bidirectional risk suggests that, as hormone-dependent organs, the malignant progression of the thyroid and mammary gland is not an isolated event, and their development is influenced by a complex regulatory network (1–6). Currently, research on their association has involved multiple dimensions, including genetic alterations, hypothalamic-pituitary axis regulation, metabolic abnormalities (e.g., diabetes, obesity), and surveillance bias. Among these, hormone-related mechanisms are particularly critical—various hormone-related factors such as thyroid hormones, sex hormones (with estrogen as the core), endocrine-disrupting chemicals, and adipokines have been confirmed to be involved in the development of both tumors (2). Sex hormones, in particular, have been identified as core regulatory factors due to their direct role in regulating cell proliferation, differentiation, and malignant transformation in thyroid and breast tissues (2). Therefore, systematically elucidating the roles and mechanisms of sex hormone metabolism-related genes associated with both cancers is of great significance for clarifying cross-cancer associations and optimizing clinical diagnosis and treatment.
The metabolic balance of sex hormones is precisely regulated by a network of genes encoding metabolic enzymes and receptors, collectively referred to as sex hormone metabolism-related genes (SMRGs). To date, the roles of SMRGs in breast cancer and thyroid cancer have been partially elucidated. In the field of breast cancer, the clinical value of SMRGs has been verified: for example, aromatase encoded by CYP19A1 is a key rate-limiting enzyme in estrogen synthesis, and its inhibitors have become standard therapeutic agents for ER-positive breast cancer by reducing estrogen production (7); both two-sample Mendelian randomization studies and epidemiological data have shown that elevated levels of total testosterone and bioavailable testosterone increase the risk of ER+ breast cancer, with this association being more pronounced in postmenopausal women (8, 9). In thyroid cancer research, the roles of SMRGs have also been gradually revealed, with the regulatory mechanisms of estrogen and its receptors receiving the most attention: estrogen can enhance the proliferation, migration, and invasion abilities of papillary thyroid carcinoma (PTC) through the ERα/KRT19 signaling axis (10); ERβ is highly expressed in PTC stem cells (PTCSCs), and its knockout can reduce the expression of stemness-related factors, decrease the ALDH+ cell population, and inhibit tumor sphere formation and growth (11); in addition, ERα36, GRP78, and GRP94 are upregulated in primary PTC tissues, and their expression levels can directly affect the malignant phenotype of PTC-derived BCPAP cells (12). Beyond estrogen, androgen receptors (AR), progesterone receptors, and prolactin receptors are also expressed in PTC lesions (13, 14), and AR activation can exert antiproliferative effects (e.g., inducing cellular senescence) in PTC cell models (15, 16). These findings collectively suggest that SMRGs hold promise as novel molecular targets for diagnostic prediction, endocrine therapy, and prognostic evaluation in both cancers.
Despite certain progress in SMRG research, significant limitations remain. Existing studies have shown that Fu et al. identified the COMP gene, which is highly expressed in both breast cancer and thyroid cancer, and confirmed that its overexpression can promote the occurrence and progression of both cancers through the estrogen signaling pathway (17); Jin et al. found that the ratio of ESR1 to ESR can serve as a prognostic marker for predicting survival in female PTC patients and has the potential to be a therapeutic target (18); Zhang et al. constructed a model by screening estrogen-related differential genes associated with THCA and found that this model is closely related to immune infiltration and therapeutic response in THCA (19). However, current research still has three shortcomings: first, the exploration of SMRGs is insufficiently systematic, with most focusing on estrogen metabolism-related genes in PTC, while neglecting the contribution of other sex hormone metabolism genes (such as those involved in androgen and progesterone metabolism) to the development of thyroid cancer; second, most studies are limited to thyroid cancer alone, failing to incorporate its bidirectional risk characteristics with breast cancer, making it difficult to reflect the cross-regulatory association of hormones between the two; third, some studies only focus on the predictive efficacy of SMRGs in PTC, lacking in-depth analysis of their underlying pathogenic mechanisms (e.g., associations with the tumor microenvironment and therapeutic sensitivity). Therefore, there is an urgent need to systematically explore the prognostic efficacy and potential mechanisms of SMRGs associated with both cancers in PTC from a cross-cancer perspective.
To address the above research gaps, this study integrates transcriptomic data of THCA and BRCA from The Cancer Genome Atlas (TCGA) database, and screens SMRGs directly or potentially associated with both cancers through cross-analysis, providing a molecular basis for deciphering cross-cancer metabolic associations. On this basis, consensus clustering and LASSO Cox regression methods are used to construct a thyroid cancer prognostic model (TBSMRPM). The innovation of this study lies in: analyzing the role of SMRGs for the first time from a cross-disease perspective, overcoming the limitations of single-cancer research; and systematically associating the model with immune infiltration, tumor mutational burden (TMB), and drug sensitivity to provide multidimensional references for clinical treatment. Ultimately, it aims to provide new insights for the prognostic evaluation and treatment strategy optimization of thyroid cancer, while laying a foundation for elucidating the role of SMRGs in cross-cancer regulation.
Materials and methodsAcquisition and processing of information
We obtained transcriptomic data for thyroid and breast cancers from the TCGA database and also downloaded clinical and survival information for thyroid cancer. The THCA data lacking DFS survival information were excluded, and the final THCA cohort used for the study included 498 tumor samples and 59 normal samples. The BRCA cohort included 1,118 tumor samples and 113 normal samples. The baseline characteristics table of the patient cohort is shown in Supplementary Table 1. We also logged on to MSigDB (Molecular Signatures Database) to obtain 401 sex hormone metabolism-related genes (20). This research was authorized by the Ethics Committee of the First Affiliated Hospital of Zhengzhou University [approval number: 2020-KY-0075-002]. The overall flowchart of this study is shown in Figure 1.
Sex hormone-related gene analysis workflow.
Flowchart depicting a research process. It starts with the TCGA-BRCA, TCGA-THCA cohorts, and sex hormone metabolism-related genes. These lead to identifying sex hormone-related genes linked to thyroid or breast cancer. The process involves screening prognosis-related genes, clinical subtyping using consensus clustering, and screening SHMRGs from subtypes. Subsequent steps include the TBSMRPM identification, nomogram analysis, immune infiltration analysis, GSEA, TMB analysis, and drug sensitivity prediction. It concludes with verifying key gene expression using qPCR.Consistent clustering based on gene features
We firstly extracted the different genes of thyroid cancer and breast cancer from the acquired transcriptome data, and then intersected with the sex hormone metabolism-related genes to extract the sex hormone metabolism-related genes associated with thyroid cancer or breast cancer. We further obtained thyroid cancer prognosis-related genes via univariate Cox analysis. By utilizing the “ConsensusClusterPlus” R package, we accomplished consensus clustering on the basis of genetic characteristics for the THCA patients (patients were classified as two subclasses, C1 and C2) (21). The specific steps are as follows: First, z-score standardization was performed on the gene expression profiles to eliminate differences in dimensions. The key parameters were set as follows: the number of repeated sampling (reps) was 500, with 80% of samples randomly selected in each sampling, while all features (pFeature=1) were retained for clustering analysis to ensure the robustness of results against data perturbations. The “km” (k-means clustering) algorithm was used for clustering, with Euclidean distance as the distance measurement method, and the range of potential clustering numbers evaluated was k=2 to k=10. The optimal number of clusters was determined through cumulative distribution function (CDF) curve analysis. And we plotted KM curves to determine the DFS difference between the clusters. After that, we also assessed the differences in immune infiltration between clusters and evaluated the stromal and immune scores utilizing the ESTIMATE algorithm.
Identification of differential genes and establishment of TBSMRPM
We identified differential genes in different clusters using the “Deseq2” package (22). A total of 498 thyroid cancer (THCA) patients were randomly divided into a training set and a test set at a ratio of 5:5, with 249 patients in each of the training set and the test set. We used the training set to screen the optimal differential genes by univariate Cox, LASSO cox and stepwise cox regression analyses, respectively, in which the “glmnet” R package was called (23). The median risk scores established by multivariate Cox regression was utilized for the differentiating of the risk groups. We also plotted risk factor linkage charts and KM curves to show the differences in population proportions, survival (DFS) status, and gene expression across risk groups. The sensitivity and specificity of the above differential genomes are reflected by receiver operating characteristic curves (ROCs).
Exploration of model-clinical correlations
To explore the correlation between TBSMRPM and clinical features, we plotted box line plots of the relationship between risk scores and age, gender, stage, extra-glandular invasion, T stage, N stage, M stage, and tumor burden. Stratified analysis further demonstrated the variability of the risk groups within each clinical subgroup.
Modeling and evaluation of the clinical factor-related nomogram
Our Cox regression analysis included TBSMRPM and clinical factors. The univariate Cox regression analysis was utilized for the screening of model components (24). Furthermore, we plotted ROC curves, calibration curves, and decision curves to fully assess the predictive efficacy, calibration, and clinical utility of the model, respectively.
Immunocorrelation analysis and gene set enrichment analysis
We assessed the correlation of immune cells with risk scores by means of seven algorithms (25–30). The immune checkpoints, immune pathways and immune cells were evaluated between the two risk groups (31). The ESTIMATE algorithm was used for the calculating of immune and stromal scores (32). In addition, changes in pathway activity between different risk groups were analyzed by GSEA (33).
Tumor mutational burden
TMB is the counting of somatic mutation sites in the tumor genome, usually described as mutations per million bases (mut/Mb). This indicator may provide some indication of a tumor’s ability to generate neoantigens and predict the efficacy of tumor immunotherapy (34). After integrating the data, we analyzed the somatic mutations associated with risk groupings using the maftools R package (35).
Drug sensitivity analysis
We screened potential therapeutic agents by correlating drug sensitivity with risk subgroups and invoked the “oncoPredict” R package to complete the calculating of drug half-maximal inhibitory concentration (IC50) (36).
Procedures of quantitative real-time PCR
Our research purchased the human papillary thyroid cancer cell lines (IHH4, KTC-1, TPC-1) and normal thyroid cell line (Nthy ori-3-1) from the Shanghai Cell Biochemical Institute (Shanghai, China). These cell lines were cultured in RPMI-1640 medium (Gemini) containing 10% fetal bovine serum (Gemini). Then the medium was placed in an incubator at 5% CO2 and 37°C for further incubation. The extraction of total RNA was performed utilizing TRIzol reagent (Invitrogen). After that, the Prime Script RT reagent kit with gDNA Eraser (Takara) was utilized for the reverse transcription process. We referred to the instructions to complete the PCR reaction by using the SYBR Green Detection kit (Takara) and utilized the 2-ΔΔCt method to complete the calculation of relative gene expression levels. Detailed primer sequences are shown in Supplementary Table 4.
Statistical analysis
All data analyses were finished by R software (version 4.3.1) and GraphPad Prism (version 8.0.2). Differences between the groups were compared using the Wilcoxon test. Kaplan-Meier curves and log-rank tests were utilized for assessment of the DFS discrepancy between the groups. P<0.05 was considered statistically significant.
ResultsConstruction of the consistent clustering
By screening tumor and normal tissue differential genes from the TCGA-THCA cohort and TCGA-BRCA cohort, and then taking the intersection with sex hormone metabolism-related genes, we obtained a total of 159 sex hormone-related genes associated with thyroid or breast cancer (Supplementary Table 2, Supplementary Figure 1C). We then extracted 26 differentially expressed genes associated with DFS by univariate Cox regression. According to the gene characteristics, the consistent clustering results showed that the curves showed a clear flattening trend and the best clustering effect when k = 2 (Figures 2A, B,Supplementary Figures 1D, E). KM curves of the clusters also demonstrated marked differences (Figure 2C, Supplementary Table 3). Enrichment scores for immune cells were assessed in the two clusters using ssGSEA analysis (Figure 2D), which showed that the cells with significant immune infiltration in the two clusters included T helper cells, aDCs, and so on. And the result demonstrated between-cluster variability in infiltration of all immune cells (P<0.05). The immune pathway result illustrated that the differential genes of the two clusters were clearly correlated with pathways of MHC class-I, HLA, and parainflammation (Figure 2E). The enrichment of the 15 immunization pathways also differed significantly across clusters (P<0.05). C1 corresponded to higher stromal and immune scores based on ESTIMATE results (Figure 2F). In addition, the immune checkpoints that were significantly expressed in C2 over C1 included TGFB1, KDR, and CD27, and those that were significantly expressed in C1 over C2 included CD274, VTCN1, and LGALS9 (Figure 2G).
Construction of the consensus clustering. (A) The result of consensus cumulative distribution function. (B) Determination of the ideal number of clusters. (C) K-M curves with between-cluster differences. (D) The ssGSEA scores of 16 kinds of immune cells. (E) Boxplot presenting significant differences in immune pathway enrichment between the two clusters. (F) Correlations between the tumor microenvironment score and the two clusters. (G) The presence of differentially expressed immune checkpoints. ns, no significance, *P<0.05; **P<0.01, ***P<0.001.
Graphical summary of clustering analysis with seven panels labeled A to G. Panel A shows a consensus cumulative distribution function plot with colored lines for clusters one to ten. Panel B is a consensus matrix heatmap for k equals two. Panel C displays a Kaplan-Meier survival plot comparing clusters C1 and C2, with statistical significance. Panels D and E are box plots showing immune infiltration across various immune cells and pathways for clusters C1 and C2. Panel F includes a box plot for immune scores, and panel G shows expression levels of various genes, all compared between clusters C1 and C2.Identification and evaluation of the TBSMRPM
After consistent clustering, we performed differential gene analysis for both clusters and screened a total of 838 differential genes (Supplementary Table 4). Univariate cox regression analysis, LASSO cox regression analysis, and stepwise cox regression analysis was used to screen for the prognostically strong related genes and to fit the optimal model. Ultimately, we screened out a total of 10 differential genes correlated with DFS to establish the TBSMRPM (Figure 3C). The formula is as follows: Risk score = (-5.548expression of SLITRK2) + (-0.450* expression of ARHGEF37) + (-0.296* expression of F3) + (0.240* expression of MMP9) + (0.313* expression of C2CD4A) + (0.449* expression of CERS1) + (0.537* expression of PLP1) + (0.635* expression of RNF223) + (0.641* expression of SLC5A1) + (1.143* expression of HORMAD2). A total of 498 THCA patients were randomly divided into a training set and a test set at a ratio of 5:5, with 249 patients in each of the training set and the test set. The expression distributions, survival status and risk scores of the 10 genes are displayed in Supplementary Figure 2A. The KM curves of the DFS status of the two groups in the training, test, and total sets showed significant differences (Supplementary Figure 2B, Supplementary Table 3), and the ROC curves demonstrated the good predictive efficacy of the model (Supplementary Figure 2C) (33), the area under the ROC curves (AUCs) of the training-focused model to predict 1-, 3-, and 5-year DFS were 0.914, 0.859, and 0.772, respectively, and in the validation-focused model, they were 0.750, 0.795, and 0.770.
Establishment of sex hormone metabolism-related genes prognostic model. (A, B) Parameterization for LASSO analysis. (C) Forest plots for stepwise cox analysis.
Panel A shows a graph of coefficients versus Log Lambda, illustrating how coefficients vary across lambda values. Panel B displays a plot of partial likelihood deviance against Log Lambda, with red dots indicating data points. Panel C is a table listing variables with hazard ratios, confidence intervals, P-values, and coefficients, alongside a forest plot for visual comparison.Correlation of TBSMRPM with clinical features
To determine the value of the model in depth, our research evaluated the correlation of TBSMRPM and important clinical features and found that the risk score was significantly correlated with stage, extra-glandular invasion, T stage, N stage, and tumor burden (Figure 4). Then, our study performed stratified analyses of clinical factors, age<60 years, age≥60 years, female, male, absence of extra-glandular invasion, presence of extra-glandular invasion, clinical stage I/II, clinical stage III/IV, T1/T2, T3/T4, N0, N1, M0, M1, absence of tumor burden, and presence of tumor burden, and other subgroup analyses showed significant variability in the DFS status-KM curves among different risk groups (Figure 5, Supplementary Table 3). These results suggest that TBSMRPM has excellent anticipation ability.
Risk scores of the patients are classified by common clinical features.
Box plots show risk scores across different clinical features: A) Age below and above 60, B) Female vs. male, C) Stage I+II vs. III+IV, D) Tumor stages T1+T2 vs. T3+T4, E) Lymph node involvement N0 vs. N1, F) Metastasis M0 vs. M1, G) No extrathyroidal extension vs. extension, H) Tumor-free vs. with tumor. Significant differences are marked with p-values.
Subgroup KM analysis on the basis of various clinical features.
Survival probability graphs show various medical analysis results. Each graph represents subgroups categorized by age, gender, cancer stage, tumor extension, and node involvement. The x-axis shows time in years, and the y-axis shows survival probability. High and low-risk groups are displayed with distinct lines. Significant differences are indicated by log-rank p-values.Validation and application of the prognostic modelIntegration of clinical factors and evaluation of nomogram performance
Our analysis results indicated that stage, T stage, N stage, M stage, and risk score were prognosis-related factors (Figure 6A). The risk score was established as an independent predictor of prognosis for DFS by multivariate Cox regression analysis (Figure 6B). Subsequently, we incorporated stage, T stage, N stage, M stage, and risk score into the model to predict DFS status of thyroid cancer patients (Figure 6C). ROC results indicated that the AUCs of the nomogram for predicting DFS status at 1, 3, and 5 years were 0.807 (Figure 6D), 0.776 (Figure 6E), and 0.752 (Figure 6F), respectively. Both decision curves (Figures 6G, I) and calibration curves (Figures 6J–L) indicated the ideality of the predictive modeling.
A nomogram model with clinical factors. (A, B) Univariate and multivariate Cox analysis results. (C) The nomogram for predicting DFS of thyroid cancer in total set. ROCs (D-F), decision curves (G-I) and calibration curves (J-L) for the nomogram and other features used to predict DFS.
Composite image showing statistical analyses and model evaluations across different panels. Panel A presents a forest plot with hazard ratios for various variables. Panel B shows a similar plot focusing on different variables. Panel C illustrates a nomogram for predicted probabilities. Panels D, E, and F feature ROC curves, indicating model performance. Panels G, H, and I display decision curve analyses depicting net benefits across risk thresholds. Panels J, K, and L compare predicted probabilities of survival with actual survival proportions, using calibration plots with confidence intervals.
The results of the seven immune algorithms indicated that the risk scores were negatively correlated with cells such as endothelial cells, CD8 T cells, and hemapoietic stem cells (coefficient < -0.3), and positively correlated with NKT cells, sebocytes, and macrophages M0 relationship (coefficient > 0.2) (Supplementary Figure 3). The ESTIMATE result indicated that the high-risk group corresponded to higher ESTIMATE scores (Figure 7A), which suggested a higher level of immune infiltration. The immune pathway analysis showed that the high-risk group was significantly enriched in Parainflammation, MHC class-I and HLA pathways (Figure 7B). Meanwhile, 12 immune pathways such as co-inhibition, APC co-stimulation, and CCR were significantly different between the TBSMRPM groups. The CD4 T cells, CD56 dim tural killer cells, etc., also showed significant differences between different risk groups, and a total of 26 immune cell infiltrations were significantly different between the two risk groups (Figure 7C). Differential expression analysis of immune checkpoints in different risk groups showed that the immune checkpoints more significantly expressed in the low-risk group included CD160, KDR, CD27, and so on. The more significantly expressed immune checkpoints in the high-risk group included LAG3, TGFBR1, LGALS9 and so on (Figure 7D).
Immunocorrelation results of the two risk groups. (A) Correlations between the tumor microenvironment score and different risk groups. (B) Boxplot for differential analysis of immune pathway enrichment. (C) Boxplot presenting the significant differences in 28 types of immune cell enrichment between different risk groups. (D) Boxplot for differential analysis of immune checkpoint expression. ns, no significance, *P<0.05; **P<0.01, ***P<0.001.
Box plot panels labeled A, B, C, and D compare immune infiltration and gene expression levels. Red boxes indicate high risk; blue boxes indicate low risk. Significant differences are marked by asterisks and “ns” for non-significant. The data categories are stromal score, immune score, and various immune cell and gene expression metrics.Correlation with tumor immune microenvironment and functional pathway enrichment
To reveal gene expression patterns associated with specific biological processes, pathways or functions, our study performed GSEA analysis and it showed that apical junction, estrogen response late, antigen processing and presentation signaling pathways and other pathways had important roles in the high-risk group (Figures 8A–H), which are important in the development and progression of malignant tumors. While signaling pathways such as fatty acid metabolism, glycerolipid metabolism, phosphatidylinositol signaling system, and steroid hormone biosynthesis were significant in the low-risk group (Figures 8I–L).
An in-depth exploration of the role of molecules. Significantly enriched signaling pathways in high-risk group (A-G) and low-risk group (I-L) in accordance to GSEA.
Twelve line graphs display running enrichment scores for various biological pathways. Each panel (A-L) represents different pathways such as HALLMARK and KEGG categories. The x-axes show the rank in ordered datasets, and the y-axes show enrichment scores. Colored bars and lines depict different genetic expressions and pathways, highlighting patterns and peak scores across the datasets.Association with tumor mutational burden and prognostic implications
The waterfall plot indicated that the rank order of frequency of somatic mutations in the low-risk group was BRAF (42%), NRAS (13%), TG (5%), HRAS (4%), and other genes were mutated in less than 4% of the genes (Figure 9A). The rank order of frequency of somatic mutations occurring in the high-risk group was BRAF (77%), TTN (5%), NRAS (4%), and other genes had a mutation frequency of less than 4% (Figure 9B). The somatic mutations in both groups were predominantly missense mutations (Supplementary Figures 4A, B). The box line plot visualized the significant variability of somatic mutations between the risk groups, and the high-risk group had higher TMB values (Figure 9C). Distinguishing the high and low TMB groups by the median somatic mutation scores of all patients showed a significant difference in KM analysis of DFS status between the two groups (Figure 9D, Supplementary Table 3). The KM curve analysis of the combined TMB group and risk group still showed a significant difference between the four groups (Figure 9E, Supplementary Table 3). In addition, Supplementary Figure 4C demonstrates that risk score and TMB were significantly correlated (R=0.25, P<0.0001). Our findings revealed that poor prognosis in the high-risk group may be associated with higher TMB.
Differential analysis of TMB. The mutations landscape of low-risk group (A) and high-risk group (B). (C) Relevance of TMB to risk groups. (D) Differences in DFS of TMB groups. (E) Comparison of KM curves for the combined TMB and risk groups.
Panel A and B: Oncoprints showing gene mutations across two sample groups, highlighting the percentage of samples altered and specific gene mutations. Panel C: Box plot comparing tumor mutation burden between high and low groups, showing significant difference (p=5.4e-06). Panel D: Survival plot stratified by high and low tumor mutation burden, with a log-rank p-value of 0.0067. Panel E: Kaplan-Meier survival curves for combinations of tumor mutation burden and risk levels, showing significant survival differences (p<0.0001).Prediction of drug sensitivity based on risk stratification
The drug sensitivity analysis revealed that Axitinib (Figure 10A), Cyclophosphamide (Figure 10B), Lapatinib (Figure 10D), fulvestrant (Figure 10G), nilotinib (Figure 10H) were more sensitive in the low-risk group, and Erlotinib (Figure 10C), ERK-6640 (Figure 10E), AZD7762 (Figure 10F), and Sapitinib (Figure 10I) were more sensitive in the high-risk group.
Drug sensitivity in different risk groups.
Violin plots labeled A to I compare estimated IC50 values for various drugs between high and low groups. Each plot shows red (high) and blue (low) distributions with associated p-values, indicating statistical significance in the differences between the groups for each drug.PCR for determining gene expression
We examined the expression of 10 genes in IHH4, KTC-1, TPC-1 human papillary thyroid carcinoma cell lines and human thyroid normal cell line Nthy ori-3-1 by PCR analysis. The results showed that C2CD4A, CERS1, MMP9, SLC5A1, and HORMAD2 were highly expressed in IHH4, KTC-1, and TPC-1 cell lines, and SLITRK2, ARHGEF37, PLP1, RNF223, and F3 were lowly expressed in IHH4, KTC-1, and TPC-1 cell lines (Figure 11). The associations between the 10 genes and sex hormone metabolism are detailed in Supplementary Table 8.
Expression of 10 genes in IHH4, KTC-1, TPC-1 human papillary thyroid cancer cell lines and human thyroid normal cell line Nthy ori-3-1. * indicates P<0.05, ** indicates P < 0.01, *** indicates P < 0.001; # indicates P<0.05, ## indicates P<0.01, ### indicates P<0.001, & indicates P<0.05, && indicates P<0.01, &&& indicates P<0.001.
Bar graphs illustrate the relative mRNA expression levels of various genes (CERS1, C2CD4A, HORMAD2, SLC5A1, MMP9, PLP1, RNF223, SLITRK2, ARHGEF37, F3) across different cell types (NTHY, IHH4, KTC-1, TPC-1). Each graph shows expression levels with statistical significance indicated by asterisks and symbols above the bars.Discussion
Numerous studies have shown that there is a potential, bidirectional pathogenic relationship between thyroid cancer and breast cancer, and that the development of one malignant tumor may increase the risk of developing the other (37–39). Of interest, although both malignancies are highly prevalent in the female population, the simultaneous occurrence of both diseases is more common in men (40). This suggests that sex hormones may serve as a vital role in the tumorigenesis and progression of thyroid cancer. However, studies on the prognostic relationship between sex hormones and thyroid cancer are still very limited (2, 41). In our research, we constructed a prognostic model for the first time on the basis of sex hormone metabolism-related genes in thyroid cancer and breast cancer, which was able to better predict the prognostic (DFS) profile of thyroid cancer patients.
We first screened for genes correlated with sex hormone metabolism. Notably, considering the potential bidirectional pathogenic relationship between thyroid and breast cancers, we could not ignore the possible indirect effects of genes with altered expression only in a single malignancy on the other tumor. Thus, we included genes related to sex hormone metabolism associated with thyroid or breast cancer together in our study. This idea is different from the study of Duan et al. (42). Subsequently, we classified thyroid cancer patients into two subtypes using consensus clustering. Since the core basis of consensus clustering lies in the similarity of expression patterns of “sex hormone metabolism-related genes”, the algorithm automatically clusters samples with more similar expression patterns into the same group. Therefore, the division of C1 and C2 inherently implies “distinguishable differences in the overall expression profiles of sex hormone-related genes”—this serves as the premise for the stable output of 2 subtypes by clustering. In addition, we utilized existing analyses (pathway enrichment, immune microenvironment, and prognosis association) to cross-validate the expression differences of sex hormone-related genes at the functional level. Furthermore, the significant differences in prognosis and immune infiltration between different subtypes can provide references for clinical staging. To further screen the more desirable prognostic genes and fit an excellent prognostic model, we performed univariate Cox, LASSO Cox and stepwise cox regression analyses after differential gene analysis of subtypes. The ROC results of both training and validation sets suggested a better predictive performance of the model. And the significant survival difference results in the clinical correlation study indicated the rationality of the TBSMRPM risk grouping based on TBSMRPM. We further established a nomogram model in accordance with clinical factors for the prognostic prediction.
In addition, understanding the tumor immune microenvironment contributes to the diagnosis and treatment of thyroid cancer (43). Various cells in the tumor microenvironment such as: endothelial cells, fibroblasts, and infiltrating immune cells serve as vital roles in the processes of tumorigenesis, invasion, and metastasis. In our results, activated dendritic cells, activated CD8+ T cells, Th1 cells, and macrophages were significantly infiltrated in the high-risk group, which may be linked to the adverse prognosis. Relevant research has shown that dendritic cells may result in the tumor progression by affecting the tumor angiogenesis and causing immune dysregulation (44), whereas significant infiltration of activated CD8+ T cells, Th1 cells suggests a higher risk of recurrence (45–47), and significant infiltration of macrophages in thyroid cancer is strongly linked to tumor invasion and short survival (48). The immune pathway analysis illustrated that the risk grouping was strongly associated with the MHC class I pathway. Angell et al. demonstrated that the absence of MHC class I expression was a vital mechanism of immune escape in thyroid malignancies, which highly overlapped with our findings (49). The immune checkpoint results demonstrated the strong correlation of risk groupings with some immune checkpoints. The study by Park et al. also confirmed that these differential genes associated with immune escape signaling were overexpressed in patients with advanced thyroid cancer, which could help to predict the recurrence of advanced thyroid cancer (47). In conclusion, our immune correlation analysis demonstrated the variability of immune function between different subtypes and risk groups in the consistent clustering. The different risk subgroups should be treated with different therapeutic strategies, and the rational formulation of immunotherapy can help to improve the efficacy of the treatment and prevent the development of drug resistance (50).
In the subsequent GSEA analysis, we found that signaling pathways such as fatty acid metabolism, glycerolipid metabolism, phosphatidylinositol signaling system, and steroid hormone biosynthesis were significantly enriched in the low-risk group. It is well known that sex hormones are important members of the steroid family, and their synthesis and metabolism are inextricably linked to these signaling pathways. Sun et al. found that increased fatty acid metabolism was linked to the decreased infiltration of T- and B-cells in male breast cancer and other cancer types, and that targeting fatty acid metabolism pathways may alleviate the immunosuppressive microenvironment in a variety of cancers (51). Lu et al. demonstrated the importance of lipid metabolism in thyroid cancer (52). Liu et al. experimentally verified that pyruvate carboxylase is closely associated with tumor aggressiveness in thyroid cancer by stimulating fatty acid synthesis (53). The phosphatidylinositol (PI) signaling system is a key cellular signaling pathway concerning the cellular processes such as cell growth, differentiation, survival, and intracellular communication, and dysregulation of the related signaling pathway is associated with thyroid cancer and breast cancer (54, 55). Furthermore, a small number of studies examining the association of other sex hormones with thyroid cancer, in addition to estrogen-related studies, have also yielded positive findings (56, 57).
After that, we performed somatic mutation and drug sensitivity analysis. And we found the high-risk group had a significantly higher frequency of BRAF gene mutations. The high-risk group had a higher TMB value and corresponded to the poor prognosis. Studies have confirmed the importance of the detection of BRAF gene mutations in the diagnosis, treatment and determination of prognosis of thyroid cancer (58). This is consistent with the results of our study. There are fewer studies related to other genes with higher mutation frequencies such as NRAS in our findings (59), which found that NRAS mutations were linked to the high risk of distant metastasis of thyroid cancer (60, 61). These findings suggest, to some extent, the potential role of this gene in thyroid cancer. Furthermore, we gained more insight into the correlation of drug sensitivity between different risk groups.
Ultimately, we screened 10 ideal prognosis-related genes. Among them, MMP9 has been shown to promote thyroid cancer incidence and progression in a large number of literatures (62, 63), which is consistent with our findings. HORMAD2 has been shown to inhibit the incidence and progression of thyroid cancer and its hypermethylation promotes the progression of thyroid cancer in a study by Lin et al. (64). Han et al. mentioned that the mRNA level of C2CD4A was dysregulated in the lung cancer (65). Rong et al. confirmed that C2CD4A was highly expressed in colorectal cancer tissues and contributed to tumor growth by inhibiting the P53 signaling pathway (66). Chen et al. verified that CERS1 was highly expressed in colorectal cancer by PCR (67). Xu et al. found that CERS1 was significantly downregulated in non-small-cell lung cancer cell lines and brain metastatic tissue and its upregulation corresponded to better prognosis (68). Lei et al. proposed that RNF223 could promote pancreatic cancer growth and migration, and identified potential protein targets and metabolism-related pathways (69). Zhang et al. demonstrated that ARHGEF37 promotes hepatocellular carcinoma outgrowth and metastasis through activation of Cdc42 (70). SLC5A1 has been confirmed to be highly expressed in many malignant tumors by several studies, including our research, and a study on pancreatic cancer demonstrated that SLC5A1 regulates cancer cell growth through AMPK/mTOR signaling (71). However, studies on the association of these genes with thyroid cancer have rarely been reported and need to be further explored.
There are some limitations to this study. First, due to limited sample acquisition, we performed PCR validation only in papillary thyroid carcinoma. Second, the model established in this study lacks validation in the external cohorts. Third, the actual role of the screened genes in thyroid cancer lacks sufficient experimental validation. In conclusion, by screening differential genes and establishing a prognostic model, we mined potential novel biomarkers related to thyroid and breast cancers, which are highly suggestive for prognostic prediction of thyroid cancer. However, the actual clinical value of the model needs to be further validated in large cohort studies, and the actual roles and specific mechanisms of these potential biomarkers in malignant tumors need to be further determined by a large number of in vivo and ex vivo experiments.
Data availability statement
The original contributions presented in the study are included in the article/Supplementary Material. Further inquiries can be directed to the corresponding author.
Ethics statement
Ethical approval was not required for the studies on humans in accordance with the local legislation and institutional requirements because only commercially available established cell lines were used.
Author contributions
LQ: Writing – original draft, Writing – review & editing. HL: Writing – review & editing, Writing – original draft. KY: Writing – original draft, Writing – review & editing. RM: Writing – review & editing, Writing – original draft. YZ: Writing – original draft, Writing – review & editing. YG: Writing – original draft, Writing – review & editing. DY: Writing – review & editing, Writing – original draft.
Funding
The author(s) declare that no financial support was received for the research and/or publication of this article.
Acknowledgments
We are very grateful to the staff in TCGA Program for their kind work in data collection and delivery. This study was funded by the General Project of Natural Science Foundation of Henan Province (No.222300420568), the Key Medical Science and Technology Project of Henan Province (No. SBGJ202101014), the Major Scientific Research Projects of Traditional Chinese Medicine in Henan Province (No.20-21ZYZD14) and the Cultivation of Young and Middle-aged Health Science and Technology Innovation Leading Talents in Henan Province (No. YXKC2020015). Thanks to the above funds for supporting this research.
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.
Generative AI statement
The author(s) declare that no Generative AI was used in the creation of this manuscript.
Any alternative text (alt text) provided alongside figures in this article has been generated by Frontiers with the support of artificial intelligence and reasonable efforts have been made to ensure accuracy, including review by the authors wherever possible. If you identify any issues, please contact us.
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.
Supplementary material
The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fonc.2025.1641195/full#supplementary-material
ReferencesNielsenSMWhiteMGHongSAschebrook-KilfoyBKaplanELAngelosP. The breast-thyroid cancer link: A systematic review and meta-analysis. Cancer Epidemiol Biomarkers Prev. (2016) 25:231–8. doi: 10.1158/1055-9965.EPI-15-0833, PMID: 26908594HaladaSCasado-MedranoVBaranJALeeJChinmayPBauerAJ. Hormonal crosstalk between thyroid and breast cancer. Endocrinology. (2022) 163. doi: 10.1210/endocr/bqac075.Cited in: Pubmed, PMID: 35587175BolfELSpragueBLCarrFE. A Linkage Between Thyroid and Breast Cancer: A Common Etiology? Cancer Epidemiol Biomarkers Prev. (2019) 28:643–9. doi: 10.1158/1055-9965.EPI-18-0877, PMID: 30541751BakosBKissAÁrvaiKSziliBDeák-KocsisBTobiásB. Co-occurrence of thyroid and breast cancer is associated with an increased oncogenic SNP burden. BMC Cancer. (2021) 21:706. doi: 10.1186/s12885-021-08377-4, PMID: 34130653AvgerinosKISpyrouNMantzorosCSDalamagaM. Obesity and cancer risk: Emerging biological mechanisms and perspectives. Metabolism. (2019) 92:121–35. doi: 10.1016/j.metabol.2018.11.001, PMID: 30445141YinD-THeHYuKXieJLeiMMaR. The association between thyroid cancer and insulin resistance, metabolic syndrome and its components: A systematic review and meta-analysis. Int J Surg. (2018) 57:66–75. doi: 10.1016/j.ijsu.2018.07.013, PMID: 30081182HamadehISPatelJNRusinSTanAR. Personalizing aromatase inhibitor therapy in patients with breast cancer. Cancer Treat Rev. (2018) 70:47–55. doi: 10.1016/j.ctrv.2018.07.014, PMID: 30086432NounuAKarSPReltonCLRichmondRC. Sex steroid hormones and risk of breast cancer: a two-sample Mendelian randomization study. Breast Cancer Res. (2022) 24:66. doi: 10.1186/s13058-022-01553-9, PMID: 36209141CoradiniDOrianaS. Impact of sex hormones dysregulation and adiposity on the outcome of postmenopausal breast cancer patients. Clin Obes. (2021) 11:e12423. doi: 10.1111/cob.12423, PMID: 33135396SongZMWangYDChaiFZhangJLvSWangJX. Estrogen enhances the proliferation, migration, and invasion of papillary thyroid carcinoma via the ERα/KRT19 signaling axis. J Endocrinol Invest. (2025) 48:653–70. doi: 10.1007/s40618-024-02473-5, PMID: 39453570LiMChaiH-FPengFMengY-TZhangL-ZZhangL. Estrogen receptor β upregulated by lncRNA-H19 to promote cancer stem-like properties in papillary thyroid carcinoma. Cell Death Dis. (2018) 9:1120. doi: 10.1038/s41419-018-1077-9, PMID: 30389909DaiY-JQiuY-BJiangRXuMLiaoL-YChenGG. Concomitant high expression of ERα36, GRP78 and GRP94 is associated with aggressive papillary thyroid cancer behavior. Cell Oncol (Dordr). (2018) 41:269–82. doi: 10.1007/s13402-017-0368-y, PMID: 29368272AhnHYSongR-YAhnHSKimHS. Expression of estrogen and progesterone receptors in papillary thyroid carcinoma in korea. Cancer Res Treat. (2021) 53:1204–12. doi: 10.4143/crt.2020.1201, PMID: 33592140CostaPCatarinoALSilvaFSobrinhoLGBugalhoMJ. Expression of prolactin receptor and prolactin in normal and Malignant thyroid: a tissue microarray study. Endocr Pathol. (2006) 17:377–86. doi: 10.1007/s12022-006-0009-x, PMID: 17525486BanuKSGovindarajuluPAruldhasMM. Testosterone and estradiol have specific differential modulatory effect on the proliferation of human thyroid papillary and follicular carcinoma cell lines independent of TSH action. Endocr Pathol. (2001) 12:315–27. doi: 10.1385/EP:12:3:315, PMID: 11740053GuptaACarnazzaMJonesMDarzynkiewiczZHalickaDO’ConnellT. Androgen receptor activation induces senescence in thyroid cancer cells. Cancers (Basel). (2023) 15. doi: 10.3390/cancers15082198, PMID: 37190127FuJHeMWuQZhangXQiXShenK. The clinical and genetic features in patients coexisting primary breast and thyroid cancers. Front Endocrinol (Lausanne). (2023) 14:1136120. doi: 10.3389/fendo.2023.1136120, PMID: 37229458YiJWKimS-JKimJKSeongCYYuHWChaiYJ. Upregulation of the ESR1 gene and ESR ratio (ESR1/ESR2) is associated with a worse prognosis in papillary thyroid carcinoma: the impact of the estrogen receptor α/β Expression on clinical outcomes in papillary thyroid carcinoma patients. Ann Surg Oncol. (2017) 24:3754–62. doi: 10.1245/s10434-017-5780-z, PMID: 28124274ZhangLZhouMGaoXXieYXiaoJLiuT. Estrogen-related genes for thyroid cancer prognosis, immune infiltration, staging, and drug sensitivity. BMC Cancer. (2023) 23:1048. doi: 10.1186/s12885-023-11556-0, PMID: 37907864LiberzonABirgerCThorvaldsdóttirHGhandiMMesirovJPTamayoP. The Molecular Signatures Database (MSigDB) hallmark gene set collection. Cell Syst. (2015) 1:417–25. doi: 10.1016/j.cels.2015.12.004, PMID: 26771021WilkersonMDHayesDN. ConsensusClusterPlus: a class discovery tool with confidence assessments and item tracking. Bioinformatics. (2010) 26:1572–3. doi: 10.1093/bioinformatics/btq170, PMID: 20427518RitchieMEPhipsonBWuDHuYLawCWShiW. limma powers differential expression analyses for RNA-sequencing and microarray studies. Nucleic Acids Res. (2015) 43:e47. doi: 10.1093/nar/gkv007, PMID: 25605792FriedmanJHastieTTibshiraniR. Regularization paths for generalized linear models via coordinate descent. J Stat Software. (2010) 33. doi: 10.18637/jss.v033.i01, PMID: 20808728IasonosASchragDRajGVPanageasKS. How to build and interpret a nomogram for cancer prognosis. J Clin Oncol. (2008) 26:1364–70. doi: 10.1200/JCO.2007.12.9791, PMID: 18323559LiTFanJWangBTraughNChenQLiuJS. TIMER: A web server for comprehensive analysis of tumor-infiltrating immune cells. Cancer Res. (2017) 77:e108–10. doi: 10.1158/0008-5472.CAN-17-0307, PMID: 29092952ChenBKhodadoustMSLiuCLNewmanAMAlizadehAA. Profiling tumor infiltrating immune cells with CIBERSORT. Methods Mol Biol. (2018) 1711:243–59. doi: 10.1007/978-1-4939-7493-1_12, PMID: 29344893PlattnerCFinotelloFRiederD. Deconvoluting tumor-infiltrating immune cells from RNA-seq data using quanTIseq. Methods Enzymol. (2020) 636:261–85. doi: 10.1016/bs.mie.2019.05.056, PMID: 32178821BechtEGiraldoNALacroixLButtardBElarouciNPetitprezF. Estimating the population abundance of tissue-infiltrating immune and stromal cell populations using gene expression. Genome Biol. (2016) 17:218. doi: 10.1186/s13059-016-1070-5, PMID: 27765066HongKCenKChenQDaiYMaiYGuoY. Identification and validation of a novel senescence-related biomarker for thyroid cancer to predict the prognosis and immunotherapy. Front Immunol. (2023) 14:1128390. doi: 10.3389/fimmu.2023.1128390, PMID: 36761753AranDHuZButteAJ. xCell: digitally portraying the tissue cellular heterogeneity landscape. Genome Biol. (2017) 18:220. doi: 10.1186/s13059-017-1349-1, PMID: 29141660CharoentongPFinotelloFAngelovaMMayerCEfremovaMRiederD. Pan-cancer immunogenomic analyses reveal genotype-immunophenotype relationships and predictors of response to checkpoint blockade. Cell Rep. (2017) 18:248–62. doi: 10.1016/j.celrep.2016.12.019, PMID: 28052254YoshiharaKShahmoradgoliMMartínezEVegesnaRKimHTorres-GarciaW. Inferring tumour purity and stromal and immune cell admixture from expression data. Nat Commun. (2013) 4:2612. doi: 10.1038/ncomms3612, PMID: 24113773SubramanianATamayoPMoothaVKMukherjeeSEbertBLGilletteMA. Gene set enrichment analysis: a knowledge-based approach for interpreting genome-wide expression profiles. Proc Natl Acad Sci U S A. (2005) 102:15545–50., PMID: 16199517ChalmersZRConnellyCFFabrizioDGayLAliSMEnnisR. Analysis of 100,000 human cancer genomes reveals the landscape of tumor mutational burden. Genome Med. (2017) 9:34. doi: 10.1186/s13073-017-0424-2, PMID: 28420421MayakondaALinD-CAssenovYPlassCKoefflerHP. Maftools: efficient and comprehensive analysis of somatic variants in cancer. Genome Res. (2018) 28:1747–56. doi: 10.1101/gr.239244.118, PMID: 30341162MaeserDGruenerRFHuangRS. oncoPredict: an R package for predicting in vivo or cancer patient drug response and biomarkers from cell line screening data. Brief Bioinform. (2021) 22. doi: 10.1093/bib/bbab260, PMID: 34260682Van FossenVLWilhelmSMEatonJLMcHenryCR. Association of thyroid, breast and renal cell cancer: a population-based study of the prevalence of second Malignancies. Ann Surg Oncol. (2013) 20:1341–7. doi: 10.1245/s10434-012-2718-3, PMID: 23263698LiSYangJShenYZhaoXZhangLWangB. Clinicopathological features, survival and risk in breast cancer survivors with thyroid cancer: an analysis of the SEER database. BMC Public Health. (2019) 19:1592. doi: 10.1186/s12889-019-7947-y, PMID: 31783815JosephKREdirimanneSEslickGD. The association between breast cancer and thyroid cancer: a meta-analysis. Breast Cancer Res Treat. (2015) 152:173–81. doi: 10.1007/s10549-015-3456-6, PMID: 26058757ZafonCObiolsGMesaJ. Second primary cancer in patients with papillary thyroid carcinoma. Anticancer Res. (2013) 33:337–40.TanHWangSHuangFTongZ. Association between breast cancer and thyroid cancer risk: a two-sample Mendelian randomization study. Front Endocrinol (Lausanne). (2023) 14:1138149. doi: 10.3389/fendo.2023.1138149, PMID: 37288296DuanJLiuCYiJWangY. Shared sex hormone metabolism-related gene prognostic index between breast and endometrial cancers. Front Endocrinol (Lausanne). (2023) 14:1126862. doi: 10.3389/fendo.2023.1126862, PMID: 36742386MenicaliEGuzzettiMMorelliSMorettiSPuxedduE. Immune landscape of thyroid cancers: new insights. Front Endocrinol (Lausanne). (2020) 11:637826. doi: 10.3389/fendo.2020.637826, PMID: 33986723QiuHHuXHeCYuBLiYLiJ. Identification and validation of an individualized prognostic signature of bladder cancer based on seven immune related genes. Front Genet. (2020) 11:12. doi: 10.3389/fgene.2020.00012, PMID: 32117435CunhaLLMarcelloMANonogakiSMorariECSoaresFAVassalloJ. CD8+ tumour-infiltrating lymphocytes and COX2 expression may predict relapse in differentiated thyroid cancer. Clin Endocrinol (Oxf). (2015) 83:246–53. doi: 10.1111/cen.12586, PMID: 25130519ChenZGuoM-LLiY-YYanKLiLShenF. Immune profiling identifies CD8+ T-cell subset signatures as prognostic markers for recurrence in papillary thyroid cancer. Front Immunol. (2022) 13:894919. doi: 10.3389/fimmu.2022.894919, PMID: 36420264ParkS-JKangYEKimJ-HParkJ-LKimS-KBaekS-W. Transcriptomic analysis of papillary thyroid cancer: A focus on immune-subtyping, oncogenic fusion, and recurrence. Clin Exp Otorhinolaryngol. (2022) 15:183–93. doi: 10.21053/ceo.2021.02215, PMID: 35255661YinHTangYGuoYWenS. Immune microenvironment of thyroid cancer. J Cancer. (2020) 11:4884–96. doi: 10.7150/jca.44506, PMID: 32626535AngellTELechnerMGJangJKLoPrestiJSEpsteinAL. MHC class I loss is a frequent mechanism of immune escape in papillary thyroid cancer that is reversed by interferon and selumetinib treatment in vitro. Clin Cancer Res. (2014) 20:6034–44. doi: 10.1158/1078-0432.CCR-14-0879, PMID: 25294906MouldRCvan VlotenJPAuYeungAWKKarimiKBridleBW. Immune responses in the thyroid cancer microenvironment: making immunotherapy a possible mission. Endocr Relat Cancer. (2017) 24:T311–29. doi: 10.1530/ERC-17-0316, PMID: 28912377SunHZhangLWangZGuDZhuMCaiY. Single-cell transcriptome analysis indicates fatty acid metabolism-mediated metastasis and immunosuppression in male breast cancer. Nat Commun. (2023) 14:5590. doi: 10.1038/s41467-023-41318-2, PMID: 37696831LuJZhangYSunMDingCZhangLKongY. Multi-omics analysis of fatty acid metabolism in thyroid carcinoma. Front Oncol. (2021) 11:737127. doi: 10.3389/fonc.2021.737127, PMID: 34976793LiuCZhouXPanYLiuYZhangY. Pyruvate carboxylase promotes thyroid cancer aggressiveness through fatty acid synthesis. BMC Cancer. (2021) 21:722. doi: 10.1186/s12885-021-08499-9, PMID: 34158007SantarpiaLEl-NaggarAKCoteGJMyersJNShermanSI. Phosphatidylinositol 3-kinase/akt and ras/raf-mitogen-activated protein kinase pathway mutations in anaplastic thyroid cancer. J Clin Endocrinol Metab. (2008) 93:278–84. doi: 10.1210/jc.2007-1076, PMID: 17989125AmiraniEHallajzadehJAsemiZMansourniaMAYousefiB. Effects of chitosan and oligochitosans on the phosphatidylinositol 3-kinase-AKT pathway in cancer therapy. Int J Biol Macromol. (2020) 164:456–67. doi: 10.1016/j.ijbiomac.2020.07.137, PMID: 32693135ZhangLJXiongYNilubolNHeMBommareddiSZhuX. Testosterone regulates thyroid cancer progression by modifying tumor suppressor genes and tumor immunity. Carcinogenesis. (2015) 36:420–8. doi: 10.1093/carcin/bgv001, PMID: 25576159O’ConnellTJDadafarinSJonesMRodríguezTGuptaAShinE. Androgen activity is associated with PD-L1 downregulation in thyroid cancer. Front Cell Dev Biol. (2021) 9:663130. doi: 10.3389/fcell.2021.663130, PMID: 34422798CaroniaLMPhayJEShahMH. Role of BRAF in thyroid oncogenesis. Clin Cancer Res. (2011) 17:7511–7. doi: 10.1158/1078-0432.CCR-11-1155, PMID: 21900390YinD-TYuKLuR-QLiXXuJLeiM. Clinicopathological significance of TERT promoter mutation in papillary thyroid carcinomas: a systematic review and meta-analysis. Clin Endocrinol (Oxf). (2016) 85:299–305. doi: 10.1111/cen.13017, PMID: 26732020MeloMGaspar da RochaABatistaRVinagreJMartinsMJCostaG. BRAF, and NRAS in primary thyroid cancer and metastatic disease. J Clin Endocrinol Metab. (2017) 102:1898–907. doi: 10.1210/jc.2016-2785, PMID: 28323937JangEKSongDESimSYKwonHChoiYMJeonMJ. NRAS codon 61 mutation is associated with distant metastasis in patients with follicular thyroid carcinoma. Thyroid. (2014) 24:1275–81. doi: 10.1089/thy.2014.0053, PMID: 24820091GuanHGuoYLiuLYeRLiangWLiH. INAVA promotes aggressiveness of papillary thyroid cancer by upregulating MMP9 expression. Cell Biosci. (2018) 8:26. doi: 10.1186/s13578-018-0224-4, PMID: 29632659OuyangXFengLYaoLXiaoYHuXZhangG. Testicular orphan receptor 4 (TR4) promotes papillary thyroid cancer invasion via activating circ-FNLA/miR-149-5p/MMP9 signaling. Mol Ther Nucleic Acids. (2021) 24:755–67. doi: 10.1016/j.omtn.2021.03.021, PMID: 33996257LinQHouSGuanFLinC. HORMAD2 methylation-mediated epigenetic regulation of gene expression in thyroid cancer. J Cell Mol Med. (2018) 22:4640–52. doi: 10.1111/jcmm.13680, PMID: 30039914HanFQianLZhangYLiuPLiRChenM. C2CD4A/B variants in the predisposition of lung cancer in the Chinese Han population. Funct Integr Genomics. (2022) 22:331–40. doi: 10.1007/s10142-022-00827-x, PMID: 35212842RongZLuoZFuZZhangPLiTZhangJ. The novel circSLC6A6/miR-1265/C2CD4A axis promotes colorectal cancer growth by suppressing p53 signaling pathway. J Exp Clin Cancer Res. (2021) 40:324. doi: 10.1186/s13046-021-02126-y, PMID: 34656159ChenLChenHLiYLiLQiuYRenJ. Endocannabinoid and ceramide levels are altered in patients with colorectal cancer. Oncol Rep. (2015) 34:447–54. doi: 10.3892/or.2015.3973, PMID: 25975960XuYPanJLinYWuYChenYLiH. Ceramide synthase 1 inhibits brain metastasis of non-small cell lung cancer by interacting with USP14 and downregulating the PI3K/AKT/mTOR signaling pathway. Cancers (Basel). (2023) 15. doi: 10.3390/cancers15071994, PMID: 37046655FengLWangJZhangJDiaoJHeLFuC. Comprehensive analysis of E3 ubiquitin ligases reveals ring finger protein 223 as a novel oncogene activated by KLF4 in pancreatic cancer. Front Cell Dev Biol. (2021) 9:738709. doi: 10.3389/fcell.2021.738709, PMID: 34722520ZhangXRenLWuJFengRChenYLiR. ARHGEF37 overexpression promotes extravasation and metastasis of hepatocellular carcinoma via directly activating Cdc42. J Exp Clin Cancer Res. (2022) 41:230. doi: 10.1186/s13046-022-02441-y, PMID: 35869555GaoH-FChenL-YChengC-SChenHMengZ-QChenZ. SLC5A1 promotes growth and proliferation of pancreatic carcinoma via glucose-dependent AMPK/mTOR signaling. Cancer Manag Res. (2019) 11:3171–85. doi: 10.2147/CMAR.S195424, PMID: 31114359