Multi-omics analysis identifies RARRES1 as a potential biomarker linked to immunosuppressive microenvironment and its radiomics prediction in ovarian cancer

Introduction

Ovarian cancer (OV) ranks as the eighth most common malignancy and the eighth leading cause of cancer-related mortality among women worldwide. In 2022, there were approximately 324,398 newly diagnosed cases and 206,839 deaths globally, with both incidence and mortality continuing to increase over the past twelve years (1). In addition to genetic predisposition, environmental and occupational exposures—such as endocrine-disrupting chemicals—have been identified as important risk factors for OV in women aged 20–49 years, suggesting that the global burden of this disease is likely to rise further (2). Despite these trends, effective early screening and diagnostic approaches for OV remain limited. Widely used biomarkers such as cancer antigen 125 (CA125) are challenged by variable sensitivity across populations, underscoring the urgent need for improved strategies (3).

Retinoic acid receptor responder 1 (RARRES1), originally identified as a retinoic acid receptor-induced membrane protein, is upregulated in skin and other tissues in response to retinoic acid and its agonists (4). In many cancers, RARRES1 expression is frequently silenced by promoter hypermethylation, and it has been implicated in the regulation of cellular metabolism, lipid biosynthesis, and glucose metabolic reprogramming. In hepatocellular carcinoma, RARRES1 inhibits cell proliferation and migration while enhancing sensitivity to lenvatinib through interaction with SPINK2 (5). In triple-negative breast cancer, integrative analyses have revealed a strong association between RARRES1 expression, tumor immune microenvironment (TIME) signatures, and responsiveness to immune checkpoint inhibitors, suggesting that RARRES1 may serve as a prognostic and predictive biomarker for immunotherapy (6). Although several recent studies in OV have suggested the potential clinical relevance of RARRES1 from the perspectives of prognostic modeling, metastasis-related genes, and the immune microenvironment, its cell type–specific association with tumor immunosuppressive states, the systematic integrative evidence linking it to the mutational landscape and co-expression networks, as well as radiomics-based strategies for non-invasive prediction of its expression status, remain to be further refined and more consistently elucidated.

In fact, the prognostic value of RARRES1 in OV and its associated biological context have been reported in multiple studies. As early as 2019, Wang et al. constructed a five-gene prognostic signature (IGF2, PEG3, DCN, LYPD1, and RARRES1) based on The Cancer Genome Atlas (TCGA) and Gene Expression Omnibus (GEO) datasets, and demonstrated that the high-risk group was enriched in WNT/β-catenin signaling and EMT-related pathways, thereby incorporating RARRES1 into a transcriptomic model for overall survival prediction (7). Furthermore, Zhang et al., in a study focusing on omental metastasis in OV, developed a six-gene OMAG prognostic model (TIMP3, FBN1, IGKC, RPL21, UCHL1, and RARRES1). Through integrated single-cell and bulk analyses, they observed expression patterns of RARRES1 in epithelial/CSC-like populations and fibroblasts, suggesting its potential involvement in metastasis-associated microenvironmental remodeling (6). In addition, stemness- and stem cell pathway-based classification studies in HGSOC have explored molecular subtypes and risk tools linked to immune scores, tumor purity, and immune escape, further supporting the significance of the “stemness-immune microenvironment” axis in disease progression (8). Beyond prognostic signatures, RARRES1 has also been implicated in shaping immune microenvironment-related phenotypes and has been subjected to functional validation. For example, studies on tumor microenvironment (TME) regulation and risk stratification have suggested associations between RARRES1, immune-related subgroups, and prognosis. In vitro experiments further showed that upregulation of RARRES1 could inhibit OV cell proliferation, migration, and invasion, indicating that it may participate in key regulatory networks influencing tumor behavior (2). On the other hand, Mendelian randomization analyses in genetic epidemiology have also evaluated the causal relationship between RARRES1 and OV risk, providing complementary clues for mechanistic investigation (9). However, most of these studies have primarily focused on exploratory associations with prognosis or immune relevance, or have been restricted to specific data types and single-layer validation. A systematic and unified framework remains lacking to elucidate the cellular origins of the RARRES1-associated immunosuppressive state, particularly the intralineage heterogeneity within macrophage populations, together with cross-omics integration involving the mutational landscape and co-expression modules. Moreover, the establishment of a radiogenomic closed loop that maps molecular states onto computed tomography (CT) imaging phenotypes to enable non-invasive prediction has yet to be comprehensively achieved.

