The gene expression, mutation, copy number variation, and clinical information of TCGA dataset were retrieved from UCSC Xena (19) and cBioPortal (20). The expression matrices of GSE14520, GSE36376, and GSE77314 were downloaded from the Gene Expression Omnibus (GEO) database, along with their corresponding survival information. The gene expression data of TCGA and GSE77314 were in the form of log2-transformed RSEM counts for further analyses. Background correction and normalization: For the data from GSE14520 and GSE36376, background correction was carried out using the robust multi-array average (RMA) for arrays from both Affymetrix and Illumina. Subsequently, normalization was performed using quantile normalization. Probe annotation: After microarray data normalization, probes were annotated using a manufacturer-provided reference file for GPL571 (Affymetrix Human Genome U133A 2.0 Array) and GPL10558 (Illumina HumanHT-12 V4.0 expression beadchip). For probes matching more than one gene, the probes were omitted. If several probes matched one gene, only the probe with the highest mean signaling intensity across samples in the dataset was retained. This process was implemented using an in-house R script for data manipulation and comparison of probe intensities.

The Gene Ontology gene list “CELL ACTIVATION INVOLVED IN IMMUNE RESPONSE” term was downloaded from the MSIGDB database (21). Cox univariate regression: In TCGA dataset, Cox univariate regression analysis was performed using the coxph function from the R package “survival”. The formula for the univariate regression was Surv(time, event) ~ gene_expression, where time represents the survival time, event indicates the event status (e.g., death), and gene_expression is the expression level of each gene in the gene list. The significance level for this analysis was set at α = 0.05. Grouped survival difference analyses: Survival data in TCGA dataset were grouped based on the median expression level of each gene. The high/low-expression group was defined as the expression level higher/lower than the median of the samples. The Kaplan–Meier survival curves were plotted for the high- and low-expression groups, and the log-rank test was used to assess the significance of survival differences between the two groups. In this step, the R package “survival” was used. Genes that were significantly associated with overall survival according to both Cox univariate regression (p < 0.05) and grouped survival difference analyses (p < 0.05) were retained as candidate genes. Panel optimization: Panel optimization was carried out by enumerating all combinations of candidate genes (gene number <7) using the combn function in R. In this step, all possible combinations of the candidate genes selected were listed and used for further analysis. For each combination, a Cox multivariate regression model was constructed in TCGA dataset using the coxph function with the formula Surv(time, event) ~ combination_of_genes. The sample datasets were equally divided into high- and low-risk groups based on the risk scores calculated from the Cox multivariate regression model. Survival differences between the high- and low-risk groups were evaluated using the log-rank test as described above, and p-values for all combinations were recorded. The combination with the smallest p-value from the log-rank test was retained as the optimized panel. Survival analysis was performed using the R package “survival”, and gene expression was visualized using “pheatmap”.

The signature scores of the GSE14520, GSE36376, and GSE77314 datasets were calculated using the same formula used for the training dataset. The samples were also equally separated into low- and high-risk samples according to the corresponding median values. For clinical information, risk score differences were calculated using Student’s t-test and visualized using the R package “vioplot”. The results of Cox multivariate regression were visualized using “survminer”. Based on the relevant data in TCGA and the model, the nomogram and calibration curve were generated and visualized using “rms”.

Non-synonymous mutations in TCGA dataset were used to construct a mutation matrix. The differences in mutations between the high- and low-risk samples were evaluated using Fisher’s exact test, and the mutational landscape was visualized using the R package “GenVisR” (22). Copy number differences were evaluated using the Wilcoxon rank sum test and visualized using “pheatmap”. The differentially expressed genes between groups were identified using the R package “limma” with an adjusted p-value <0.01 and a |log2-fold change| > 2. The significantly enriched Gene Ontology (molecular function, cellular component, and biological process) and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway terms were evaluated and visualized using the R package “clusterProfiler” (23). For gene set enrichment analyses, the genes were ranked by fold change and visualized using default parameters.

