CREM regulates the tumor immune microenvironment and predicts prognosis in thyroid carcinoma

Introduction

Thyroid cancer (THCA) is among the most prevalent malignant tumors of the endocrine system worldwide, with a steadily rising incidence (1), particularly in papillary thyroid carcinoma (PTC), which accounts for over 80% of all THCA cases (2,3). Epidemiological studies highlight a higher incidence in women, and men have a worse prognosis (4,5). PTC is generally characterized by indolent growth yet shows a notable tendency for regional lymph node metastasis, particularly in pediatric populations (6,7). Surgical resection, typically followed by radioactive iodine therapy, is standard, yet challenges in early detection and risk stratification underscore the need for new biomarkers (8,9). Prognostic factors associated with unfavorable outcomes include male sex, tumor size ≥2 cm, distant metastasis, and elevated pre-treatment thyroid-stimulating hormone (TSH) levels (10,11). Genetic mutations and epigenetic modifications are closely associated with thyroid tumorigenesis and progression (12,13). Additionally, the tumor microenvironment (TME) has emerged as a critical component in THCA pathophysiology. Infiltration by tumor-associated macrophages (TAMs) and other immune cells has been associated with enhanced tumor aggressiveness and metastatic potential (14,15). For example, infiltration by TAMs, regulatory T cells (Tregs), and cytotoxic CD8⁺ T cells has been found to be closely linked to THCA progression (16). TAMs often exhibit an M2-like phenotype that promotes tumor growth and immune evasion through the secretion of immunosuppressive cytokines (17). Tregs contribute to local immunosuppression by inhibiting effector T cell function, while decreased infiltration or dysfunction of CD8⁺ T cells correlates with impaired anti-tumor immunity and poor clinical outcomes (18,19). Therefore, elucidating the roles of TME-related biomarkers may uncover novel therapeutic targets and inform immunomodulatory treatment strategies.