Recent studies have increasingly implicated RARRES1 in OV prognosis, immune contexture, stemness-related programs, and metastatic behavior, suggesting that its association with disease progression and immune modulation is supported by accumulating evidence. Accordingly, the present study does not aim to redefine RARRES1 as a novel prognostic marker, but rather to consolidate and extend these findings within an integrated multi-omics framework. Specifically, we combined bulk RNA-seq, single-cell transcriptomics, somatic mutation profiling, and co-expression network analysis to refine the cellular and molecular context associated with RARRES1 expression, with particular attention to macrophage heterogeneity and immunosuppressive programs. Most importantly, we developed and validated a CT-based radiomics machine-learning model to non-invasively predict RARRES1 expression status. By linking imaging phenotypes with molecular and immune features, this radiomics-driven approach represents the key original contribution of the present work and provides a potential pathway for clinical translation beyond transcriptomic association studies. We present this article in accordance with the TRIPOD reporting checklist (available at https://tcr.amegroups.com/article/view/10.21037/tcr-2025-1-2593/rc).

Methods

Data acquisition and preprocessing

Transcriptomic (RNA-seq) and somatic mutation data for 428 OV samples were downloaded from TCGA-OV project. Single-cell RNA-seq data from seven OV tumor samples (GSE184880) were obtained from the GEO database. For radiomics analysis, CT images and corresponding region-of-interest (ROI) segmentations were acquired from The Cancer Imaging Archive (TCIA, https://www.cancerimagingarchive.net/). Only pre-treatment CT scans with complete DICOM series and manual gross tumor volume (GTV) annotations were included. ROIs were reviewed by two board-certified radiologists to ensure annotation accuracy, and discrepancies were resolved by consensus. DICOM data were converted to NIfTI format for preprocessing and feature extraction.

Prognostic and differential expression analysis

Kaplan-Meier (KM) survival analysis was performed to evaluate the prognostic significance of RARRES1 in OV. Patients were divided into high- and low-expression groups based on the median RARRES1 expression value. Differentially expressed genes (DEGs) were identified using the DESeq2 package with thresholds set at |log₂FC| ≥1 and false discovery rate (FDR) ≤0.05. The top 30 upregulated and downregulated genes were visualized in a z-score normalized heatmap.

scRNA-seq data analysis

Basic processing and cell type annotation

scRNA-seq data were processed using the Seurat package. Raw unique molecular identifier (UMI) counts were normalized to 10,000 reads per cell using the LogNormalize method. The top 2,000 highly variable genes were identified using the vst method, and all genes were scaled and centered with ScaleData. Batch effects were corrected using the Harmony algorithm. Principal component analysis (PCA) was performed, and the top 10 principal components (PCs) were selected based on elbow plot inspection for downstream analysis. A k-nearest neighbor (KNN) graph was constructed using Harmony-corrected embeddings, and clustering was performed across a resolution gradient (0.01–3.0), with resolution=1 selected for final clustering. Two-dimensional embeddings were generated using UMAP and t-SNE. Cell types were annotated in two steps: (I) automatic annotation using the SingleR package with the Human Primary Cell Atlas as a reference, and (II) manual refinement based on canonical marker genes and published marker gene lists, resulting in the identification of nine major cell types, including T cells, NK cells, macrophages, fibroblasts, monocytes, epithelial cells, B cells, tissue stem cells, and endothelial cells. No hepatocyte population was retained after manual refinement. Macrophages were further subclustered into distinct subpopulations.

Trait-associated cell scoring with scPagwas

The scPagwas framework was applied to integrate GWAS and scRNA-seq data in order to identify trait-relevant cell populations. GWAS data were obtained from the UK Biobank VCF file (ukb-b-18798.vcf.gz) and converted into a summary statistics format containing chrom, pos, rsid, beta, se, and maf. Single-cell transcriptomic data were provided as a Seurat object. Cell type contributions and single-cell trait relevance scores (TRS) were calculated based on KEGG pathway annotations, hg37 genomic block annotation, and linkage disequilibrium (LD) reference data. Bootstrap resampling was performed to evaluate the stability of the results. The derived TRS values were subsequently used for downstream visualization and subpopulation-level analyses.

Cell-cell communication analysis

Cell-cell communication networks were inferred using CellChat with the Secreted Signaling database. Overexpressed signaling genes and ligand-receptor pairs were identified with identifyOverExpressedGenes and identifyOverExpressedInteractions, and protein-protein interaction (PPI) networks were integrated to enhance reliability. Communication probabilities were calculated using computeCommunProb (min.cells =3), aggregated to the pathway level, and visualized using circle plots, bubble plots, chord diagrams, and centrality analysis.

Bayesian deconvolution

Macrophage subpopulation proportions in TCGA-OV bulk RNA-seq samples were estimated using BayesPrism. Preprocessed single-cell count matrices and cell type annotations were used as reference data, retaining only protein-coding genes. Outlier genes were removed based on bulk-single-cell consistency plots, and final proportions were computed using the run.prism and get.fraction functions.

GSEA and transcription factor (TF) activity analysis

DEGs among macrophage subpopulations were ranked by average log₂ fold change and subjected to gene set enrichment analysis (GSEA) using KEGG pathway gene sets and the clusterProfiler package. The top three significantly enriched pathways were visualized with gseaplot2, and all significant pathways were displayed using ridge plots. TF regulatory networks were analyzed using the DoRothEA package with high-confidence human regulons (A/B/C levels). TF activity scores were calculated with the VIPER algorithm, stored in the “dorothea” assay of the Seurat object, and visualized through PCA, neighbor graph construction, and clustering. Average TF activity per cell type was calculated after z-score normalization.

Somatic mutation and WGCNA analysis

Somatic mutation profiles were analyzed using maftools. Samples were divided into RARRES1-high and RARRES1-low groups. Mutation landscape summaries, including the top 30 mutated genes, mutation types, and transition/transversion ratios, were generated using plotmafSummary. Differential mutation frequencies were assessed with mafCompare, and odds ratios (ORs) with 95% confidence intervals were visualized. For WGCNA, genes were first ranked by mean expression across the TCGA-OV cohort. In parallel, differential expression analysis was performed using DESeq2 by comparing RARRES1-high versus RARRES1-low tumors. Genes showing significant differential expression (|log2FC| ≥1, FDR ≤0.05) were preferentially retained, and the final gene set for network construction consisted of the top 10,000 genes derived from the combined consideration of mean expression and DESeq2 results. Soft-thresholding power was determined using pickSoftThreshold, and adjacency matrices were converted to topological overlap matrices (TOMs). Initial modules were merged at a correlation threshold of 0.25, resulting in 13 modules. Module-trait correlations were calculated using Pearson’s correlation, with hub genes defined as |MM| >0.8 and within the top 5% of gene significance scores.

Radiomics feature extraction and machine learning model construction

Feature preprocessing and selection were performed in multiple steps. First, features with very low variance were removed to reduce redundancy. Next, an F-test was applied to rank features by relevance to the outcome, and non-linear expansions were tested on top-ranked features. Data were then split into training and testing sets with a 6:4 ratio under stratified sampling. For model construction, four machine learning classifiers were compared: random forest (RF), gradient boosting (GB), XGBoost (XGB), and KNN. Hyperparameters for each model were optimized via fivefold cross-validated randomized search, using the area under the receiver operating characteristic curve (ROC-AUC) as the primary metric. Model performance was evaluated on both training and independent testing sets using accuracy, precision, recall, F1-score, and ROC-AUC. Calibration curves, precision-recall curves, and confusion matrices were generated to assess discrimination and calibration. Feature importance was quantified for tree-based models (RF, GB, XGB), and SHapley Additive exPlanations (SHAP) values were further computed to visualize the contribution of individual radiomics features to the final predictions.

Quantitative reverse transcription polymerase chain reaction (qRT-PCR)

Fresh primary ovarian carcinoma tissues were collected from 12 patients in the Department of Oncology, Xianning Central Hospital, between June 2024 and April 2025. Total RNA was extracted from human OV tissues using a total RNA extraction kit (Beyotime, Shanghai, China). mRNA was reverse transcribed into cDNA with the PrimeScript One Step RT-PCR Kit (Takara, Kyoto, Japan). qRT-PCR was performed on a LightCycler 96 system (Roche, Basel, Switzerland) using the TB Green Premix Ex Taq II Kit (Takara). Glyceraldehyde 3-phosphate dehydrogenase (GAPDH) was used as an internal reference, and the relative expression levels of CD8A, CD68, CD163 and RARRES1 were calculated using the 2^−ΔΔCT method. The names and sequences of the primers used in the experiment are shown in Table 1. The study was conducted in accordance with the Declaration of Helsinki and its subsequent amendments. The study was approved by Xianning Central Hospital ethics committee (No. 2025HU24), and all participants provided informed consent.

Western blot

Frozen human OV tissues were lysed in RIPA buffer containing phenylmethylsulfonyl fluoride (PMSF) (Beyotime, Shanghai, China). After centrifugation at 4 ℃ (12,000 ×g, 10 min), the supernatants were mixed with 5× loading buffer (1/5 of the supernatant volume) and boiled in a 100 ℃ water bath for 5–10 min. Equal amounts of protein were separated by 10% SDS-PAGE and transferred to polyvinylidene fluoride (PVDF) membranes. The membranes were incubated with specific primary antibodies against CD8A (rabbit mAb, CST, #85336, 1:1,000), CD68 (rabbit mAb, ABclonal, A20803, 1:1,000), CD163 (rabbit mAb, CST, #68922, 1:1,000), RARRES1 (rabbit pAb, Proteintech, 27491-1-AP, 1:1,000) and β-actin (rabbit mAb, CST, #4970, 1:2,000) (all from the indicated manufacturers), followed by incubation with HRP-conjugated goat anti-rabbit IgG (H&L) secondary antibody (Thermo Fisher, A16110, 1:5,000). The bands were visualized on a chemiluminescence imaging system using an enhanced chemiluminescence exposure solution.

Statistical analysis

Survival differences were assessed using the log-rank test. TME metrics and immune-related gene expression differences were visualized using boxplots. All expression values were log₂-transformed. Between-group comparisons were conducted using two-sided Wilcoxon rank-sum tests, with significance thresholds set at P<0.05 (*P<0.05, **P<0.01, ***P<0.001; ns, not significant).

Results

Identification of OV trait-associated macrophage subpopulations using scPagwas

Figure 1 presents the technical roadmap of this paper. Quality control (QC) was performed on scRNA-seq data derived from seven OV tissue samples. The distributions of nFeature_RNA, nCount_RNA, mitochondrial gene percentage (percent.mt), and ribosomal gene percentage (percent.rb) indicated overall good data quality without significant outliers (Figure 2A). After batch effect correction and integration using the Harmony algorithm, cells exhibited uniform distribution in the two-dimensional UMAP space (Figure 2B). Based on canonical marker genes, nine major cell types were annotated, including T cells, NK cells, macrophages, fibroblasts, monocytes, epithelial cells, B cells, tissue stem cells, and endothelial cells (Figure 2C and Figure S1A) (10-13). Notably, the proportions of these cell types varied substantially across OV samples (Figure 2D), suggesting considerable inter-patient heterogeneity within the TME.

Figure 1 The technical roadmap of this paper. CT, computed tomography; FFPE, formalin fixation and paraffin-embedding; GSEA, gene set enrichment analysis; HGSOC, high-grade serous ovarian cancer; IL-17, interleukin 17; NF-κb, nuclear factor-kappaB; OV, ovarian cancer; PM, proteostasis-associated macrophage; qPCR, quantitative polymerase chain reaction; RARRES1, retinoic acid receptor responder 1; ROC, receiver operating characteristic; SAM, stress-activated macrophage; SHAP, SHapley Additive exPlanations; TAM, tumor-associated macrophage; TCGA, The Cancer Genome Atlas; TCIA, The Cancer Imaging Archive; TLR, toll like receptor; TNF, tumor necrosis factor; WB, western blot; WGCNA, weighted gene co-expression network analysis.

Figure 2 Basic analysis of scRNA-seq data. (A) Distributions of nFeature_RNA, nCount_RNA, percent.mt, and percent.rb after quality control of OV samples. (B) Two-dimensional UMAP plot after Harmony-based batch effect correction and integration. (C) Cell type annotation based on canonical markers. (D) Proportions of cell types across OV samples. (E) Distribution of macrophage subtypes (left), UMAP visualization of scPagwas.TRS.Score (middle), and violin plots of TRS scores across macrophage subtypes (right). (F) Proportional distribution of macrophage subtypes across OV samples. (G) PCC heatmap of macrophage subpopulation transcriptomes with hierarchical clustering. (H) Bootstrap resampling analysis of macrophage subpopulation scores. OV, ovarian cancer; PCC, Pearson correlation coefficient; TRS, trait relevance scores; UMAP, uniform manifold approximation and projection.

Among all immune cell subsets, macrophages displayed remarkable population diversity. Further subdivision revealed four major macrophage subpopulations: pro-inflammatory M1 macrophages (M1), M2-like tumor-associated macrophages (TAMs), stress-activated macrophages (SAMs), and proteostasis-associated macrophages (PMs) (Figure 2E and Figure S1B, left) (14-18). UMAP visualization based on the scPagwas.TRS.Score demonstrated distinct functional activities among these subpopulations (Figure 2E, middle), while violin plots further quantified the TRS score distributions across the four subsets (Figure 2E, right).

Analysis of the relative proportions of macrophage subpopulations across different OV samples revealed notable variability (Figure 2F). Pearson correlation analysis of transcriptional profiles showed specific correlation patterns between functional subtypes (Figure 2G), with hierarchical clustering clearly separating them into distinct clusters. To further explore the dynamic relationships among macrophage subpopulations, pseudotime trajectory analysis was performed, revealing that M1, TAM, SAM, and PM cells were distributed along a shared trajectory with gradual divergence, consistent with a continuum of activation states rather than fixed lineages (Figure S2). Finally, bootstrap resampling analysis indicated that M1 macrophages exhibited consistently higher scores compared to other subpopulations, underscoring their functional significance within the OV microenvironment (Figure 2H).

TAM-centered inflammatory programs and transcriptional regulators associated with RARRES1

Pathway and regulatory analyses anchored the radiogenomic signal in a TAM-centered, immuno-inflammatory program. GSEA on TAM-associated genes highlighted significant enrichment of inflammation and innate/adaptive immune pathways—NF-κB, IL-17, Toll-like receptor, and TNF signaling—together with tumor-relevant processes such as lipid metabolism and cytokine–receptor interactions (Figure 3A), consistent with the canonical role of NF-κB/TLR–TNF axes in immune activation and inflammatory transcriptional control. Transcription-factor activity profiling inferred with DoRothEA (coupled to VIPER-style footprinting) showed higher activity of NFKB1, RELA, STAT3, and CEBPA within TAMs, coherently placing TAMs at the intersection of inflammatory and myeloid-lineage regulatory circuits (Figure 3B).

Figure 3 Mechanistic support linking the predefined RARRES1 state to immune context. (A) GSEA enrichment results of TAMs. (B) Heatmap of TF activities across subtypes. (C) Comparison of stromal, immune, ESTIMATE, and tumor purity scores between low- and high-RARRES1 expression groups, along with representative top DEGs. The x-axis represents sample groups; the y-axis represents normalized expression values or scores (log₂-transformed). Four TME indicators are shown: Stromal⁺ (stromal infiltration score), Immune⁻(immune infiltration score), Estimate⁻ (ESTIMATE score; lower scores indicate lower tumor purity), and Purity⁺ (tumor purity score). (D) Module-trait correlation heatmap, showing associations between modules and macrophage subtypes (red = positive, green = negative; P values in parentheses). ***, P<0.001. DEGs, differentially expressed genes; GSEA, gene set enrichment analysis; RARRES1, retinoic acid receptor responder 1; TAMs, tumor-associated macrophages; TF, transcription factor; TME, tumor microenvironment.

Bulk-level immune context in TCGA-OV further supported this pattern: compared with the low-expression group, high RARRES1 tumors displayed lower stromal, lower immune, and lower ESTIMATE scores but higher tumor purity, indicating a reduction of non-malignant components in the tissue microenvironment (Figure 3C). Finally, WGCNA module-trait correlation heatmaps linked specific co-expression modules to macrophage subtypes, delineating systems-level covariation between the TAM/M1 programs and trait readouts (Figure 3D). Collectively, Figure 3 ties TAM-centered signaling pathways and TF activity to bulk immune composition and network modules, providing compact mechanistic support that bridges the predefined RARRES1 state with the immune context of the OV microenvironment.

Functional enrichment analysis of macrophage-associated gene modules further clarified their biological roles (Figure S3A,S3B). Genes within the magenta module, which showed strong positive correlation with M1 macrophages, were significantly enriched in immune activation-related pathways, including antigen processing and presentation, interferon signaling, leukocyte-mediated immunity, and ribosome biogenesis, consistent with a transcriptionally active, pro-inflammatory and immune-stimulatory phenotype. In contrast, genes in the pink module, preferentially associated with TAMs, were enriched in pathways related to extracellular matrix organization, cell-matrix adhesion, lysosome and phagosome function, lipid metabolism, and cytokine-cytokine receptor interaction, reflecting a macrophage program oriented toward tissue remodeling, metabolic adaptation, and immune regulation rather than cytotoxic immune activation.

Tissue‑level validation of RARRES1 and immune‑marker expression

To provide exploratory tissue-level support for the transcriptomic and single-cell findings, we analyzed 12 fresh primary ovarian carcinoma specimens using qRT-PCR and Western blotting (Figure 4A). Based on median RARRES1 mRNA expression, samples were stratified into a RARRES1-low group (n=6) and a RARRES1-high group (n=6). As expected, RARRES1 expression was higher in the RARRES1-high group at both mRNA and protein levels (Figure 4B-4D). Notably, tumors in the RARRES1-high group showed lower CD8A expression, together with increased CD68 and CD163 levels, suggesting reduced cytotoxic T-cell infiltration and enrichment of macrophage populations with M2-like features. Given the limited sample size, these results should be interpreted as preliminary and descriptive, but they are directionally consistent with the immune-cold, macrophage-dominant TME inferred from the TCGA-OV bulk analyses and scRNA-seq based macrophage profiling.

Figure 4 Tissue expression of RARRES1 and immune‑markers in RARRES1‑stratified ovarian carcinoma cohorts. (A) Sample grouping: 12 primary ovarian carcinoma specimens stratified into RARRES1‑low (n=6) vs. RARRES1-high (n=6) groups. (B) qRT-PCR results for RARRES1, CD8A, CD68 and CD163. (C) Representative Western blot for RARRES1, CD8A, CD68 and CD163. (D) Quantification of protein expression. Data are presented as mean ± standard deviation; ***, P<0.001. qRT-PCR, quantitative reverse transcription polymerase chain reaction; RARRES1, retinoic acid receptor responder 1.

CT‑derived radiomics signature of RARRES1 high vs. low in ovarian carcinoma

Using pre-treatment CT images from the TCGA-OV radiogenomics cohort, we identified a radiomics signature that robustly discriminated tumors with high versus low RARRES1 expression (Figure 5A). Among the extracted features, first-order intensity-based descriptors consistently emerged as the most predictive components of the final random-forest model (Figure 5B,5C). Specifically, gray-level distribution features reflecting central tendency and lower-intensity tails, including the median intensity, mean intensity, and 10th-percentile gray-level values, showed the highest global importance scores (Figure 5D). At the population level, RARRES1-high tumors were characterized by systematically altered intratumoral attenuation profiles, suggesting differences in tissue density and compositional heterogeneity captured on routine CT imaging. These intensity-based features were more informative than higher-order texture metrics, indicating that global voxel intensity characteristics rather than fine-grained spatial textures primarily underlie the radiogenomic association with RARRES1 expression. Model interpretability analysis using SHAP further demonstrated that these first-order features contributed coherently at the individual-patient level. In representative cases, elevated median and mean gray-level values increased the probability of classification as RARRES1-high, whereas lower-intensity distributions favored RARRES1-low predictions (Figure 5E). This consistency between global feature importance and local explanations supports the biological relevance of the identified radiomics signature. Importantly, RARRES1-high status was independently associated with poor overall survival and an immune-cold TME in the transcriptomic analyses. Therefore, the radiomics features predictive of RARRES1 expression may also carry indirect prognostic relevance, serving as non-invasive surrogates of an immunosuppressive and macrophage-dominant tumor phenotype. Collectively, these results establish that simple, interpretable CT-derived intensity features are sufficient to capture clinically meaningful molecular and immune states, linking radiologic phenotype to RARRES1-associated tumor biology.

Figure 5 Locked radiomics model discriminates predefined RARRES1 expression. (A) Test-set CT exemplars with tumor ROI overlays; inset shows the time-based 6:4 split → model lock → held-out evaluation workflow. (B) Test ROC (random forest). (C) Test confusion matrix. (D) SHAP global feature importance. (E) SHAP case-level waterfall explanations. Radiomics features were extracted with PyRadiomics using IBSI-aligned settings; data from TCIA TCGA-OV Radiogenomics (HGSOC). AUC, area under the curve; CT, computed tomography; HGSOC, high-grade serous ovarian cancer; IBSI, Image Biomarker Standardisation Initiative; OV, ovarian cancer; ROC, receiver operating characteristic; ROI, region-of-interest; SHAP, SHapley Additive exPlanations; TCGA, The Cancer Genome Atlas; TCIA, The Cancer Imaging Archive.

Discussion

In this study, we systematically investigated the expression characteristics, immunological roles, and clinical significance of RARRES1 in OV from a multi-omics and multi-scale perspective. Analyses of the TCGA-OV cohort demonstrated that high RARRES1 expression was closely associated with poor prognosis and accompanied by decreased stromal and immune components as well as increased tumor purity within the TME. These findings suggest that RARRES1 may accelerate tumor progression by reshaping the TME, reducing immune cell infiltration, and promoting immune evasion. Differential immune gene expression analysis further revealed that elevated RARRES1 expression was associated with downregulation of antigen processing and presentation, chemokine, interferon, and natural killer (NK) cell cytotoxicity pathways.

Notably, existing literature suggests that RARRES1 may exert context-dependent roles across tumor types and biological settings. Early studies characterized RARRES1 as a tumor suppressor frequently silenced by promoter hypermethylation, with inhibitory effects on proliferation, migration, and metabolic reprogramming in hepatocellular carcinoma and lymphoma models (19,20). Similarly, functional experiments in OV cell lines reported that RARRES1 overexpression suppressed tumor cell growth and invasiveness, suggesting an intrinsic anti-tumor role (21). In contrast, several recent transcriptomic and systems-level studies—including prognostic risk models, stemness-related signatures, and single-cell analyses—have consistently associated higher RARRES1 expression with poor survival, immune suppression, or aggressive molecular subtypes in OV and triple-negative breast cancer (6,7,22). Our results align more closely with this latter body of evidence, supporting the view that RARRES1 expression in patient tumors may reflect an unfavorable immune microenvironment rather than acting solely as a classical tumor suppressor.

Single-cell RNA-seq analyses provided additional cellular insights into RARRES1-associated immune suppression. Macrophages exhibited marked heterogeneity and were subdivided into four functional subpopulations, including M1 macrophages, TAMs, SAMs, and PMs. Notably, TAMs occupied a central position in the intercellular communication network and showed strong enrichment of inflammatory signaling pathways such as NF-κB, IL-17, Toll-like receptor, and TNF signaling. At first glance, this enrichment of inflammatory pathways appears counterintuitive given the immunosuppressive role of TAMs. However, gene set enrichment analysis reflects the steady-state expression of downstream effector genes, whereas TF activity inferred by DoRothEA/VIPER captures upstream regulatory potential rather than final pathway output. In TAMs, elevated activity of inflammatory and stress-responsive TFs such as NFKB1, RELA, STAT3, and CEBPA likely reflects chronic inflammatory stimulation within the TME. Importantly, in cancer contexts these transcriptional programs are frequently repurposed to drive immune suppression, metabolic adaptation, and tissue remodeling rather than effective anti-tumor immunity. Importantly, this apparent pro-inflammatory profile does not contradict the classical M1/M2 paradigm but instead highlights the critical distinction between inflammatory signaling and effective anti-tumor immunity within the TME. The inflammatory programs observed in TAMs are more consistent with a form of chronic, tumor-associated inflammation, which is known to promote immune suppression, tissue remodeling, and tumor progression rather than cytotoxic immune activation. By contrast, M1 macrophages displayed relatively stronger features related to antigen presentation and immune activation, consistent with their canonical anti-tumor roles. Importantly, this TAM-centered inflammatory program was closely associated with high RARRES1 expression, which correlated with reduced CD8⁺ T-cell infiltration, increased CD68⁺/CD163⁺ macrophage abundance, lower immune and stromal scores, and higher tumor purity. Collectively, these findings support the concept that RARRES1 marks an inflammatory but immune-suppressive TAM phenotype, providing a mechanistic explanation for why enhanced inflammatory signaling coexists with an immune-cold microenvironment and poor clinical outcome in RARRES1-high OV.

Several of the identified immunoregulatory factors have previously been implicated in OV biology. TGF-β promotes OV development through induction of TGFBI, a macrophage-secreted protein that fosters an immunosuppressive TME and represents a potential therapeutic target (23). MIF, a known biomarker candidate, enhances early OV diagnosis when combined with CA125, CA19-9, or EGFR detection (24). SPP1 is closely linked to T cell exhaustion and poor prognosis in OV, and its knockdown induces tumor cell apoptosis, highlighting its therapeutic potential (14). NF-Κb-related pathways mediate POSTN-driven OV progression by enhancing tumor invasion and recruitment of TAMs and CAFs (25). IL-17 and Th17 cells modulate OV immune responses, with IL-17A levels in ascites correlating with improved survival. TLR2 expression is elevated in late-stage OV and may serve as a prognostic marker (26). STAT3 is frequently overexpressed, promoting tumor growth, angiogenesis, stemness, and therapeutic resistance (27). CEBPB facilitates OV cell proliferation and invasion via SOS1-ERK1/2 activation, and may contribute to PARP inhibitor resistance (28). Collectively, these factors converge on TAM-mediated immune suppression, supporting our findings.

Mutation profiling further highlighted that OV was dominated by missense mutations and SNPs, with TP53, TTN, and CSMD3 among the most frequently mutated genes. Importantly, RARRES1 expression groups exhibited distinct mutational landscapes, with certain mutations enriched in the low-expression group. WGCNA identified modules strongly correlated with M1 macrophages, enriched in ribosome biogenesis and protein synthesis pathways, potentially reflecting the elevated translational demand during pro-inflammatory immune responses. At the tissue level, our qRT-PCR and Western-blot validation strengthens this multi-omics story: tumors in the RARRES1-high group showed significantly lower CD8 expression and markedly higher CD68 and CD163 levels compared with the RARRES1-low group, consistent with an immune-cold, macrophage-rich microenvironment.

Radiomics analysis demonstrated that an RF-based machine learning model could accurately predict RARRES1 expression levels, achieving high discriminative performance and predictive consistency across both training and testing cohorts. These findings suggest that radiomics may serve as a non-invasive surrogate biomarker for RARRES1 status, facilitating patient stratification and treatment decision-making. This molecular-imaging integration strategy offers a novel approach for precision oncology.

Several limitations of this study should be acknowledged. First, the tissue-level validation was conducted on a small number of clinical specimens (n=6 vs. 6) and should therefore be regarded as exploratory rather than definitive. These experiments were designed to provide proof-of-concept support for the multi-omics findings, and their results require confirmation in larger, independent cohorts with sufficient statistical power. Second, although we observed a robust association between high RARRES1 expression and macrophage-dominant, immune-cold features across bulk transcriptomics, single-cell analyses, and preliminary tissue assays, direct causal mechanism—such as whether RARRES1 actively drives macrophage polarization—were not addressed experimentally. Future studies incorporating functional assays will be essential to elucidate the mechanistic role of RARRES1 in TME remodeling.

Conclusions

This study comprehensively elucidated the multifaceted roles of RARRES1 in OV. Our findings suggest that a subgroup of HGSOC with high RARRES1 expression may have a more aggressive immune-cold phenotype, and thus RARRES1 holds potential as a stratification biomarker for progression and immunotherapy-guidance. Furthermore, the radiomics-based prediction model demonstrated that RARRES1 expression can be reliably estimated through non-invasive imaging features, supporting its application as a predictive biomarker. Collectively, these results provide new insights into the mechanistic involvement of RARRES1 in OV and offer valuable evidence for its clinical translation in personalized therapeutic strategies.

Acknowledgments

None.

Footnote

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

Data Sharing Statement: Available at https://tcr.amegroups.com/article/view/10.21037/tcr-2025-1-2593/dss

Peer Review File: Available at https://tcr.amegroups.com/article/view/10.21037/tcr-2025-1-2593/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-1-2593/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. The study was approved by Xianning Central Hospital ethics committee (No. 2025HU24), and all participants provided informed consent.

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. Zhang R, Zhao J, Du P, et al. Global burden and trends in ovarian cancer attributable to environmental risks and occupational risks in females aged 20-49 from 1990 to 2021, with projections to 2050: a cross-sectional study. BMC Public Health 2025;25:2547. [Crossref] [PubMed]
2. Zhao A, Zhou S, Yang X, et al. Transcription factor networks and novel immune biomarkers reveal key prognostic and therapeutic insights in ovarian cancer. Discov Oncol 2025;16:309. [Crossref] [PubMed]
3. Momenimovahed Z, Mazidimoradi A, Allahqoli L, et al. The Role of CA-125 in the Management of Ovarian Cancer: A Systematic Review. Cancer Rep (Hoboken) 2025;8:e70142. [Crossref] [PubMed]
4. Wang Z, He K, Liu M, et al. Enhanced mitochondrial biogenesis facilitates the development of cutaneous squamous cell carcinoma. Cancer Lett 2025;618:217623. [Crossref] [PubMed]
5. Guo Y, Chai B, Zhang H, et al. RARRES1 inhibits hepatocellular carcinoma progression and increases its sensitivity to lenvatinib through interaction with SPINK2. Biol Direct 2024;19:15. [Crossref] [PubMed]
6. Yu Z, Liu H, Ye J, et al. Integrative analysis identifies cancer cell-intrinsic RARRES1 as a predictor of prognosis and immune response in triple-negative breast cancer. Front Genet 2024;15:1360507. [Crossref] [PubMed]
7. Wang R, Ye XH, Zhao XL, et al. Development of a five-gene signature as a novel prognostic marker in ovarian cancer. Neoplasma 2019;66:343-9. [Crossref] [PubMed]
8. Wu H, Li D, Sun L, et al. Determination of high-grade serous ovarian cancer stem cell-based subtypes and prognostic model and identification of highly expressed VSIG4 and STAB1 in macrophages. J Ovarian Res 2025;18:159. [Crossref] [PubMed]
9. Liu S, Lin H, Zhang K, et al. Potential drug targets for ovarian cancer identified through Mendelian randomization and colocalization analysis. J Ovarian Res 2025;18:32. [Crossref] [PubMed]
10. Tirosh I, Izar B, Prakadan SM, et al. Dissecting the multicellular ecosystem of metastatic melanoma by single-cell RNA-seq. Science 2016;352:189-96. [Crossref] [PubMed]
11. Aran D, Looney AP, Liu L, et al. Reference-based analysis of lung single-cell sequencing reveals a transitional profibrotic macrophage. Nat Immunol 2019;20:163-72. [Crossref] [PubMed]
12. Zhang X, Lan Y, Xu J, et al. CellMarker: a manually curated resource of cell markers in human and mouse. Nucleic Acids Res 2019;47:D721-8. [Crossref] [PubMed]
13. Zhang D, Lu W, Cui S, et al. Establishment of an ovarian cancer omentum metastasis-related prognostic model by integrated analysis of scRNA-seq and bulk RNA-seq. J Ovarian Res 2022;15:123. [Crossref] [PubMed]
14. Wang J, Zhu N, Su X, et al. Novel tumor-associated macrophage populations and subpopulations by single cell RNA sequencing. Front Immunol 2023;14:1264774. [Crossref] [PubMed]
15. Revel M, Sautès-Fridman C, Fridman WH, et al. C1q+ macrophages: passengers or drivers of cancer progression. Trends Cancer 2022;8:517-26. [Crossref] [PubMed]
16. Xu L, Chen Y, Liu L, et al. Tumor-associated macrophage subtypes on cancer immunity along with prognostic analysis and SPP1-mediated interactions between tumor cells and macrophages. PLoS Genet 2024;20:e1011235. [Crossref] [PubMed]
17. Bahrami S, Drabløs F. Gene regulation in the immediate-early response process. Adv Biol Regul 2016;62:37-49. [Crossref] [PubMed]
18. Zhang W, Yang J, Wang B, et al. HRS mediates tumor immune evasion by regulating proteostasis-associated interferon pathway activation. Cell Rep 2023;42:113352. [Crossref] [PubMed]
19. Maimouni S, Lee MH, Sung YM, et al. Tumor suppressor RARRES1 links tubulin deglutamylation to mitochondrial metabolism and cell survival. Oncotarget 2019;10:1606-24. [Crossref] [PubMed]
20. Patel J, Xun D, Creswell K, et al. Loss of RARRES1 function Promotes Follicular Lymphomagenesis and Inhibits B cell Differentiation in Mice. Int J Biol Sci 2022;18:2670-82. [Crossref] [PubMed]
21. Li H, Wang Y, Li G, et al. Integrative analysis of cuproptosis-associated genes for predicting immunotherapy response in single-cell and multi-cohort studies. J Gene Med 2024;26:e3600. [Crossref] [PubMed]
22. Li N, Yu K, Huang D, et al. Identifying a Novel Eight-NK Cell-related Gene Signature for Ovarian Cancer Prognosis Prediction. Curr Med Chem 2024;31:1578-94. [Crossref] [PubMed]
23. Lecker LSM, Berlato C, Maniati E, et al. TGFBI Production by Macrophages Contributes to an Immunosuppressive Microenvironment in Ovarian Cancer. Cancer Res 2021;81:5706-19. [Crossref] [PubMed]
24. Ghose A, McCann L, Makker S, et al. Diagnostic biomarkers in ovarian cancer: advances beyond CA125 and HE4. Ther Adv Med Oncol 2024;16:17588359241233225. [Crossref] [PubMed]
25. Lin SC, Liao YC, Chen PM, et al. Periostin promotes ovarian cancer metastasis by enhancing M2 macrophages and cancer-associated fibroblasts via integrin-mediated NF-κB and TGF-β2 signaling. J Biomed Sci 2022;29:109. [Crossref] [PubMed]
26. Sobstyl M, Niedźwiedzka-Rystwej P, Hrynkiewicz R, et al. Toll-like Receptor 2 as a Marker Molecule of Advanced Ovarian Cancer. Biomolecules 2021;11:1205. [Crossref] [PubMed]
27. Guo S, Ramar V, Guo AA, et al. TRPM7 transactivates the FOSL1 gene through STAT3 and enhances glioma stemness. Cell Mol Life Sci 2023;80:270. [Crossref] [PubMed]
28. Tan J, Wang D, Tu A, et al. CEBPB regulates ERK1/2 activity through SOS1 and contributes to ovarian cancer progression. Med Oncol 2025;42:242. [Crossref] [PubMed]