Infiltration was estimated using different algorithms based on the transcriptome of TCGA dataset, including XCELL (24), CIBERSORT (25), TIMER (26), QUANTISEQ (27), EPIC (28), and TIMER2.0. Differences in infiltration between groups were quantified using Student’s t-test and visualized using the R packages “vioplot” and “pheatmap”. Similarly, the IC50 values of each TCGA sample were quantified using gene expression values and a default training dataset with the R package “OncoPredict” (29). The differences in the IC50 values were also determined and visualized.

The single-cell sequencing R data for 14 processed HCC patients were retrieved from the GEO database (accession number GSE156625). The data were processed according to a previously reported method (30), with the annotation results. The data were processed using the R package “Seurat” (31). Pseudo RNA-seq data were generated using Seurat::AggregateExpression, the expression data were log2 transformed to calculate the risk score values of these 14 samples, and the cell proportions were subsequently visualized. The R packages “CellChat” (32) and “Patchwork” were used for intercellular communication analysis with default parameters.

The cells were treated with trypsin, resuspended in a complete medium, and counted. For each group, 2,000 cells per well in 100 µL of medium were seeded in triplicate. The cells were allowed to fully settle, the cell density was observed, and the cells were cultured. Starting from the second day after plating, 10 µL of CCK-8 solution was added to each well. After 1–3 hours, the plate was gently shaken for 2–5 minutes, and a spectrophotometer was used to measure the optical density (OD) at 450 nm.

Chambers were placed into an empty 24-well plate, 100 µL of serum-free medium was added to the chambers, and the chambers were hydrated for 1–2 hours. The cells were digested with trypsin, suspended, and counted. The medium was removed, and 600 µL of medium containing 30% fetal bovine serum (FBS) was added to the lower chamber. The cell suspension was diluted with serum-free medium, and 100 µL of the suspension was added to each chamber. The chambers were transferred into a medium containing 30% FBS and incubated for 4–24 hours. The medium was removed; the cells were fixed and incubated at room temperature for 10–30 minutes. The upper chamber was washed with 1× Phosphate Buffered Saline (PBS), the cells were stained by immersing them in a staining solution for 5–10 minutes, the excess stain was absorbed, and non-migratory cells were removed, washed again, air-dried, and photographed.

Chromatin immunoprecipitation (ChIP) was conducted utilizing a ChIP assay kit (Millipore Sigma, Burlington, MA, USA) following the manufacturer’s instructions. The eluted DNA was immunoprecipitated with either an IgG or an anti-RORC antibody (14-6988-82, eBioscience, San Diego, CA, USA). Quantification of each immunoprecipitated DNA sample was performed via quantitative PCR (qPCR) using primers specifically designed to amplify the proximal promoter region of the CDC6 gene, encompassing predicted RORC binding sites. All assays were conducted in triplicate, and the results were normalized to the input DNA.

The CDC6 3′UTR was cloned into the psi-CHECK2 vector and transfected into the cells using Lipofectamine 3000 under standard conditions at 37°C for 4 hours. Luciferase activity was assessed after 48 hours of incubation using the dual-luciferase reporter assay system (Promega, Madison, WI, USA) with measurements taken at 490 nm. Firefly luciferase activity was normalized by calculating the ratio of firefly to Renilla luciferase activity. The primer sequences were listed as follows: RORC-F: GTAACGCGGCCTACTCCTG, RORC-R: GTCTTGACCACTGGTTCCTGT; GAPDH-F: ACAACTTTGGTATCGTGGAAGG, GAPDH-R: GCCATCACGCCACAGTTTC; CDC6-F: AAGCCCTG, CDC6-R: TCAAATACCAATCTTCGTCCC. Sequences of siCDC6 were 5′-AGACUAUAACUCUACAGAUUGUGdAdA-3′. The siRORC sequences were siRORC-1 (UUAAAGUGCACAUAAGAGGUUCCUCUUAUGUGCACUUUAAAG), siRORC-2 (UUUAAAGUGCACAUAAGAGGUCUCUUAUGUGCACUUUAAAGA), siRORC-3 (UCUUUAAAGUGCACAUAAGAGCUUAUGUGCACUUUAAAGAUA).

