HOXD9 prioritization within developmental regulator-related genes suggests stage-adjusted prognostic relevance in hepatocellular carcinoma

Introduction

Hepatocellular carcinoma (HCC) is far more than just a daunting statistic in global cancer mortality; for clinicians and researchers alike, it remains a notoriously unpredictable disease (1,2). Even as treatment options ranging from systemic therapies to multimodal management continue to expand, patient outcomes remain highly variable (3). This is not merely a matter of clinical variation; it reflects the molecular heterogeneity that characterizes HCC.

In HCC, clinical decision-making still heavily relies on traditional tumor staging and histological features, which remain the cornerstone for prognosis and treatment planning (4,5). Despite advances in molecular and genetic profiling, these established pathological criteria provide essential information about tumor size, number, vascular invasion, and differentiation that correlate with patient outcomes (6,7). However, the heterogeneity of HCC in terms of etiology, tumor biology, and patient characteristics limits the predictive accuracy of staging systems alone (8). Emerging biomarkers and imaging techniques show promise but are yet to replace or significantly augment the prognostic power of conventional histopathology in routine clinical practice (9). Molecular biomarkers such as programmed death-ligand 1 (PD-L1) expression, tumor mutational burden, and immune cell profiles provide insights into tumor biology and immunotherapy response but are not yet standard for prognosis (10). Integrating new molecular insights with classical staging and histological assessment is essential for advancing personalized management of HCC (11). This integration aims to combine the proven robustness of traditional methods with the detailed biological information provided by emerging molecular markers, enhancing risk stratification and treatment tailoring (12,13). This highlights a specific biological niche that deserves further scrutiny: the developmental regulatory pathways. These are not just random markers; they are the original architects of cell fate and tissue patterning (14). The reawakening of embryonic gene programs, particularly those involving the homeobox (HOX) gene family, is increasingly recognized as a hallmark of aggressive cancer behavior (15). As master regulators of embryogenesis and cell fate determination, these developmental genes can become aberrantly reactivated in tumors, endowing cancer cells with stem cell-like properties, enhanced plasticity, and greater metastatic potential (16). This plasticity may allow tumor cells to adapt to environmental pressures, evade therapy, and contribute to tumor heterogeneity and progression (17). In HCC, such lineage-related dysregulation is emerging as a biologically meaningful feature of tumor progression (18). Rather than screening the entire transcriptome in an unconstrained manner, we therefore focused on a pre-specified developmental regulator-related candidate space to identify truly prioritized prognostic markers (19). This candidate space was defined according to developmental regulator-related gene-family annotation and biological rationale, rather than by visual inspection of The Cancer Genome Atlas liver hepatocellular carcinoma (TCGA-LIHC) differential expression results.