The cAMP response element modulator (CREM) is a transcription factor regulated by the cAMP signaling pathway and belongs to the CREB/ATF (cAMP response element-binding protein/activating transcription factor) family (20-22). CREM modulates gene expression by binding to cAMP response elements (CREs) in the promoters of target genes and is involved in a broad range of physiological and pathological processes (23,24). CREM also plays a key role in immune regulation, particularly in T cell biology, where it influences cytokine production and T cell activation. Given that THCA is frequently accompanied by immune infiltration and chronic inflammation, including the recruitment of Tregs and dysfunctional effector T cells, CREM may play a yet uncharacterized role in modulating immune responses within the thyroid TME. Liu et al. demonstrated that during the immunosuppressive effects of proton pump inhibitors, upregulation of CREM in T cells leads to the inhibition of pro-inflammatory cytokines such as IL-2 and IFN-γ, thereby dampening T cell activation (25). In tumorigenesis, CREM has emerged as an important regulator of the TME, a complex milieu composed of tumor cells, immune cells, and stromal components. In hepatitis B/C virus-associated hepatocellular carcinoma (HCC), elevated CREM expression correlates with increased infiltration of M2-polarized macrophages, which facilitate tumor progression by secreting immunosuppressive cytokines such as CCL18 (26). CREM has been implicated in immune modulation in liver cancer by promoting M2 macrophage infiltration; a similar mechanism may underlie its influence on the immune landscape in THCA. These macrophages contribute to immune evasion and inflammation-driven tumor growth, highlighting the role of CREM in shaping the immunological landscape of the TME. This study aimed to investigate the role of CREM in THCA progression and its association with immune cell infiltration, with the goal of evaluating its potential as a prognostic biomarker and therapeutic target in TME modulation. We present this article in accordance with the TRIPOD reporting checklist (available at https://tcr.amegroups.com/article/view/10.21037/tcr-2025-1673/rc).

Methods

Transcriptomic profiling and computational analysis in THCA

RNA-sequencing and clinical data for THCA were obtained from The Cancer Genome Atlas (TCGA) via the Genomic Data Commons (GDC) portal (https://portal.gdc.cancer.gov). A total of 566 samples were included, comprising 507 primary tumor tissues and 59 normal thyroid tissues. Clinical parameters such as age, sex, histological subtype, TNM stage, and progression-free interval (PFI) were incorporated into downstream analyses. The study was conducted in accordance with the Declaration of Helsinki and its subsequent amendments. Raw RNA-seq data were aligned using the STAR pipeline and normalized to transcripts per million (TPM). Log₂ transformation [log₂(TPM + 1)] was applied, and batch effects were corrected using the ComBat function from the “sva” R package. Quality control involved removal of duplicated entries, samples with missing clinical information, or low total read counts. Genes with zero expression in >50% of samples were filtered out, and outlier samples were identified by principal component analysis (PCA).

Survival analysis

To assess the prognostic relevance of CREM expression in THCA, survival analyses were conducted using TCGA clinical data. Patients were stratified into high and low expression groups based on the median (50th percentile) CREM expression value. Kaplan-Meier curves for PFI were generated using the “survival” package (version 3.6) and visualized with “ggplot2”. The expression value of the target gene CREM was extracted, and clinical variables (T stage, N stage, pathological stage, pathological type, CREM expression, etc.) were integrated. The samples were divided into different subgroups according to clinical variables. The Kruskal-Wallis test was used to compare the expression differences of CREM among clinical subgroups, and the significance threshold was set to (P<0.05).

Receiver operating characteristic (ROC) curve analysis

To evaluate the diagnostic utility of CREM expression, ROC curve analysis was performed using the “pROC” package (version 1.18.0). The area under the curve (AUC) was calculated to assess discriminative accuracy between tumor and normal samples. All plots were generated with “ggplot2”.

Immune infiltration profiling

The CIBERSORT algorithm (based on linear support vector regression) was used to perform deconvolution analysis on the infiltration level of immune cells in tumor tissue. An LM22 reference matrix (containing 547 characteristic genes of 22 immune cell subtypes) (27) was used as input, along with a TPM-normalized expression matrix transformed by log₂. Parameters were set as follows: perm =100 permutation test count (for P value calculation), and QN = FALSE for data normalization (since the TPM data was already normalized). Samples with P<0.05 were retained, and the relative proportions of the 22 immune cell types in each sample were obtained. Immune infiltration results were visualized using “ggplot2”.

Least absolute shrinkage and selection operator (LASSO)-Cox regression modeling for prognostic gene signature development

This study utilized LASSO-Cox regression analysis to identify the most valuable and succinct gene combinations from a vast array of genes linked to CREM expression, thereby mitigating the risk of overfitting associated with high-dimensional data in constructing a multi-gene predictive feature. Initially, utilizing transcriptome data from the complete TCGA-THCA cohort, we computed the Pearson correlation coefficients between all genes and CREM, selecting 46 genes with |R| >0.3 and P<0.05 as the preliminary candidate gene set for LASSO analysis.

The LASSO-Cox regression analysis was conducted utilizing the “glmnet” package (version 4.1.7) in R (version 4.2.0). This model does feature selection by implementing an L1 penalty on the regression coefficients, thereby reducing the coefficients of insignificant variables to zero. The analysis utilized patient survival time and survival status as dependent variables, whereas the expression levels of the previously indicated candidate genes served as independent variables. We conducted 10-fold cross-validation to ascertain the appropriate regularization value λ. In determining the value of λ, we employed the λ.1se criterion, which identifies the maximum λ value within one standard deviation of the minimal partial likelihood deviation. This criterion seeks to build a more straightforward and resilient model, even with a reduced number of variables.

Definition and validation of prognostic indicators

The final model produces a predictive feature including N genes and their associated LASSO coefficients. The risk score for each patient is determined by the subsequent formula: Risk Score = (Gene 1 expression level × Coefficient₁) + (Gene 2 expression level × Coefficient₂) + ... + (Gene N expression level × Coefficient_N). Due to the absence of an independent external validation cohort in this investigation, we employed the following internal validation indicators to assess the efficacy of this predictive feature:

• Discrimination: the C-index was employed to assess the alignment between the model’s predictions and actual observations.
• Stratification capability: patients were categorized into high-risk and low-risk categories according to the median risk score. Survival curves were generated employing the Kaplan-Meier method, and the survival disparities between the two groups were assessed utilizing the Log-rank test.

Protein-protein interaction (PPI) network and co-expression analysis

PPI networks for CREM, IL37, PSG6, and RCAN1 were constructed using the STRING database (https://string-db.org/), with the organism set to Homo sapiens. Interaction evidence was limited to experimentally validated sources, with a minimum interaction confidence score of 0.150 (low confidence) and a maximum of 50 first-shell interactors per query gene.

To further investigate transcriptional co-expression patterns, the gene correlation module of GEPIA2 was used to identify the top 200 genes most strongly associated with CREM, IL37, PSG6, and RCAN1 across TCGA datasets. Pairwise correlations were assessed using Pearson correlation coefficients, and statistical significance was determined via P values. Overlapping gene sets among the four genes were identified and visualized using Venn diagrams to explore potential shared regulatory networks.

Statistical analyses

All statistical analyses were conducted in R (version 4.2.1). Non-parametric Kruskal-Wallis tests were used to compare CREM expression levels between THCA and normal thyroid tissues. Gene expression cutoffs were defined based on median values. Associations between gene expression levels were assessed using Spearman’s rank correlation coefficients. Prognostic relevance was evaluated through univariate and multivariate Cox proportional hazards regression analyses. Covariates included age, sex, histological subtype, pathological T stage, overall pathological stage, and CREM expression group (high vs. low). Hazard ratios (HRs) and 95% confidence intervals (CIs) were reported. A two-tailed P<0.05 was considered statistically significant.

Results

CREM expression in THCA

To explore the potential role of CREM in thyroid carcinogenesis, transcriptomic data from 507 THCA samples and 59 normal thyroid tissues were retrieved from the TCGA database. A pan-cancer analysis revealed variable CREM expression across multiple tumor types (Figure 1A,1B), suggesting a context-dependent role in oncogenesis. Specifically in THCA, comparative analysis demonstrated a significant downregulation of CREM mRNA in tumor tissues relative to adjacent normal controls (P<0.001, Figure 1C). A paired comparison between tumor and matched normal tissues further confirmed this trend, with significantly reduced CREM expression in the tumor group (P<0.001, Figure 1D).

Figure 1 CREM expression and prognostic potential in THCA. (A) Box plot showing CREM expression in normal versus tumor tissues from the TCGA database. (B) Radar chart illustrating differential CREM expression across various cancers. (C) Comparison of CREM expression between thyroid cancer patients and normal individuals. (D) The difference in CREM expression between thyroid cancer patients and paired adjacent normal samples. (E) AUC analysis of the prognostic value of CREM in THCA. (F) ROC curve analysis of the diagnostic ability of CREM to distinguish thyroid cancer from normal tissues. *, P<0.05; **, P<0.01; ***, P<0.001. AUC, area under the curve; CI, confidence interval; CREM, cAMP responsive element modulator; FPR, false positive rate; ROC, receiver operating characteristic; TCGA, The Cancer Genome Atlas; THCA, thyroid carcinoma; TPM, transcripts per million; TPR, true positive rate.

To assess the predictive value of CREM expression in THCA, we performed time-dependent ROC analysis using the “timeROC” package. The resulting AUC values over a 1–5-year follow-up period indicated favorable prognostic performance, implicating CREM as a potential protective factor in THCA progression (Figure 1E). In parallel, diagnostic ROC analysis using the “pROC” package yielded a characteristic curve for CREM expression, supporting its clinical applicability as a diagnostic biomarker (Figure 1F). Collectively, these findings indicate that reduced CREM expression is significantly associated with THCA and may hold both diagnostic and prognostic relevance in clinical settings.

CREM expression is correlated with clinicopathological features and prognosis in THCA

To elucidate the clinical significance of CREM in THCA, we systematically analyzed its expression in relation to key clinicopathological parameters using data derived from the TCGA database. Our results revealed statistically significant associations between CREM expression and multiple clinical features, including pathological stage (P<0.05), histological subtype (P<0.001), and PFI events (P<0.05) (Figure 2A). These findings suggest that reduced CREM expression constitutes a clinically meaningful molecular alteration with potential implications for tumor behavior and disease progression.

Figure 2 CREM expression and clinical pathological features in THCA. (A) CREM expression levels in tumor tissues across various clinical characteristics from the TCGA database: T stage, N stage, pathological stage, histological type, and PFI events. (B) The correlation between CREM expression and PFI in THCA patients. (C) Prognostic subgroup analysis of CREM expression in THCA patients with different clinical characteristics. (D) The prognosis of different THCA subtypes varied under different CREM expression levels, with CREM expression having the most significant impact on follicular subtypes (P=0.006). ns, P>0.05; *, P<0.05; ***, P<0.001. CI, confidence interval; CREM, cAMP responsive element modulator; HR, hazard ratio; PFI, progression-free interval; TCGA, The Cancer Genome Atlas; THCA, thyroid carcinoma; TPM, transcripts per million.

Because THCA has a relatively good prognosis, progression-free survival (PFS) events clearly encompass recurrence, distant metastasis, and death from any cause. It is a comprehensive indicator that reflects disease control status early on; therefore, we used PFS as the survival analysis indicator in this study. PFS is a widely used composite endpoint in oncology that measures the time from treatment initiation (or randomization) to the first occurrence of any of the following events: (I) disease progression, defined as a ≥20% increase in the sum of target lesion diameters or the appearance of new measurable lesions; (II) recurrence: this encompasses local recurrence (the reemergence of malignancies at the source location or regional lymph nodes) and distant metastasis (the dissemination of tumors to remote organs such as the lungs and bones); (III) mortality from Any Etiology: PFS events encompass fatalities from all causes, rather than solely those attributable to specific diseases.

We further assessed the prognostic value of CREM using Kaplan-Meier survival analysis, with PFI selected as the primary outcome measure due to the generally favorable overall survival rates in THCA. High CREM expression was associated with a significantly reduced risk of progression (HR =0.53, P=0.02), indicating a potential protective role (Figure 2B). This substantial risk reduction is clinically meaningful and comparable to effect sizes of established prognostic biomarkers in THCA and other cancers. Elevated CREM expression therefore could be a potentially robust favorable prognostic marker in THCA and may help stratify patients into low-risk and high-risk groups to guide personalized follow-up strategies. Subgroup analyses stratified by clinicopathological characteristics provided additional insights (Figure 2C): patients with stage I disease (HR =0.34, P=0.03) demonstrated particularly strong associations between low CREM expression and poorer clinical outcomes. Although not statistically significant, a similar trend was observed in patients with advanced tumor stages (T3/T4; HR =1.65, P=0.07). There were differences in CREM expression levels among different THCA subtypes, with both classical and follicular subtypes showing significant differences compared to tall cell variant (TCV) (Figure S1). Subgroup prognostic analysis also showed that the prognosis of different THCA subtypes varied under different CREM expression levels, with CREM expression having the most significant impact on follicular subtypes (P=0.006) (Figure 2D).

To refine prognostic assessments, we performed Cox proportional hazards regression modeling, using PFI as the dependent variable (Table 1). Univariate analysis identified several prognostic indicators: tumor stage (T1 vs. T2–T4: HR =4.684, 95% CI: 1.688–12.993, P=0.003), pathological stage (stage I vs. II–IV: HR =2.454, 95% CI: 1.394–4.318, P=0.002), age group (≤45 vs. >45 years: HR =1.638, 95% CI: 0.935–2.869, P=0.08), and CREM expression (high vs. low: HR =0.530, 95% CI: 0.301–0.933, P=0.02). In multivariate analysis, only T stage (adjusted HR =3.195, 95% CI: 1.092–9.350, P=0.03) and CREM expression (adjusted HR =0.562, 95% CI: 0.317–0.997, P=0.049) remained independently significant, whereas pathological stage lost statistical significance (HR =2.601, 95% CI: 0.750–9.024, P=0.13).

Differential expression and functional enrichment analysis of CREM-associated genes in THCA

To elucidate the transcriptional landscape associated with CREM expression in THCA, we conducted a differential gene expression analysis using the Limma package, restricting the dataset to protein-coding transcripts. Samples were stratified into high and low CREM expression groups based on the median expression value. Applying stringent cutoffs (|log₂ fold change| >1 and adjusted P value <0.05), a total of 548 differentially expressed genes (DEGs) were identified, comprising 326 upregulated and 222 downregulated genes in the high CREM expression group (Figure 3A). Subsequent co-expression analysis in the TCGA-THCA cohort identified genes most strongly correlated with CREM. The top 25 positively and negatively correlated transcripts are visualized via heatmaps (Figure 3B,3C), highlighting key co-expression patterns that may underlie CREM-mediated regulatory pathways.

Figure 3 Differential gene expression is correlated with CREM in THCA. (A) Volcano plot of differentially expressed genes. (B) Heatmap of the top 25 genes positively correlated with CREM in THCA. (C) Heatmap of the top 25 genes negatively correlated with CREM in THCA. ***, P<0.001. CREM, cAMP responsive element modulator; THCA, thyroid carcinoma; TPM, transcripts per million.

For functional annotation, we performed Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) enrichment analyses on the top 180 genes positively correlated with CREM expression (Table S1). GO biological process (BP) enrichment revealed significant involvement in RNA splicing, protein localization to organelles and the nucleus, protein folding, responses to misfolded and unfolded proteins, Golgi vesicle-mediated transport, endoplasmic reticulum (ER) stress responses, and cellular responses to DNA damage (Figure 4A). Parallel KEGG analysis indicated enrichment in pathways related to protein processing in the ER, Hippo signaling, spliceosome assembly, AMPK signaling, mRNA surveillance, and the TGF-β signaling cascade (Figure 4A).

Figure 4 Functional enrichment of CREM-associated genes. (A) BP and KEGG pathway enrichment analysis of CREM-positive co-expressed genes. (B) BP and KEGG enrichment analysis of CREM-negative co-expressed genes. (C) EMAP visualization of BP and KEGG enrichment results for CREM-positive co-expressed genes. (D) EMAP visualization of BP and KEGG enrichment results for CREM-negative co-expressed genes. BP, biological process; CREM, cAMP responsive element modulator; EMAP, enrichment map; KEGG, Kyoto Encyclopedia of Genes and Genomes.

To further delineate the biological consequences of CREM downregulation, enrichment analyses were extended to the top 100 genes negatively correlated with CREM expression (Table S2). KEGG analysis highlighted significant associations with immune-related disorders and signaling pathways, including human T-cell leukemia virus type 1 infection, cell adhesion molecules (CAMs), antigen processing and presentation, type 1 diabetes mellitus, cytokine-cytokine receptor interaction, autoimmune thyroid disease, allograft rejection, and transplant immunology. GO enrichment further highlighted the immunological impact of low CREM expression, revealing strong involvement in T cell activation, immunoglobulin-mediated responses, antigen presentation (particularly endogenous peptide antigen processing), leukocyte adhesion and intercellular signaling, cytokine production, and lymphocyte activation cascades (Figure 4B). To enhance the interpretability of these findings, enrichment results from both KEGG and GO analyses were visualized using enrichment map (EMAP)-based structural representations, providing a global overview of the functional networks perturbed by CREM expression variation (Figure 4C,4D).

Association between CREM expression and immune landscape in THCA

To elucidate the immunological implications of CREM expression in THCA, we conducted a comprehensive immune infiltration analysis using transcriptomic data from the TCGA-THCA cohort. Patients were stratified into high and low CREM expression groups based on median expression values. Pan-cancer analysis revealed that CREM expression is positively correlated with immune cell infiltration across various tumor types. However, in THCA, a contrasting trend was observed, with elevated CREM expression exhibiting significant inverse associations with multiple immune cell subsets (Figure 5A).

Figure 5 Correlation between CREM expression and immune cell infiltration in THCA. (A) Correlation between CREM expression and 24 immune cell types in normal and THCA tissues from the TCGA database. (B) Correlation between the relative abundance of immune cell types and CREM expression in THCA tissues. (C) Infiltration patterns of Tregs, T cells CD4 resting, B cells naive, and eosinophils in high vs. low CREM expression groups. (D) Correlation between CREM expression and Treg functional markers FOXP3, CTLA-4, GITR, ICOS, as well as the Treg-specific chemokines CCL22, CCL17, and CCL5. (E) The expression of CREM in relation to other immunological checkpoints. ns, P>0.05; *, P<0.05; **, P<0.01; ***, P<0.001. CREM, cAMP responsive element modulator; TCGA, The Cancer Genome Atlas; THCA, thyroid carcinoma; TPM, transcripts per million; Treg, regulatory T cells.

In THCA, higher CREM expression was most negatively correlated with the infiltration levels of Tregs (Figure 5B). Comparative subgroup analysis between high and low CREM expression cohorts confirmed statistically significant differences in the abundance of these immune populations (P<0.001) (Figure 5C). These findings suggest that CREM may exert an immunomodulatory effect on the TME in THCA. Elevated CREM expression is associated with reduced infiltration of immunosuppressive and antigen-presenting cells, indicating a potential role in immune evasion or disrupted immune surveillance. However, further experimental studies are needed to clarify the mechanisms underlying CREM-mediated immune regulation in THCA. We also examined the link between CREM and the Treg functional markers FOXP3, CTLA-4, GITR, ICOS, as well as the Treg-specific chemokines CCL22, CCL17, and CCL5, discovering a negative correlation among all variables (Figure 5D). Simultaneously, we examined the expression of CREM in relation to immunological checkpoints, revealing a negative correlation between CREM and multiple immune checkpoints (Figure 5E).

Identification of CREM-associated predictive markers and functional network analysis in THCA

To elucidate potential CREM-related biomarkers predictive of THCA susceptibility, we employed LASSO Cox regression modeling. This penalized regression approach shrinks less informative coefficients toward zero, retaining only the most predictive features based on optimal λ regularization. The analysis identified four genes with non-zero coefficients: IL37, AL670729.3, PSG6, and RCAN1 (Figure 6A,6B). For functional interpretation, due to significant expression differences of IL37, PSG6, and RCAN1 in THCA. IL37, PSG6, and RCAN1 were selected from LASSO-Cox regression based on non-zero coefficients, indicating their strong association with PFI. As shown in results, their expression is significantly correlated with the expression of CREM (Figure S2). The subsequent analyses focused on these three protein-coding genes: IL37, PSG6, and RCAN1. To explore the functional network of CREM, we performed a combined PPI and co-expression analysis. Using the STRING database, we retrieved 50 known CREM-interacting proteins. Concurrently, the GEPIA2 platform was used to identify the top 200 genes co-expressed with CREM across TCGA tumors. Intersection analysis across four relevant gene sets (Figure 6C; Table S3) yielded a set of candidate genes for constructing an integrated PPI network (Figure 6D). A chord diagram (Figure 6E) was generated to visualize key hub molecules (node degree >10; Table S4), revealing a densely interconnected functional network potentially modulated by CREM. To assess the prognostic relevance of IL37, PSG6, and RCAN1, we conducted Kaplan-Meier survival analyses using PFI as the primary outcome. High expression of IL37 (HR =2.30, P=0.006) and PSG6 (HR =2.41, P=0.003), both inversely correlated with CREM, was associated with significantly poorer clinical outcomes (Figure 6F). In contrast, elevated RCAN1 expression (HR =0.41, P=0.003) was positively correlated with CREM, thereby improving prognosis. These findings suggest potential antagonistic interactions between CREM and IL37/PSG6, and a synergistic relationship with RCAN1, implicating these genes in CREM-mediated regulatory pathways in THCA pathogenesis.

Figure 6 LASSO model and network analysis of CREM-related genes. (A) LASSO coefficient curves for variable selection. (B) Visualization of coefficient changes in the LASSO model. (C) Venn diagram showing overlap of CREM, IL37, PSG6, and RCAN1 co-expressed genes. (D) Interaction network of CREM, IL37, PSG6, RCAN1, and related genes based on STRING analysis. (E) Key molecules in the CREM-related network (node degree >10). (F) Prognostic analysis of IL37, PSG6, and RCAN1 expression in relation to PFI in THCA patients. CI, confidence interval; CREM, cAMP responsive element modulator; HR, hazard ratio; LASSO, least absolute shrinkage and selection operator; PFI, progression-free interval; THCA, thyroid carcinoma.

IL37, PSG6, and RCAN1 expression and its correlation with immune characteristics in THCA

Building upon our earlier observation of an inverse relationship between CREM expression and immune cell infiltration in THCA, we investigated the immunoregulatory profiles of three CREM-associated genes (IL37, PSG6, and RCAN1) to evaluate their potential roles in modulating the tumor immune microenvironment. Using TCGA-derived datasets, we performed comprehensive immune infiltration analyses stratified by gene expression levels. IL37 exhibited significant positive correlations with multiple immunosuppressive and antigen-presenting cell subsets, including Tregs, DCs (both conventional and activated), macrophages, and T helper cells (Th1 and Th2), while showing a negative correlation with Th17 cells (Figure 7A). PSG6 exhibited a similar immunological profile, demonstrating positive associations with DCs, macrophages, Tregs, and Th1/Th2 cells, and a negative association with Th17 cells. In contrast, RCAN1 displayed an inverse immunophenotype: its expression negatively correlated with Tregs, DCs, cytotoxic T lymphocytes, and macrophages, while positively associating with natural killer (NK) cell subsets. Comparative heatmap analysis (Figure 7B) revealed that CREM and RCAN1 share similar immune correlation profiles, distinct from the concordant patterns observed for IL37 and PSG6. Collectively, these data support a model in which CREM-associated gene networks differentially shape immune cell composition in the THCA’s TME.

Figure 7 IL37, PSG6, and RCAN1 expression is correlated with immune cell infiltration in THCA. (A) Correlation between immune cell abundance and IL37, PSG6, and RCAN1 expression in THCA tissues. The point size reflects the absolute magnitude of the Spearman correlation coefficient (|R|), indicating the strength of the association. (B) Relationship between CREM, IL37, PSG6, RCAN1 expression levels and immune cell infiltration in normal and tumor tissues. ns, P>0.05; *, P<0.05; **, P<0.01; ***, P<0.001. CREM, cAMP responsive element modulator; THCA, thyroid carcinoma.

Discussion

The CREM gene has been implicated in the pathogenesis of various malignancies (28); however, its specific functional and prognostic roles in THCA remain insufficiently elucidated. In this study, we assessed the diagnostic and prognostic relevance of CREM in THCA using transcriptomic and clinical data from the TCGA database. CREM was found to be significantly downregulated in thyroid tumor tissues compared to adjacent normal thyroid tissues, consistent with its putative role as a tumor suppressor (29). The HR is a fundamental metric in survival analysis, utilized to evaluate the risk of an endpoint event (such as PFS) between two cohorts. In the context of this investigation, an HR of 0.53 indicates that patients exhibiting elevated CREM expression possess a 0.53 times risk of disease progression, relapse, or mortality compared to those with low expression. The P value of 0.028 shows a strong negative connection between CREM expression and prognosis. This differential expression pattern supports its potential diagnostic utility, as evidenced by a robust AUC of 0.751 in ROC analysis. Furthermore, lower CREM expression was significantly linked to advanced pathological T staging and higher overall tumor stage, indicating its role in disease progression. These observations are concordant with findings from Ye et al. (30), who reported CREM downregulation in association with poor clinical features across multiple cancer types. Kaplan-Meier survival analysis further demonstrated that elevated CREM expression is associated with prolonged PFI, confirming its prognostic value. While Kaplan-Meier analysis suggests that CREM is a favorable prognostic marker, future multivariate Cox regression incorporating age, stage, and subtype is needed to confirm whether this association is independent of established clinical factors. Collectively, these results highlight CREM as a potential biomarker for both early detection and risk stratification in THCA.

The results of our immune cell-related analyses provide significant insights into the role of CREM in THCA, particularly with respect to its interactions with various immune cell populations. Notably, Tregs exhibited a negative correlation with CREM expression. Tregs are essential for maintaining immune tolerance and suppressing anti-tumor immunity, which can contribute to tumor progression (31-33). The negative correlation between CREM expression and Treg infiltration indicates that lower levels of CREM may create an immunosuppressive TME, which could weaken effective anti-tumor immune responses. This conclusion was confirmed after our analysis of T cell regulatory markers such as FOXP3, CTLA-4, GITR, and ICOS. Furthermore, we observed a negative correlation between CREM expression and Treg-specific chemokines CCL22, CCL17, and CCL5 and immunological checkpoints. We suggest that CREM, functioning as a transcriptional regulator, presumably directly suppresses the transcription of genes including CCL22, CCL17, and CCL5, based on current evidence. This may immediately result in a decrease in Treg infiltration. Furthermore, weaker recruitment signals hinder the successful recruitment of Tregs expressing specific receptors, such as CCR4, into the TME, hence diminishing the total level of Treg infiltration. Concurrently, CREM may also reduce the expression of molecules such as FOXP3 and CTLA-4 in other cells (or Tregs themselves) inside the TME, leading to diminished immunosuppressive action despite the presence of a restricted number of Tregs. This elucidates the inverse correlation between the expression of Treg functional markers and CREM. This finding also indicates a ubiquitous and significant negative connection between CREM and the immunosuppressive TME in THCA. CREM is probably a nuclear negative regulator of the immune microenvironment in THCA, and its diminished expression induces a pervasive immunosuppressive state by synergistically enhancing Treg recruitment and various immune checkpoint pathways, resulting in tumor progression and unfavorable prognosis.

Moreover, the enrichment analysis revealed that genes downregulated and co-expressed with CREM are primarily involved in immune system responses, including the regulation of T cell activation and antigen processing and presentation. These pathways are integral to the adaptive immune response, which is essential for tumor recognition and elimination. For example, somatic recombination of immune receptors, a process vital for generating a diverse repertoire of T cell receptors, plays a critical role in effective immune surveillance against tumors (34). This suggests that CREM may be a key regulator of the immune landscape in THCA, potentially influencing the efficacy of therapeutic strategies, such as immune checkpoint inhibitors.

In addition to CREM, the IL37 gene has emerged as a critical mediator of immune responses, particularly in cancer. As a member of the interleukin family, IL37 is known for its anti-inflammatory properties and its role in modulating immune responses (35-37). Recent studies emphasize the significance of regulatory B cells (Bregs) in humans. These cells primarily suppress autoimmune diseases via IL37 and play a crucial role in maintaining immune tolerance (38). In our analysis, we observed a negative correlation between IL37 expression and CREM levels, suggesting that reduced CREM expression may be associated with elevated IL37 levels. This relationship suggests that CREM might protect against THCA by enhancing immune evasion, which could contribute to tumor progression. PSG6, another gene identified in our study, is part of the pregnancy-specific glycoprotein family and has been implicated in various immune processes. PSG6 has been shown to exert immunoregulatory effects by inducing human monocytes to secrete IL-10, IL-6, and TGF-β1 (39). Our findings reveal a negative correlation between PSG6 expression and CREM, suggesting that higher levels of PSG6 may contribute to a more immunosuppressive TME. This association is particularly relevant, as PSG6’s role in promoting immune tolerance could facilitate tumor growth and metastasis.

In contrast, RCAN1 demonstrated a positive correlation with CREM expression. Previous studies have identified RCAN1, especially its isoform RCAN1.4, as a suppressor of metastatic progression (MPS) in various cancers, including THCA, which often experiences metastatic dormancy (40,41). The positive correlation between RCAN1 and CREM implies that reduced CREM levels may lower RCAN1 expression, potentially influencing tumor behavior. Since RCAN1 is linked to the suppression of THCA metastasis, its relationship with CREM indicates that higher CREM expression may help slow the progression of THCA. Understanding the interactions among these genes could provide valuable insights into the molecular mechanisms underlying THCA and highlight potential therapeutic targets for intervention. We propose that CREM, in conjunction with IL37, PSG6, and RCAN1, forms a bidirectional regulatory network defined by “inhibition-synergy”, collectively impacting THCA proliferation through two avenues: modification of the immune microenvironment and tumor cell biological behavior. The essential processes encompass immunosuppressive imbalance, signaling pathway intersection, and prognostic association. The elevated expression of CREM and RCAN1 correlates with favorable PFS, whereas diminished CREM expression alongside heightened IL37/PSG6 expression forecasts adverse prognosis, thereby affirming the substantial influence of the synergistic interactions among these four genes on tumor progression. This study offers preliminary insights; the precise process necessitates additional experimental validation.

The identified and confirmed network demonstrates that CREM, as a pivotal node, ultimately impedes THCA advancement by concurrently diminishing pro-tumor influences and augmenting tumor-suppressive mechanisms. This offers a novel perspective for comprehending the biological behavior of THCA, especially in subgroups with favorable prognoses. CREM, in conjunction with IL37, PSG6, and RCAN1, establishes a bidirectional regulatory network characterized by “inhibition-synergy”, collectively impacting THCA growth across two dimensions: immune microenvironment modification and tumor cell biological behavior. The fundamental mechanisms pertain to immunosuppressive imbalance, signaling pathway intersection, and prognostic correlation. Elevated CREM and RCAN1 expression correlates with good PFS, but diminished CREM coupled with elevated IL37/PSG6 indicates unfavorable prognosis, so affirming the significant impact of the synergistic interaction among these four genes on tumor progression. The precise mechanisms necessitate additional experimental verification.

Although our research reveals the key role of CREM in regulating the immune microenvironment in the context of THCA, numerous questions remain and require further investigation. First, this study was based exclusively on bioinformatics analyses of publicly available transcriptomic and clinical datasets from the TCGA-THCA cohort. Future studies will be necessary to confirm protein-level expression (e.g., via IHC or Western blot) and to assess the biological effects of CREM on THCA cell proliferation, migration, and immune modulation through functional assays. Second, immune heterogeneity, as inferred from single-sample Gene Set Enrichment Analysis (ssGSEA) or pan-immune gene signatures, was not explicitly incorporated as covariates in survival models. Future studies should evaluate whether CREM remains prognostically significant when adjusting for detailed immune landscape features. Third, the reliance on TCGA bulk RNA-seq data, which lacks cell-type resolution and may be confounded by tumor heterogeneity or stromal contributions. Single-cell RNA-seq and spatial transcriptomics could provide finer insights into CREM’s role in specific cell populations. Although the present study is based on comprehensive bioinformatics analyses, experimental validation is essential to establish causality. Future studies should conduct CREM knockdown and overexpression in TPC-1 and BCPAP cell lines to evaluate proliferation, migration, and secretion of immunosuppressive chemokines (CCL2/CCL22). In addition, tumor cell-macrophage/Treg co-culture experiments should be performed to confirm direct regulation of M2 polarization and Treg recruitment. ChIP-seq and luciferase reporter assays could be used to verify direct transcriptional regulation of IL37, PSG6, and RCAN1 by CREM. These experiments will further validate the mechanistic role of CREM and support its therapeutic target.

Conclusions

In conclusion, our study identifies CREM as a key biomarker for predicting THCA prognosis and highlights its potential as a therapeutic target in immunotherapy, particularly due to its association with immune cell infiltration. The interactions among CREM, IL37, PSG6, and RCAN1 show that more research is needed to understand their roles in THCA. Investigating these molecular interactions could lead to important insights for developing targeted therapeutic strategies that improve patient outcomes.

Acknowledgments

The authors would like to express their gratitude to EditSprings (https://www.editsprings.cn) for the expert linguistic services provided.

Footnote

Reporting Checklist: The authors have completed the TRIPOD reporting checklist. Available at https://tcr.amegroups.com/article/view/10.21037/tcr-2025-1673/rc

Peer Review File: Available at https://tcr.amegroups.com/article/view/10.21037/tcr-2025-1673/prf

Funding: None.

Conflicts of Interest: All authors have completed the ICMJE uniform disclosure form (available at https://tcr.amegroups.com/article/view/10.21037/tcr-2025-1673/coif). The authors have no conflicts of interest to declare.

Ethical Statement: The authors are accountable for all aspects of the work in ensuring that questions related to the accuracy or integrity of any part of the work are appropriately investigated and resolved. The study was conducted in accordance with the Declaration of Helsinki and its subsequent amendments.

Open Access Statement: This is an Open Access article distributed in accordance with the Creative Commons Attribution-NonCommercial-NoDerivs 4.0 International License (CC BY-NC-ND 4.0), which permits the non-commercial replication and distribution of the article with the strict proviso that no changes or edits are made and the original work is properly cited (including links to both the formal publication through the relevant DOI and the license). See: https://creativecommons.org/licenses/by-nc-nd/4.0/.

References

1. Seib CD, Sosa JA. Evolving Understanding of the Epidemiology of Thyroid Cancer. Endocrinol Metab Clin North Am 2019;48:23-35. [Crossref] [PubMed]
2. Ravella L, Lopez J, Descotes F, et al. DICER1 mutated, solid/trabecular thyroid papillary carcinoma in an 11-year-old child. Ann Pathol 2018;38:316-320. [Crossref] [PubMed]
3. Ma LJ, Wang SX, Zhu GQ. Multiple papillary carcinoma of thyroid with extensive metastasis of the lateral neck and parapharyngeal: a case report. Lin Chuang Er Bi Yan Hou Tou Jing Wai Ke Za Zhi 2019;33:675-6. [Crossref] [PubMed]
4. Megwalu UC, Moon PK. Thyroid Cancer Incidence and Mortality Trends in the United States: 2000-2018. Thyroid 2022;32:560-70. [Crossref] [PubMed]
5. Pizzato M, Li M, Vignat J, et al. The epidemiological landscape of thyroid cancer worldwide: GLOBOCAN estimates for incidence and mortality rates in 2020. Lancet Diabetes Endocrinol 2022;10:264-72. [Crossref] [PubMed]
6. Hogan AR, Zhuge Y, Perez EA, et al. Pediatric thyroid carcinoma: incidence and outcomes in 1753 patients. J Surg Res 2009;156:167-72. [Crossref] [PubMed]
7. Lebbink CA, van den Broek MFM, Kwast ABG, et al. Opposite Incidence Trends for Differentiated and Medullary Thyroid Cancer in Young Dutch Patients over a 30-Year Time Span. Cancers (Basel) 2021;13:5104. [Crossref] [PubMed]
8. Kebebew E, Clark OH. Differentiated thyroid cancer: "complete" rational approach. World J Surg 2000;24:942-51. [Crossref] [PubMed]
9. Kluijfhout WP, Pasternak JD, Drake FT, et al. Application of the new American Thyroid Association guidelines leads to a substantial rate of completion total thyroidectomy to enable adjuvant radioactive iodine. Surgery 2017;161:127-33. [Crossref] [PubMed]
10. Yin X, Lu C, Sun D, et al. Stimulating thyroglobulin to TSH ratio predict long-term efficacy of (131)I therapy in patients with differentiated thyroid cancer after total thyroidectomy: a retrospective study. Endocrine 2024;84:1064-71. [Crossref] [PubMed]
11. Jiang YJ, Han ZJ, Hu YX, et al. Family history of malignant or benign thyroid tumors: implications for surgical procedure management and disease-free survival. Front Endocrinol (Lausanne) 2023;14:1282088. [Crossref] [PubMed]
12. Huang G, Chen J, Zhou J, et al. Epigenetic modification and BRAF gene mutation in thyroid carcinoma. Cancer Cell Int 2021;21:687. [Crossref] [PubMed]
13. Catalano MG, Fortunati N, Boccuzzi G. Epigenetics modifications and therapeutic prospects in human thyroid cancer. Front Endocrinol (Lausanne) 2012;3:40. [Crossref] [PubMed]
14. Liu Q, Sun W, Zhang H. Roles and new Insights of Macrophages in the Tumor Microenvironment of Thyroid Cancer. Front Pharmacol 2022;13:875384. [Crossref] [PubMed]
15. Qin L, Chen C, Gui Z, et al. IRX5's influence on macrophage polarization and outcome in papillary thyroid cancer. Front Oncol 2024;14:1399484. [Crossref] [PubMed]
16. Fang D, Zhou L, Zheng B. Research Progress on the Immunological Correlation Between Papillary Thyroid Carcinoma and Hashimoto's Thyroiditis. J Immunol Res 2025;2025:7192808. [Crossref] [PubMed]
17. Huang R, Kang T, Chen S. The role of tumor-associated macrophages in tumor immune evasion. J Cancer Res Clin Oncol 2024;150:238. [Crossref] [PubMed]
18. Ueyama A, Nogami W, Nashiki K, et al. Immunotherapy Targeting CCR8+ Regulatory T Cells Induces Antitumor Effects via Dramatic Changes to the Intratumor CD8+ T Cell Profile. J Immunol 2023;211:673-82. [Crossref] [PubMed]
19. Mao FY, Kong H, Zhao YL, et al. Increased tumor-infiltrating CD45RA(-)CCR7(-) regulatory T-cell subset with immunosuppressive properties foster gastric cancer progress. Cell Death Dis 2017;8:e3002. [Crossref] [PubMed]
20. Molina CA, Foulkes NS, Lalli E, et al. Inducibility and negative autoregulation of CREM: an alternative promoter directs the expression of ICER, an early response repressor. Cell 1993;75:875-86. [Crossref] [PubMed]
21. Mioduszewska B, Jaworski J, Kaczmarek L. Inducible cAMP early repressor (ICER) in the nervous system--a transcriptional regulator of neuronal plasticity and programmed cell death. J Neurochem 2003;87:1313-20. [Crossref] [PubMed]
22. Seidl MD, Nunes F, Fels B, et al. A novel intronic promoter of the Crem gene induces small ICER (smICER) isoforms. FASEB J 2014;28:143-52. [Crossref] [PubMed]
23. Chang JH, Vuppalanchi D, van Niekerk E, et al. PC12 cells regulate inducible cyclic AMP (cAMP) element repressor expression to differentially control cAMP response element-dependent transcription in response to nerve growth factor and cAMP. J Neurochem 2006;99:1517-30. [Crossref] [PubMed]
24. Servillo G, Della Fazia MA, Sassone-Corsi P. Coupling cAMP signaling to transcription in the liver: pivotal role of CREB and CREM. Exp Cell Res 2002;275:143-54. [Crossref] [PubMed]
25. Liu W, Jakobs J, Rink L. Proton-Pump Inhibitors Suppress T Cell Response by Shifting Intracellular Zinc Distribution. Int J Mol Sci 2023;24:1191. [Crossref] [PubMed]
26. Song G, Shi Y, Zhang M, et al. Global immune characterization of HBV/HCV-related hepatocellular carcinoma identifies macrophage and T-cell subsets associated with disease progression. Cell Discov 2020;6:90. [Crossref] [PubMed]
27. Bindea G, Mlecnik B, Tosolini M, et al. Spatiotemporal dynamics of intratumoral immune cells reveal the immune landscape in human cancer. Immunity 2013;39:782-95. [Crossref] [PubMed]
28. Yu K, Kuang L, Fu T, et al. CREM Is Correlated With Immune-Suppressive Microenvironment and Predicts Poor Prognosis in Gastric Adenocarcinoma. Front Cell Dev Biol 2021;9:697748. [Crossref] [PubMed]
29. Fenton MS, Marion KM, Hershman JM. Identification of cyclic adenosine 3',5'-monophosphate response element modulator as an activator of the human sodium/iodide symporter upstream enhancer. Endocrinology 2008;149:2592-606. [Crossref] [PubMed]
30. Ye M, Wu H, Mei Y, et al. High expression of CREM is associated with poor prognosis in gastric cancer patients. Nan Fang Yi Ke Da Xue Xue Bao 2024;44:1776-82. [Crossref] [PubMed]
31. Muth S, Klaric A, Radsak M, et al. CD27 expression on Treg cells limits immune responses against tumors. J Mol Med (Berl) 2022;100:439-49. [Crossref] [PubMed]
32. Kim JH, Kim BS, Lee SK, Regulatory T. Cells in Tumor Microenvironment and Approach for Anticancer Immunotherapy. Immune Netw 2020;20:e4. [Crossref] [PubMed]
33. Ke X, Wang J, Li L, et al. Roles of CD4+CD25(high) FOXP3+ Tregs in lymphomas and tumors are complex. Front Biosci 2008;13:3986-4001. [Crossref] [PubMed]
34. Kim JW, Kim S, Yang SY, et al. T-Cell Receptor Repertoire Characteristics Associated with Prognostic Significance in High-Grade Serous Ovarian Carcinoma. Genes (Basel) 2023;14:785. [Crossref] [PubMed]
35. Sarhan D, Hippen KL, Lemire A, et al. Adaptive NK Cells Resist Regulatory T-cell Suppression Driven by IL37. Cancer Immunol Res 2018;6:766-75. [Crossref] [PubMed]
36. Nold MF, Nold-Petry CA, Zepp JA, et al. IL-37 is a fundamental inhibitor of innate immunity. Nat Immunol 2010;11:1014-22. [Crossref] [PubMed]
37. Zhao M, Li Y, Guo C, et al. IL-37 isoform D downregulates pro-inflammatory cytokines expression in a Smad3-dependent manner. Cell Death Dis 2018;9:582. [Crossref] [PubMed]
38. Wang L, Christodoulou MI, Jin Z, et al. Human regulatory B cells suppress autoimmune disease primarily via interleukin-37. J Autoimmun 2025;153:103415. [Crossref] [PubMed]
39. Snyder SK, Wessner DH, Wessells JL, et al. Pregnancy-specific glycoproteins function as immunomodulators by inducing secretion of IL-10, IL-6 and TGF-beta1 by human monocytes. Am J Reprod Immunol 2001;45:205-16. [Crossref] [PubMed]
40. Nairon KG, Nigam A, Khanal T, et al. RCAN1.4 regulates tumor cell engraftment and invasion in a thyroid cancer to lung metastasis-on-a-chip microphysiological system. Biofabrication 2024;
41. Wang C, Saji M, Justiniano SE, et al. RCAN1-4 is a thyroid cancer growth and metastasis suppressor. JCI Insight 2017;2:e90651. [Crossref] [PubMed]