Primary anti-CDC6 antibody (1:500; 11640‐1‐AP; Proteintech, Chicago, IL, USA) and anti-RORC antibody (1:500, 14-6988-82, eBioscience) were utilized for immunohistochemical (IHC) staining. At least three pathologists not informed of the clinical, pathological, or clinical outcome of all patients participated in quantifying the signal intensity of the samples. The samples were deparaffinized, hydrated, blocked, and added into the primary anti-CDC6 rabbit antibody (diluted 1:500) and cultured for 10 h at 4°C.

The correlation between the expression level of the immune cell activation gene and overall survival in TCGA dataset was estimated using Cox univariate regression, and survival differences between the high-expression and low-expression groups were determined. In total, 22 genes were identified. To facilitate clinical practice and remove redundant information, combinations of genes with fewer than six genes were enumerated. Signatures were constructed by applying Cox multivariate regression on different gene combinations. The combination with the smallest p-value was identified as the optimized panel and signature. As a result, a signature based on six genes (SLC11A1, RORC, NKG7, ITM2A, DNASE1L3, and CLCF1) was constructed. The signature values were calculated as follows: signature = (0.07577395 * SLC11A1) + (−0.15675603 * RORC) + (−0.18440305 * NKG7) + (−0.01589646 * ITM2A) + (−0.08330426 * DNASE1L3) + (0.07315706 * CLCF1). The samples were divided into the high- and low-risk groups, and the median value across TCGA samples was used as the cutoff. The high-risk samples were characterized by significantly shorter survival (median survival: 23.8 months, 95% CI: 20.6–41.8) than the low-risk samples (median: 83.2, 95% CI: 69.6–NA, p < 0.001) ( Figure 1A ). SLC11A1 and CLCF1 were overexpressed in the high-risk samples, while RORC, NKG7, ITM2A, and DNASE1L3 were more highly expressed in the low-risk samples ( Figure 1B ). Consistent with this finding, the high-risk samples also had a significantly worse progression-free survival rate than the low-risk samples ( Figure 1C ). The 3-year survival ROC curve was visualized ( Figure 1D ), and the area under the curve (AUC) was calculated. Compared to clinical indicators, including age, sex, and stage, the signature performed better. Collectively, these results suggest that the signature is a valuable biomarker for predicting the prognosis of HCC patients.

To verify the robustness of the model, it was validated across three entirely independent GEO datasets: GSE14520, GSE36376, and GSE77314. The signature score was calculated according to the expression values and the corresponding coefficients in the above formula. The samples in each dataset were also divided into the high- and low-risk groups based on the median value in each dataset. The high-risk samples had a significantly shorter survival period than the low-risk samples in the GSE14520 ( Figure 2A ), GSE36376 ( Figure 2B ), and GSE77314 ( Figure 2C ) datasets. Consistently, the gene expression patterns of the candidate genes also resembled those in the training dataset, TCGA ( Figures 2A–C , bottom panel). In summary, the performance of the signature was robust and reproducible across datasets instead of being a result of overfitting.

The relationships between the signature and clinical indicators were evaluated. The signature was significantly associated with pathological stage but was independent of age and sex ( Figure 3A ). Cox multivariate regression revealed that the signature was significantly associated with survival (p < 0.001), while age and sex were not ( Figure 3B ). The immune status and immunotherapy biomarker PD-L1 were significantly negatively associated with the signature ( Figure 3C ), while PD-1 was not. A nomogram was constructed to predict the 3-year survival of patients with HCC in TCGA dataset ( Figure 3D ). The risk score had a broader range of points (from 0 to 100), indicating a greater impact on the total score compared to age, sex, and stage, which have narrower ranges. The calibration curve was also visualized ( Figure 3E ). These results indicate that the signature reflects some clinical indicators but contains more information for predicting the prognosis of HCC patients.