In this study, we did not aim to simply add another gene to the growing list of HCC-associated biomarkers. Instead, we established a stepwise and rigorous prioritization framework in the TCGA-LIHC cohort to identify developmental regulators with genuine clinical relevance. Through this process, HOXD9 emerged as the most compelling candidate. Beyond identifying HOXD9, we further asked whether it could provide incremental prognostic value beyond pathologic stage alone. By additionally validating the findings in the LIRI-JP (Liver Cancer-RIKEN, Japan) cohort, we sought to move beyond a purely descriptive bioinformatics analysis and assess whether a development-focused strategy could refine prognostic stratification in HCC in a clinically meaningful way. We present this article in accordance with the TRIPOD reporting checklist (available at https://tcr.amegroups.com/article/view/10.21037/tcr-2026-0825/rc).

Methods

Data sources and cohort construction

The TCGA-LIHC cohort was used as the primary discovery cohort. Raw RNA sequencing (RNA-seq) count data and corresponding clinical information were obtained from the Genomic Data Commons (GDC) portal. The full tissue set was used for the initial transcriptome-wide differential expression analysis, whereas a dedicated survival cohort was constructed for clinicopathological association analyses and prognostic modeling.

For external validation, the LIRI-JP cohort was obtained from International Cancer Genome Consortium (ICGC) resources. Donor-level expression and clinical data were cleaned and merged, and samples with the required expression, stage, and survival information were retained for external validation. The study was conducted in accordance with the Declaration of Helsinki and its subsequent amendments.

Differential expression analysis and visualization

Differential expression analysis between tumor and normal samples was performed using the differential expression analysis for sequence count data 2 (DESeq2) package in R. Differentially expressed genes (DEGs) were defined using an adjusted P value <0.05 and |log₂ fold change| >1. Principal component analysis (PCA), volcano plots, and heatmaps were used to visualize global transcriptional differences between tumor and normal tissues. DESeq2-derived results were also used to obtain gene-level differential expression statistics for subsequent candidate gene screening.

Pre-specified candidate gene set and stepwise prioritization strategy

To reduce arbitrary post hoc single-gene selection, we first defined a developmental regulator-related candidate gene set before inspection of the TCGA-LIHC volcano plot and before downstream survival modeling. The candidate genes were included based on gene-family annotation and developmental biology rationale, including HOX/HOX-related transcription factors and other developmental transcriptional regulators. The rationale for including each of the 11 genes is provided in Table S1. The TCGA-LIHC volcano plot was used only to visualize the differential expression status of these pre-specified candidates within the transcriptome-wide differential expression landscape and was not used to define the candidate list.

The pre-specified candidate genes were then evaluated through a stepwise framework integrating tumor vs. normal differential expression, clinicopathological association analysis, Kaplan-Meier survival analysis, and Cox regression. Genes showing more consistent evidence across these analytical layers were prioritized for further prognostic evaluation. Within this framework, HOXD9 and MESP2 emerged as the two leading candidate genes. HOXD9 was subsequently identified as the final independent prognostic core gene, whereas MESP2 was retained as a strong supporting gene.

Expression validation and clinicopathological association analysis

Candidate gene expression was further evaluated in relation to tumor-normal status and clinicopathological characteristics. We validated their expression by comparing tumor and normal tissues in the TCGA-LIHC cohort. Within the survival cohort, we compared expression levels across different pathologic stages and tumor grades using the Kruskal-Wallis test. Associations between candidate gene expression and vital status were evaluated using the Wilcoxon rank-sum test. For all these comparisons, we used variance-stabilizing transformation (VST)-normalized values derived from the RNA-seq count matrix to ensure data stability.

Survival analysis and Cox regression

We chose overall survival (OS) as our primary endpoint. Kaplan-Meier survival analysis was performed for visualization and descriptive survival stratification between high- and low-expression groups for each candidate gene. For these Kaplan-Meier plots, patients were stratified into high- and low-expression groups based on the median VST-normalized value of each gene and were compared using the log-rank test. The median cutoff was used as a simple, non-outcome-optimized threshold for visualization and was not used to define the main Cox-based prognostic models. For the survival screening step across the 11 pre-specified candidate genes, P values from Kaplan-Meier log-rank tests were adjusted using the Benjamini-Hochberg false discovery rate (FDR) method. Both nominal and FDR-adjusted P values were reported.

For the regression work, we started with a univariate Cox model to screen available candidate clinicopathological variables and key candidate genes, including age, sex, pathologic stage, tumor grade, HOXD9 expression, and MESP2 expression. The multivariable models were constructed using variables available with sufficient completeness in the TCGA-LIHC survival cohort. Gene expression values were modeled as continuous VST-normalized variables in univariate Cox regression, multivariable Cox regression, and prognostic model construction. For univariate Cox screening across the 11 pre-specified candidate genes, P values were also adjusted using the Benjamini-Hochberg FDR method. Any variable of interest was then entered into a multivariable Cox proportional hazards regression. This allowed us to identify independent prognostic factors, for which we reported hazard ratios (HRs), 95% confidence intervals (CIs), and two-sided P values. Survival analyses were performed using the survival package in R.

Prognostic model development and internal evaluation

Based on the multivariable Cox regression results, a main prognostic model integrating pathologic stage and HOXD9 expression was constructed. An extended model additionally including MESP2 expression was also evaluated. A nomogram based on the main model was constructed to estimate 1-year and 3-year OS probabilities.

Model performance was evaluated using calibration analysis and time-dependent receiver operating characteristic (ROC) analysis. Calibration analysis was used to compare predicted and observed survival probabilities, and time-dependent ROC analysis was used to calculate the area under the curve (AUC) at 1 and 3 years. These were performed using the rms and timeROC packages, respectively.

Formal comparison between the stage-only model and the stage plus HOXD9 model

A formal comparison between the stage-only model and the stage plus HOXD9 model was performed using the same complete-case dataset. We used the likelihood ratio test to compare model fit and Harrell’s concordance index (C-index) for discrimination. Bootstrap resampling with 2,000 iterations was used to estimate the 95% CI and P value for the difference in C-index. Time-dependent AUC values at 1 and 3 years were also compared to assess the incremental prognostic value of HOXD9 beyond pathologic stage.

External validation in the LIRI-JP cohort

External validation was performed in the LIRI-JP cohort. After cleaning and merging the donor-level data, we first re-evaluated HOXD9 and MESP2 using Kaplan-Meier analysis. The main external validation model, including tumor-node-metastasis (TNM) stage and HOXD9 expression, was evaluated using multivariable Cox regression and time-dependent ROC analysis. We also explored the extended model in this external set, reporting 1-year and 3-year AUCs.

Statistical analysis

All analyses were conducted in R. Differential expression analysis was performed using DESeq2. Survival analyses and Cox regression were performed using the survival package. Nomogram construction and calibration analysis were performed using rms, and time-dependent ROC analysis was performed using timeROC. Benjamini-Hochberg FDR adjustment was applied to the candidate-gene survival screening analyses as described above. Unless specified otherwise, all tests were two-sided, and we used P<0.05 as the threshold for statistical significance.

Results

The overall analytical workflow of the study is summarized in Figure 1. Briefly, TCGA-LIHC RNA-seq and clinical data were collected to construct tumor-normal cohorts for differential expression analysis and a survival cohort for prognostic evaluation. Developmental regulator-related candidate genes were then prioritized through a stepwise framework integrating tumor vs. normal expression differences, clinicopathological associations, Kaplan-Meier survival analysis, and univariate Cox regression. Candidate clinicopathological variables and key candidate genes were subsequently entered into multivariable Cox regression to identify independent prognostic factors. Based on the final screening results, a prognostic model integrating pathologic stage and HOXD9 expression was developed and assessed using nomogram construction, calibration analysis, time-dependent ROC analysis, and formal comparison with the stage-only model. External validation was further performed in the LIRI-JP cohort.

Figure 1 Overall study workflow of candidate gene prioritization, prognostic modeling, and external validation in LIHC. The study workflow is shown schematically. TCGA-LIHC RNA-seq and clinical data were collected to construct tumor-normal cohorts for differential expression analysis and a survival cohort for prognostic evaluation. Developmental regulator-related candidate genes were prioritized through a stepwise framework integrating tumor vs. normal differential expression, clinicopathological association analysis, Kaplan-Meier survival analysis, and univariate Cox regression. Candidate clinicopathological variables and key candidate genes were then entered into multivariable Cox regression to identify independent prognostic factors. A prognostic model integrating pathologic stage and HOXD9 expression was subsequently developed and assessed using nomogram construction, calibration analysis, time-dependent ROC analysis, and formal comparison with the stage-only model. External validation was further performed in the LIRI-JP cohort. DESeq2, differential expression analysis for sequence count data 2; LIHC, liver hepatocellular carcinoma; LIRI-JP, Liver Cancer-RIKEN, Japan; PCA, principal component analysis; RNA-seq, RNA sequencing; ROC, receiver operating characteristic; TCGA, The Cancer Genome Atlas.

Baseline characteristics and transcriptome-wide differential expression in the TCGA-LIHC cohort

A total of 370 patients were included in the TCGA-LIHC survival cohort. The median age was 61.0 years [interquartile range (IQR), 51.2–69.0 years], and 249 patients (67.3%) were male. Among these patients, 130 (35.1%) had died during follow-up. Pathologic stage I, II, III, and IV disease was observed in 168 (45.4%), 87 (23.5%), 86 (23.2%), and 6 (1.6%) cases, respectively, whereas 23 cases (6.2%) had unknown stage information. Tumor grade G1, G2, G3, and G4 accounted for 14.9%, 47.8%, 32.7%, and 3.2% of the cohort, respectively. The median OS time was 587.5 days (IQR, 327.2–1,085.0 days) (Table 1).

Transcriptome-wide differential expression analysis showed a clear separation between tumor and normal samples in the TCGA-LIHC cohort. PCA demonstrated distinct clustering of tumor and normal tissues, indicating substantial global transcriptional differences between the two groups (Figure 2A). Consistently, the heatmap of the top DEGs showed markedly different expression patterns between tumor and normal tissues (Figure 2B). The volcano plot further demonstrated a large number of significantly dysregulated genes. The pre-specified developmental regulator-related candidate genes were annotated on the volcano plot only to visualize their positions within the transcriptome-wide differential expression landscape. All 11 candidate genes, including HOXD9, MESP2, HOXA13, HOXA10, SIX2, SIX4, HOXD4, HOXD3, OSR2, IRX5, and VAX2, were located on the upregulated side of the distribution (Figure 2C).

Figure 2 PCA, volcano plot, and heatmap of differential expression analysis in the TCGA-LIHC cohort. (A) PCA showing the global transcriptional separation between tumor and normal tissues in the TCGA-LIHC cohort. (B) Heatmap of representative DEGs between tumor and normal tissues. (C) Volcano plot of transcriptome-wide differential expression analysis. Upregulated genes are shown on one side of the distribution and downregulated genes on the other. Pre-specified developmental regulator-related candidate genes are labeled to visualize their positions within the transcriptome-wide differential expression landscape. DEG, differentially expressed gene; PCA, principal component analysis; TCGA-LIHC, The Cancer Genome Atlas liver hepatocellular carcinoma.

Expression validation and clinicopathological associations of candidate developmental regulator-related genes

Expression validation in tumor and normal tissues showed that all 11 pre-specified developmental regulator-related candidate genes were significantly upregulated in LIHC tissues compared with normal liver tissues (all P<0.001; Figure 3, Table S2). Among these genes, HOXA13, SIX2, HOXA10, HOXD4, and HOXD9 showed particularly large positive log₂ fold changes, whereas IRX5 showed a smaller but still significant increase in tumor tissues (Table S2).

Figure 3 Expression validation of candidate developmental regulator-related genes in normal and tumor tissues of the TCGA-LIHC cohort. Boxplots showing the expression distributions of candidate developmental regulator-related genes in normal liver tissues and LIHC tumor tissues in the TCGA-LIHC cohort. Expression values are presented as VST-normalized values. P values were calculated for tumor vs. normal comparisons. All displayed genes were significantly differentially expressed between the two groups. LIHC, liver hepatocellular carcinoma; TCGA-LIHC, The Cancer Genome Atlas liver hepatocellular carcinoma; VST, variance stabilizing transformation.

Clinicopathological association analysis showed that all candidate genes were significantly associated with pathologic stage and tumor grade. For pathologic stage, all 11 genes were significant, with P values ranging from <0.001 to 0.002. For tumor grade, all 11 genes were also significant, and all P values were <0.001 (Table 2). In contrast, associations with vital status were more selective. HOXD9 expression was significantly associated with vital status (P<0.001), and MESP2 and IRX5 also showed significant differences between alive and dead patients (P=0.004 and P=0.045, respectively), whereas the remaining candidate genes were not significantly associated with vital status (Table 2). These results indicated that developmental regulator-related candidate genes were broadly associated with tumor progression features, while HOXD9 and MESP2 showed stronger links with survival status.

Survival screening supported HOXD9 as the most robust prognostic candidate gene

Kaplan-Meier survival analysis showed that patients with high HOXD9 expression had shorter OS than those with low HOXD9 expression (nominal P<0.001; Figure 4A). A similar pattern was observed for MESP2, with high MESP2 expression also associated with shorter OS (nominal P=0.005; Figure 4B). These high- and low-expression groups were defined by the median VST-normalized expression value and were used for Kaplan-Meier visualization and descriptive survival stratification only. Subsequent Cox regression and prognostic modeling treated gene expression as continuous VST-normalized variables.

Figure 4 Kaplan-Meier survival curves of HOXD9 and MESP2 in the TCGA-LIHC survival cohort. (A) Kaplan-Meier OS curves stratified by HOXD9 expression. (B) Kaplan-Meier OS curves stratified by MESP2 expression. High- and low-expression groups were defined using the median VST-normalized expression value of each gene. Median-based grouping was used for Kaplan-Meier visualization and descriptive survival stratification only, and was not used to define the Cox-based prognostic models. P values shown in the figure are nominal log-rank P values. FDR-adjusted survival screening results are provided in Table S3. FDR, false discovery rate; OS, overall survival; TCGA-LIHC, The Cancer Genome Atlas liver hepatocellular carcinoma; VST, variance stabilizing transformation.

The broader candidate-gene survival summary showed that HOXD9 had the strongest and most consistent prognostic signal among the candidate genes. In univariate Cox regression, higher HOXD9 expression was associated with worse OS and remained significant after Benjamini-Hochberg correction (HR =1.212, 95% CI: 1.091–1.346; nominal P<0.001; FDR =0.003811). HOXD9 also remained significant in Kaplan-Meier screening after multiple-testing correction (nominal P<0.001; FDR =0.006689). MESP2 showed a significant Kaplan-Meier association after FDR correction (nominal P=0.005; FDR =0.030039), whereas its univariate Cox association did not remain significant after FDR adjustment (HR =1.204, 95% CI: 1.036–1.398; nominal P=0.02; FDR =0.084246). IRX5 showed nominal Kaplan-Meier significance before correction (nominal P=0.03), but this association did not remain significant after FDR adjustment (FDR =0.111633), and its univariate Cox result was not significant (HR =1.156, 95% CI: 0.992–1.347; nominal P=0.06; FDR =0.234950). The remaining candidate genes did not show significant survival associations after multiple-testing correction (Table S3). Based on these results, HOXD9 was prioritized as the most robust survival-associated candidate gene, while MESP2 was retained as a supporting candidate for multivariable evaluation.

Multivariable Cox regression identified pathologic stage and HOXD9 expression as independent prognostic factors

To determine whether the key candidate genes retained prognostic value after adjustment for available clinicopathological variables, multivariable Cox regression analysis was performed including pathologic stage, tumor grade, age, sex, HOXD9 expression, and MESP2 expression. In this model, only pathologic stage (HR =1.565, 95% CI: 1.266–1.935; P<0.001) and HOXD9 expression (HR =1.165, 95% CI: 1.035–1.312; P=0.01) remained significantly associated with OS. In contrast, MESP2 expression (HR =1.122, 95% CI: 0.953–1.321; P=0.17), tumor grade (HR =1.191, 95% CI: 0.912–1.556; P=0.20), age (HR =1.008, 95% CI: 0.993–1.022; P=0.31), and sex (HR =0.904, 95% CI: 0.612–1.336; P=0.61) were not significant (Table 3, Figure 5).

Figure 5 Forest plot of multivariable Cox regression analysis in the TCGA-LIHC survival cohort. Forest plot showing HRs, 95% CIs, and P values from the multivariable Cox proportional hazards model. Variables included pathologic stage, HOXD9 expression, MESP2 expression, tumor grade, age, and sex. Pathologic stage and HOXD9 expression remained significantly associated with OS in the final multivariable model. CI, confidence interval; HR, hazard ratio; OS, overall survival; TCGA-LIHC, The Cancer Genome Atlas liver hepatocellular carcinoma.

The corresponding univariate and multivariable Cox summary also showed that pathologic stage and HOXD9 expression were the only retained independent factors in the final model, whereas MESP2 lost significance after adjustment (Table S3). These findings defined pathologic stage and HOXD9 expression as the components of the final main prognostic model.

Construction and performance of the main prognostic model

A nomogram integrating pathologic stage and HOXD9 expression was constructed to estimate 1-year and 3-year OS probabilities in the TCGA-LIHC cohort (Figure 6). The main prognostic model yielded a 1-year AUC of 0.707 and a 3-year AUC of 0.701. By comparison, the extended model incorporating pathologic stage, HOXD9 expression, and MESP2 expression yielded a 1-year AUC of 0.714 and a 3-year AUC of 0.699 (Table 4).

Figure 6 Nomogram integrating pathologic stage and HOXD9 expression for prediction of OS in the TCGA-LIHC cohort. Nomogram constructed from the main prognostic model integrating pathologic stage and HOXD9 expression in the TCGA-LIHC cohort. The nomogram was developed to estimate 1-year and 3-year OS probabilities. OS, overall survival; TCGA-LIHC, The Cancer Genome Atlas liver hepatocellular carcinoma.

Calibration analysis of the main model showed general agreement between predicted and observed survival probabilities at both 1 year and 3 years (Figure 7A,7B). Time-dependent ROC analysis showed AUC values of 0.707 at 1 year and 0.701 at 3 years for the main model (Figure 7C,7D), consistent with the values summarized in Table 4.

Figure 7 Calibration and time-dependent ROC curves of the main prognostic model in the TCGA-LIHC cohort. (A) Calibration curve for 1-year OS in the main prognostic model. (B) Calibration curve for 3-year OS in the main prognostic model. (C) Time-dependent ROC curve for 1-year OS in the main prognostic model. (D) Time-dependent ROC curve for 3-year OS in the main prognostic model. The main prognostic model included pathologic stage and HOXD9 expression. AUC, area under the curve; OS, overall survival; ROC, receiver operating characteristic; TCGA-LIHC, The Cancer Genome Atlas liver hepatocellular carcinoma.

Training-cohort analysis suggested incremental prognostic value of HOXD9 beyond pathologic stage

To further assess whether HOXD9 improved prognostic stratification beyond pathologic stage alone, a main Cox model integrating pathologic stage and HOXD9 expression was evaluated separately. In this model, both pathologic stage (HR =1.530, 95% CI: 1.249–1.874; P<0.001) and HOXD9 expression (HR =1.186, 95% CI: 1.056–1.332; P=0.004) were significantly associated with OS (Table 5).

Formal comparison between the stage-only model and the stage plus HOXD9 model showed statistically detectable improvement after inclusion of HOXD9. The C-index increased from 0.601 to 0.651, corresponding to a ΔC-index of 0.050 (95% CI: 0.026–0.099; P=0.004). The 1-year AUC increased from 0.656 to 0.707, with a ΔAUC of 0.051 (P=0.02). The 3-year AUC increased from 0.659 to 0.701, with a ΔAUC of 0.042 (P=0.06). Model fit also improved, as indicated by the likelihood ratio test (χ²=8.149, P=0.004) (Table 6). These results indicated statistically significant but numerically modest improvement in model discrimination after adding HOXD9 to pathologic stage. The detailed comparison results are provided in Table S4.

External validation in the LIRI-JP cohort

External validation was performed in the LIRI-JP cohort. A total of 117 analysis rows corresponding to 117 unique donors were available, with 106 cases having non-missing HOXD9 expression and 114 cases having non-missing MESP2 expression. TNM stage and OS data were complete for all 117 cases. Baseline clinicopathological characteristics of the LIRI-JP external validation cohort are summarized in Table S5. The median age was 69.0 years (IQR, 61.0–75.0), and 87 patients (74.4%) were male. TNM stage I, II, III, and IV disease was observed in 18 (15.4%), 56 (47.9%), 34 (29.1%), and 9 (7.7%) cases, respectively. The median OS/follow-up time was 1,050.0 days (IQR, 750.0–1,290.0 days). Etiologic background was not available in the analyzed donor-level data. Kaplan-Meier analysis in the external cohort showed that HOXD9 expression significantly stratified OS (P=0.01; Figure 8). In contrast, MESP2 did not show a significant Kaplan-Meier difference in the external cohort (P=0.22; Table S6).

Figure 8 Kaplan-Meier survival curve of HOXD9 in the LIRI-JP external validation cohort. Kaplan-Meier OS curves stratified by HOXD9 expression in the external LIRI-JP cohort. High- and low-expression groups were defined according to the prespecified expression grouping strategy used in the external validation analysis. P values were calculated using the log-rank test. LIRI-JP, Liver Cancer-RIKEN, Japan; OS, overall survival.

For the main external validation model including TNM stage and HOXD9 expression, TNM stage remained significantly associated with OS (HR =2.453, 95% CI: 1.452–4.142; P<0.001), whereas HOXD9 expression was not significant after adjustment (HR =1.027, 95% CI: 0.763–1.383; P=0.86) (Table 7). The main external validation model yielded a 1-year AUC of 0.819 and a 3-year AUC of 0.739 (Table S6). In the extended external validation model including TNM stage, HOXD9 expression, and MESP2 expression, TNM stage remained significant (HR =2.371, 95% CI: 1.391–4.041; P=0.002), whereas HOXD9 and MESP2 were not significant (Tables S6,S7).

Discussion

HCC is characterized by substantial clinical and molecular heterogeneity, which continues to limit prognostic stratification based on conventional clinicopathological factors alone (20,21). In the present study, developmental regulator-related candidate genes were systematically screened in TCGA-LIHC, and HOXD9 and MESP2 emerged as the two most consistently supported candidates (22). Among them, only HOXD9 remained independently associated with OS after adjustment for clinicopathological variables in the TCGA-LIHC cohort. On this basis, a prognostic model integrating pathologic stage and HOXD9 expression was established (23). This model showed moderate discriminatory ability and performed better than pathologic stage alone in the training cohort. However, external validation in the LIRI-JP cohort did not reproduce the independent stage-adjusted association of HOXD9, although HOXD9 retained survival-stratification ability in Kaplan-Meier analysis.

The candidate gene set was defined on the basis of developmental regulator-related gene-family annotation and biological rationale before inspection of the TCGA-LIHC volcano plot and before downstream prognostic modeling (24). Thus, Figure 2C was used to display the differential expression positions of the pre-specified candidates, rather than to define the candidate list. In this study, candidate genes were prioritized through a stepwise framework that combined tumor vs. normal differential expression, clinicopathological association, Kaplan-Meier analysis, and Cox regression. Importantly, the median-based high- and low-expression grouping was retained only for Kaplan-Meier visualization and descriptive survival stratification, whereas the main prognostic inference was based on Cox regression and prognostic modeling using continuous VST-normalized expression values. Within this process, HOXD9 showed the most consistent overall pattern across analytical layers and remained significant in both Kaplan-Meier and univariate Cox survival screening after Benjamini-Hochberg correction. MESP2 retained a supportive but weaker profile, with FDR-adjusted significance in Kaplan-Meier analysis but not in univariate Cox regression. This association also did not persist in the final multivariable model. This distinction is important for the interpretation of the current findings. HOXD9 should be viewed as the final independent core prognostic gene in the TCGA-LIHC cohort, whereas MESP2 is more appropriately regarded as a strong supporting gene rather than an equivalent endpoint marker (25).

The current findings also need to be interpreted in the context of existing literature. HOXD9 has already appeared in multigene prognostic signatures in HCC/LIHC-related studies, which means that its presence in liver cancer transcriptomic analyses is not entirely unprecedented (26). The value of the current study, therefore, does not lie in presenting HOXD9 as a completely novel observation, but in showing that HOXD9 remained the most consistently prioritized candidate within a developmental regulator-focused screening framework and retained significance after multivariable adjustment (27). In addition, the present work extends beyond simple candidate identification by formally assessing whether HOXD9 contributes prognostic information beyond pathologic stage. This point is particularly relevant because many public-dataset prognostic studies report statistically significant genes without directly testing whether they improve risk discrimination over an established clinicopathological baseline (28).

This more formal model comparison represents one of the main strengths of the study. After HOXD9 was added to pathologic stage, the likelihood ratio test showed improved model fit, the C-index increased, and the 1-year time-dependent AUC was significantly higher than that of the stage-only model. The 3-year AUC also increased numerically, although the difference did not reach conventional statistical significance. Taken together, these findings indicate that HOXD9 contributed additional prognostic information beyond pathologic stage alone in the TCGA-LIHC cohort. However, the absolute improvement in discrimination was modest, with ΔAUC values of 0.051 at 1 year and 0.042 at 3 years. Therefore, the statistical significance of these increments should not be equated with established clinical utility. At the same time, the overall AUC values remained in a moderate range. For this reason, the model should be interpreted as an exploratory and parsimonious risk-stratification model rather than a stand-alone high-accuracy clinical prediction tool. This interpretation is further supported by the absence of several HCC-specific clinical variables from the model, including liver functional reserve, etiologic background, alpha-fetoprotein (AFP), vascular invasion, and treatment-related factors. Whether this modest improvement is sufficient to alter clinical decision-making requires further validation in larger cohorts with richer clinical annotation. This more restrained interpretation is more consistent with the actual performance observed in the training cohort.

The external validation results also merit cautious interpretation. In the LIRI-JP cohort, HOXD9 significantly stratified OS by Kaplan-Meier analysis, and the main model achieved acceptable 1-year and 3-year AUC values. However, HOXD9 did not remain independently significant after adjustment for TNM stage in the external multivariable model (29). The adjusted HR for HOXD9 was close to null, indicating that the independent prognostic signal observed in TCGA-LIHC was not reproduced in LIRI-JP. A similar pattern was observed for MESP2, which also failed to retain significance in the extended external model. Accordingly, the incremental prognostic value of HOXD9 beyond pathologic stage should be interpreted as a training-cohort finding rather than as externally confirmed evidence (30). This pattern may reflect differences in cohort composition, staging structure, clinical annotation, and sample size across public datasets (31). The added baseline summary of the LIRI-JP cohort provides clinical context for interpreting the external validation results, including the cohort’s age distribution, sex distribution, TNM stage structure, follow-up time, and incomplete etiologic annotation. Therefore, the LIRI-JP analysis supports the unadjusted survival-stratification relevance of HOXD9 but does not validate its independent prognostic contribution beyond stage (32).

Several strengths of this study should be noted. First, candidate gene prioritization was performed through a layered analytical framework rather than unrestricted single-gene selection from the full DEG landscape. Second, the study integrated clinicopathological association analysis, survival analysis, multivariable adjustment, prognostic model development, and formal incremental-value testing into a single workflow. Third, an independent external cohort was used for validation. These aspects improved the internal consistency of the analysis and reduced the arbitrariness that often weakens purely descriptive biomarker studies. Nonetheless, the present study also has important limitations. It was retrospective and based on public datasets, and the range of available clinical covariates was limited. Several clinically important prognostic variables in HCC, including liver functional reserve measures such as Child-Pugh score or Model for End-stage Liver Disease (MELD), etiologic background, AFP level, vascular invasion, and treatment information, were not consistently available and therefore could not be incorporated into the final model. As a result, the observed association between HOXD9 expression and OS may be affected by residual confounding from unmeasured clinical factors, and the pathologic stage plus HOXD9 model should not be interpreted as a complete real-world clinical prediction model. In addition, the study was designed to evaluate transcriptomic prognostic associations and model performance rather than biological mechanism, and no experimental validation was performed. Although median-based high- and low-expression grouping was used only for Kaplan-Meier visualization, this cutoff should not be interpreted as a validated biological or clinical threshold. Finally, although HOXD9 showed incremental prognostic value beyond pathologic stage in TCGA-LIHC, this independent stage-adjusted signal was not reproduced in the LIRI-JP multivariable model. This limits the strength of the external validation and indicates that larger, better-annotated cohorts are required before HOXD9 can be considered an externally validated independent prognostic biomarker.

In summary, this study identified HOXD9 as the most consistently supported developmental regulator-related prognostic candidate in TCGA-LIHC. HOXD9, together with pathologic stage, formed a concise prognostic model with moderate discriminatory performance in the TCGA-LIHC training cohort. Formal comparison with the stage-only model further showed that HOXD9 added incremental prognostic information within this cohort. In the LIRI-JP cohort, HOXD9 retained Kaplan-Meier survival-stratification ability, but its independent association with OS was not reproduced after TNM-stage adjustment. Overall, these findings support HOXD9 as a candidate biomarker for further prognostic evaluation, rather than an externally confirmed independent prognostic marker, and justify further biological and clinical validation.

Conclusions

In conclusion, this study systematically prioritized developmental regulator-related candidate genes in HCC and identified HOXD9 as the most consistently supported independent prognostic candidate in the TCGA-LIHC cohort. HOXD9, together with pathologic stage, formed a concise prognostic model with moderate discriminatory ability in the TCGA-LIHC training cohort, and formal model comparison showed that HOXD9 provided incremental prognostic information beyond pathologic stage alone within this cohort. In the LIRI-JP external validation cohort, HOXD9 stratified OS in Kaplan-Meier analysis but did not remain independently associated with OS after TNM-stage adjustment. Overall, HOXD9 may represent a candidate biomarker requiring further validation for prognostic stratification in HCC, rather than an externally confirmed independent prognostic marker at the current stage.

Acknowledgments

The authors acknowledge The Cancer Genome Atlas (TCGA) and the International Cancer Genome Consortium (ICGC) for providing publicly available data resources.

Footnote

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

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

Funding: This study was supported by the Project of Ningbo Leading Medical & Health Discipline (Traditional Chinese Medicine) (No. 2026-Z13).

Conflicts of Interest: All authors have completed the ICMJE uniform disclosure form (available at https://tcr.amegroups.com/article/view/10.21037/tcr-2026-0825/coif). All authors report funding support for the present work from the Project of Ningbo Leading Medical & Health Discipline (Traditional Chinese Medicine) (No. 2026-Z13). The authors have no other 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. Yang JD, Hainaut P, Gores GJ, et al. A global view of hepatocellular carcinoma: trends, risk, prevention and management. Nat Rev Gastroenterol Hepatol 2019;16:589-604. [Crossref] [PubMed]
2. Mauro E, de Castro T, Zeitlhoefler M, et al. Hepatocellular carcinoma: Epidemiology, diagnosis and treatment. JHEP Rep 2025;7:101571. [Crossref] [PubMed]
3. Cabibbo G, Enea M, Attanasio M, et al. A meta-analysis of survival rates of untreated patients in randomized clinical trials of hepatocellular carcinoma. Hepatology 2010;51:1274-83. [Crossref] [PubMed]
4. Amin MB, Greene FL, Edge SB, et al. The Eighth Edition AJCC Cancer Staging Manual: Continuing to build a bridge from a population-based to a more "personalized" approach to cancer staging. CA Cancer J Clin 2017;67:93-9.
5. Mauro E, Forner A. Barcelona Clinic Liver Cancer 2022 update: Linking prognosis prediction and evidence-based treatment recommendation with multidisciplinary clinical decision-making. Liver Int 2022;42:488-91. [Crossref] [PubMed]
6. Kreitmann L, D'Souza G, Miglietta L, et al. A computational framework to improve cross-platform implementation of transcriptomics signatures. EBioMedicine 2024;105:105204. [Crossref] [PubMed]
7. Shinkawa H, Tanaka S, Kabata D, et al. The Prognostic Impact of Tumor Differentiation on Recurrence and Survival after Resection of Hepatocellular Carcinoma Is Dependent on Tumor Size. Liver Cancer 2021;10:461-72. [Crossref] [PubMed]
8. Chan YT, Zhang C, Wu J, et al. Biomarkers for diagnosis and therapeutic options in hepatocellular carcinoma. Mol Cancer 2024;23:189. [Crossref] [PubMed]
9. Echle A, Rindtorff NT, Brinker TJ, et al. Deep learning in cancer pathology: a new generation of clinical biomarkers. Br J Cancer 2021;124:686-96. [Crossref] [PubMed]
10. Zhang N, Yang X, Piao M, et al. Biomarkers and prognostic factors of PD-1/PD-L1 inhibitor-based therapy in patients with advanced hepatocellular carcinoma. Biomark Res 2024;12:26. [Crossref] [PubMed]
11. Bruix J, Han KH, Gores G, et al. Liver cancer: Approaching a personalized care. J Hepatol 2015;62:S144-56. [Crossref] [PubMed]
12. Reig M, Forner A, Rimola J, et al. BCLC strategy for prognosis prediction and treatment recommendation: The 2022 update. J Hepatol 2022;76:681-93. [Crossref] [PubMed]
13. Cheng CCY, Cheung MF, Lee AY, et al. Multi-omic analysis of hepatocellular carcinoma reveals aberrant cis-regulatory changes and dysregulated retrotransposons with prognostic potential. Commun Biol 2025;8:1792. [Crossref] [PubMed]
14. Segal J, Cronk J, Ball B, et al. Parallels in Canonical Developmental Signaling Pathways between Normal Development and the Tumor Microenvironment. Cold Spring Harb Perspect Med 2025;a041609.
15. Li B, Huang Q, Wei GH. The Role of HOX Transcription Factors in Cancer Predisposition and Progression. Cancers (Basel) 2019;11:528. [Crossref] [PubMed]
16. Ben-Porath I, Thomson MW, Carey VJ, et al. An embryonic stem cell-like gene expression signature in poorly differentiated aggressive human tumors. Nat Genet 2008;40:499-507. [Crossref] [PubMed]
17. Bhat GR, Sethi I, Sadida HQ, et al. Cancer cell plasticity: from cellular, molecular, and genetic mechanisms to tumor heterogeneity and drug resistance. Cancer Metastasis Rev 2024;43:197-228. [Crossref] [PubMed]
18. Li MM, Tang YQ, Gong YF, et al. Development of an oncogenic dedifferentiation SOX signature with prognostic significance in hepatocellular carcinoma. BMC Cancer 2019;19:851. [Crossref] [PubMed]
19. Meng Q, Zhou Q, Chen X, et al. Prognostic hub gene CBX2 drives a cancer stem cell-like phenotype in HCC revealed by multi-omics and multi-cohorts. Aging (Albany NY) 2023;15:12817-51. [Crossref] [PubMed]
20. Ahn JC, Teng PC, Chen PJ, et al. Detection of Circulating Tumor Cells and Their Implications as a Biomarker for Diagnosis, Prognostication, and Therapeutic Monitoring in Hepatocellular Carcinoma. Hepatology 2021;73:422-36. [Crossref] [PubMed]
21. Chen J, Kaya NA, Zhang Y, et al. A multimodal atlas of hepatocellular carcinoma reveals convergent evolutionary paths and 'bad apple' effect on clinical trajectory. J Hepatol 2024;81:667-78. [Crossref] [PubMed]
22. Wang J, Ding ZW, Chen K, et al. A predictive and prognostic model for hepatocellular carcinoma with microvascular invasion based TCGA database genomics. BMC Cancer 2021;21:1337. [Crossref] [PubMed]
23. Wang Z, Zhu J, Liu Y, et al. Development and validation of a novel immune-related prognostic model in hepatocellular carcinoma. J Transl Med 2020;18:67. [Crossref] [PubMed]
24. Toshner M, Dunmore BJ, McKinney EF, et al. Transcript analysis reveals a specific HOX signature associated with positional identity of human endothelial cells. PLoS One 2014;9:e91334. [Crossref] [PubMed]
25. Ahmadi SE, Rahimi S, Zarandi B, et al. MYC: a multipurpose oncogene with prognostic and therapeutic implications in blood malignancies. J Hematol Oncol 2021;14:121. [Crossref] [PubMed]
26. Chen C, Liu YQ, Qiu SX, et al. Five metastasis-related mRNAs signature predicting the survival of patients with liver hepatocellular carcinoma. BMC Cancer 2021;21:693. [Crossref] [PubMed]
27. Dong S, Wang R, Wang H, et al. HOXD-AS1 promotes the epithelial to mesenchymal transition of ovarian cancer cells by regulating miR-186-5p and PIK3R3. J Exp Clin Cancer Res 2019;38:110. [Crossref] [PubMed]
28. Wang T, She Y, Yang Y, et al. Radiomics for Survival Risk Stratification of Clinical and Pathologic Stage IA Pure-Solid Non-Small Cell Lung Cancer. Radiology 2022;302:425-34. [Crossref] [PubMed]
29. van Leeuwen FD, Steyerberg EW, van Klaveren D, et al. Instability of the AUROC of Clinical Prediction Models. Stat Med 2025;44:e70011. [Crossref] [PubMed]
30. Wasenang W, Chaiyarit P, Proungvitaya S, et al. Serum cell-free DNA methylation of OPCML and HOXD9 as a biomarker that may aid in differential diagnosis between cholangiocarcinoma and other biliary diseases. Clin Epigenetics 2019;11:39. [Crossref] [PubMed]
31. Luo S, Liu L, Sun Y, et al. Spatial heterogeneity reveals an evolutionary signature predicting therapeutic response and clinical outcomes in hepatocellular carcinoma. Front Bioinform 2025;5:1669236. [Crossref] [PubMed]
32. Collins GS, Dhiman P, Ma J, et al. Evaluation of clinical prediction models (part 1): from development to external validation. BMJ 2024;384:e074819. [Crossref] [PubMed]