As mentioned, multiple genes/omic signatures may reflect various biological statuses of cancer, and the relationships between the signature and mutation, copy number variation (CNV), and transcriptome were investigated next. The mutation rate of TP53 was significantly different between the low- and high-risk groups, and the mutation rate was significantly higher in the high-risk group ( Figure 4A , right panel). At the CNV level, variations were more frequent ( Figure 4B ) in the high-risk samples. Differentially expressed genes between the high- and low-risk groups were identified, and the enriched Gene Ontology (including molecular function, cellular component, and biological processes) was also identified ( Figures 4C–E ). Cancer progression-related processes and pathways were also identified. Gene Set Enrichment Analysis (GSEA) revealed that canonical cancer pathways were enriched in the high-risk group ( Figure 4F ). In addition to the cell cycle, cell adhesion, and drug metabolism pathways ( Figure 4G ), immune response-related pathways, including cytokine–cytokine interaction, chemokine signaling, and NK cell-mediated cytotoxicity, were also significantly enriched ( Figure 4H ). Overall, the signature reflected the biological status of HCC at both the genome and transcriptome levels.

The signature was significantly associated with various immune pathways and PD-L1 expression, and the signature was constructed with immune cell activation-related genes, which prompted us to analyze immune cell infiltration and the signature. Immune cell infiltration was estimated with multiple algorithms (including CIBERSORT, TIMER, XCELL, and EPIC) using the transcriptome of TCGA dataset. A high proportion of immune cells was differentially infiltrated according to these algorithms ( Figure 5A ). CD8+ T cells, M0 macrophages, and NK cells were differentially infiltrated between the low- and high-risk groups ( Figure 5B ). Compared to the candidate genes, the signature more significantly and comprehensively reflected immune infiltration ( Figure 5C ). Compared to the candidate genes, the signature was significantly correlated with all immune cell types listed, while the candidate genes were not. Collectively, the results above indicate that the signature represents the immune infiltration of HCC.

The immune infiltration algorithms were based on the expression of cell lineage-specific genes, but the single-cell landscape is still vague. Single-cell sequencing data from 14 HCC samples (GSE156625) were used for analyses. After preprocessing and annotation ( Figure 6A ), the samples’ signature values were calculated using aggregated expression to generate pseudo bulk RNA-seq data. The signature values were calculated for 14 samples ( Figure 6B ), and these samples were divided into the high- and low-risk groups. Almost all cell types were distributed in the high- and low-risk groups ( Figure 6C ), but the cell proportions were different ( Figures 6D, E ). Importantly, immune activation cells, including Tregs and CD4+ and CD8+ T cells, were highly enriched in the low-risk group, while basal cells (Epithelial cells (ECs), fibroblasts, and hepatocytes) were more abundant in the high-risk group. The candidate genes used for signature construction were expressed in different cell types ( Figure 6F ). NKG7 was expressed in NK/NKT cells, ITM2A had the highest expression in CD8+ T cells, and DNASE1L3 was detected in ECs. These results indicate that the candidate genes originated from different cell types and that the signature reflected immune infiltration at the single-cell level.

Since the candidate genes were immune cell activation-related genes, the differences in cell–cell interactions reflected by the signature were investigated next. Overall, the interaction number and strength were similar between the high- and low-risk groups ( Figure 7A ), while the interactions between cell types were different ( Figure 7B ). The differential interactions among cell types were estimated and visualized in Figures 7C, D . The most altered interactions between groups were those involving fibroblasts, NKT cells, and myeloid cells. The differential pathways were also identified ( Figure 7E ). PECAM2, GAP, BAFF, and CDH5 were detected only in the low-risk samples, while COLLAGEN, FN1, NOTCH, THBS, ADGRA, VEGF, APP, and VTN signaling were detected specifically in the high-risk samples ( Figures 7F, G ). The differences in fibroblast interactions were mostly between mast cells and ECs, while the differences in NKT cell interactions were mostly between myeloid cells and mast cells ( Figure 7H ). Taken together, these results suggest that the signature reflects the status of immune cell–cell interactions, especially those involving fibroblasts, myeloid cells, and NK cells.

Since the signature reflected various HCC statuses and these statuses were significantly associated with drug response, the IC50 value of each drug for each cancer sample was evaluated using the “oncopredict” algorithm. Differences in IC50 values between the high-risk and low-risk groups for each drug were estimated. The low-risk samples were sensitive to most differential drugs ( Figures 8A, B ), while the high-risk samples were sensitive to only a few drugs. Correlations between candidate genes/the signature and drug IC50 values were evaluated, and NKG7 and ITM2A contributed the most to drug sensitivity ( Figure 8C ). Since it is difficult to alter the signature value, RORC was knocked down, and the drug sensitivity was quantified by estimating the IC50. As expected, the siRORC group was significantly more resistant to both paclitaxel and docetaxel in both the HepG2 and Huh7 cell lines ( Figures 8D, E ). Collectively, these results indicate that candidate gene RORC and the signature predict chemotherapy sensitivity, especially to paclitaxel and docetaxel.

To validate the roles of the candidate genes, the growth of the HCC cell lines HepG2 and Huh7 was assessed following the knockdown of RORC ( Figure 9A ). The growth rate was measured using the CCK-8 assay. Significant changes in the proliferation of both cancer cell lines were observed ( Figure 9B ). Consistently, the downregulation of RORC also enhanced the migration ability of both cell lines ( Figures 9C, D ), indicating the role of RORC in inhibiting the proliferation and migration phenotypes of HCC cell lines.

The RORC was reported as a transcription factor, while the detailed downstream genes are still not clear. To investigate the role of RORC in HCC, we analyzed several cell lines knocking down and/or overexpressing RORC, and we identified the differentially expressed genes in all these cell lines between control and overexpression groups ( Figure 10A ). Of these genes, CDC6 was the frequently reported prognostic gene in HCC. The expression of CDC6 was also significantly associated with survival in TCGA dataset ( Figures 10B, C ). To validate this finding, we knocked down RORC in the HepG2 and Huh7 cell lines and detected significantly increased CDC6 expression ( Figure 10D ). Additionally, ChIP-PCR revealed that RORC was enriched in the promoter region of CDC6 ( Figure 10E ), and the luciferase report indicated that the RORC promoter significantly decreased the luciferase signal intensity ( Figure 10F ). Since the role of CDC6 has been widely reported, we next investigated whether RORC inhibits proliferation via CDC6. After knocking down CDC6 in the siRORC group, the proliferation rate of the HepG2 and Huh7 cell lines was quantified using CCK-8 assay, as the previous result. CDC6 knockdown significantly restored the proliferation increase caused by RORC knockdown in both cell lines ( Figure 10F ). Taken together, these results indicate that RORC suppresses proliferation via CDC6 inhibition.

To further address the prognostic role of RORC and CDC6, we collected 100 primary HCC samples from our affiliation and quantified the expression value with immune histochemical staining microarray. As a result, higher expression of RORC resulted in a prolonged survival period compared to the RORC-low group ( Figure 11A ), while elevated CDC6 was correlated with a worse survival rate ( Figure 11B ), which is consistent with TCGA dataset ( Figures 11C , 10B, C ). It was noted the RORC was significantly downregulated in tumor tissues compared to the normal tissue ( Figure 11D ), while CDC6 was upregulated ( Figure 11E ). In addition, CDC6 and RORC were significantly associated with overall survival using Cox multivariate regression ( Figures 11F, G , p < 0.05). Collectively, these results indicate that RORC and its downstream target CDC6 serve as prognostic markers for HCC.