The HuH-7 and PLC-PRF-5 human liver cancer cell lines were obtained from FuHeng Biology. HuH-7 cells were cultured in DMEM (cat. no. C11995500BT; Gibco; Thermo Fisher Scientific, Inc.) and PLC-PRF-5 cells were cultured in Minimum Essential Medium (MEM; cat. no. C11095500BT; Gibco; Thermo Fisher Scientific, Inc.). Cell media were supplemented with 10% FBS (cat. no. FBS500-S; Ausgenex Pty, Ltd.) and 1% penicillin/streptomycin (cat. no. 15140122; Gibco; Thermo Fisher Scientific, Inc.). Cells were kept at 37˚C in a humidity-controlled incubator with 5% CO₂ supplementation. All cell lines underwent rigorous verification for mycoplasma contamination and were confirmed to be free of contamination. The cells were authenticated using short tandem repeat analysis.

siRNA transfection was conducted using Lipofectam ine^™ 2000 transfection reagent (cat. no. 52887; Invitrogen; Thermo Fisher Scientific, Inc.). All siRNA transfections were performed for 48 h at room temperature. HuH-7 and PLC-PRF-5 cells were cultured in 6-well plates to allow for adhesion and proliferation overnight. Before the introduction of the transfection agent, 1 ml serum-free DMEM/MEM was used to replace the medium. The scrambled control, CCCTC binding factor (CTCF) and nuclear respiratory factor 1 (NRF1) siRNA molecules were chemically synthesized by Nanjing GenScript Biotech Co., Ltd. For 6-well plates, 10 µl siRNA (20 µM) was dissolved in 5 µl siRNA transfection reagent and incubated for 5 min at room temperature. The aforementioned mixture was then introduced to the cells. After 12 h, 1 ml DMEM/MEM supplemented with 10% FBS was added to each well. Cells were seeded at a density of 6x10⁵ cells/well in 6-well culture plates. After 24 h, the cells were transfected with the siRNA for 48 h and then RNA was extracted. The siRNA sequences used were as follows: NRF1, 5'-GGAAACUUCGAGCCACGUU-3'; CTCF, 5'-GCGAAAGCAGCAUUCCUAUAU-3'; and control, 5'-UUCUCCGAACGUGUCACGU-3'.

Total RNA was extracted from HCC cells using RNeasy Kits (cat. no. 74104; Qiagen GmbH) according to the manufacturer's instructions. cDNA synthesis was conducted with the Primescript RT-reagent kit (cat. no. RR047A; Takara Bio, Inc.). The temperature and duration of reverse transcription were: 37˚C for 15 min and 85˚C for 5 sec. qPCR was performed using SYBR Premix Ex Taq (Takara Biotechnology Co., Ltd.) on an ABI7500 system (Applied Biosystems; Thermo Fisher Scientific, Inc.). The 2^-ΔΔCq method was utilized for data analysis, using β-actin for normalization (19,20). The thermocycling conditions were as follows: 95˚C for 5 min; 40 cycles of 95˚C for 15 sec and 60˚C for 30 sec; 1 cycle of 95˚C for 15 sec, 60˚C for 60 sec and 95˚C for 15 sec. The primer sequences used for qPCR were as follows: NRF1 forward, 5'-CCGGAAGAGGCAACAAACAC-3' and reverse, 5'-CTTGCTGTCCCACACGAGTAGT-3'; CTCF forward, 5'-CATCCAGCATCAGAAGTCACACA-3' and reverse, 5'-GCCTCTCCTGTCTACAAGCGTAA-3'; PAXIP1 forward, 5'-CCAGCTGTACGGACACTGAGG-3' and reverse, 5'-TTGTATGTCCCTGCTGGCTGT-3'; and β-actin forward, 5'-CACTCTTCCAGCCTTCCTTC-3' and reverse, 5'-GTACAGGTCTTTGCGGATGT-3'.

ChIP-seq data for PAXIP1, MYB proto-oncogene like 2 (MYBL2) and FOXO1 were retrieved from the ChIP-Atlas database (http://chip-atlas.org/). The cutoff of broad peak call was q<1x10^-5 (transcription start site ±1 kb). These datasets (GSE32465 and GSE104247) (21,22) underwent further analysis using the ChIP-seq pipeline, mainly using the open-source BEDTools and deepTools software suites. BEDTools (version 2.29.2) (23) was employed for the genome arithmetic. The computeMatrix program within deepTools (version 3.4.3) (24) facilitated the calculation of scores across genome regions and generated intermediate files for subsequent visualization with plotHeatmap in deepTools. For genome annotation, R (version 4.1.0; http://www.r-project.org/) was used to analyze ChIPseeker (25). Gene Ontology (GO) analysis was implemented using Database for Annotation, Visualization and Integrated Discovery (DAVID) functional annotation tools (https://david.ncifcrf.gov/) (26). Metascape (https://metascape.org/) (27) was also employed to perform DisGeNET and PaGenBase enrichment analyses for PAXIP1 target genes.

HCCDB (http://lifeome.net/database/hccdb/) is a comprehensive HCC expression atlas encompassing 15 publicly available HCC gene expression datasets, which collectively include 3,917 samples (28). This repository integrates data from prominent sources such as Gene Expression Omnibus, The Cancer Genome Atlas (TCGA) Liver HCC Project (TCGA-LIHC) and Liver Cancer-RIKEN, Japan Project from the International Cancer Genome Consortium. The HCCDB provides a platform for visualizing outcomes of various computational analyses, including differential expression analysis, as well as tissue- and tumor-specific expression assessments. The 15 HCC datasets were searched with the keyword ‘PAXIP1’.

The mRNA and protein expression profiles were downloaded from the Cancer Cell Line Encyclopedia (https://depmap.org/portal/interactive/) and the HPA (https://www.proteinatlas.org/) (29,30). cBioPortal (https://www.cbioportal.org/), an open-source cancer genomics data platform, was used to analyze the mutations, copy-number alterations and gene expression of PAXIP1 in patients with HCC (31,32). Liver studies were selected and the keyword ‘PAXIP1’ was searched on the query page of the cBioPortal website (33-42).

Analysis of PAXIP1 expression in TCGA LIHC cohort was performed using the LinkedOmics database (http://www.linkedomics.org/) (43). Statistical analysis of PAXIP1 co-expression was conducted using Pearson's correlation coefficient, with results visualized using volcano plots, heatmaps and scatter plots. The functional module of the platform facilitates the examination of GO biological processes, Kyoto Encyclopedia of Genes and Genomes (KEGG) pathways, kinase-target enrichment, microRNA (miRNA)-target enrichment and TF-target enrichment via gene set enrichment analysis (GSEA).

KM survival analysis and plotting were performed using the R packages survival and survminer (https://CRAN.R-project.org/package=survminer). Analysis was performed using default parameters. The area under the curve (AUC) was analyzed using the R package timeROC (https://CRAN.R-project.org/package=timeROC). Based on prognostic clinical indicators and the survival analysis of the Cox regression model, age and pT_stage were entered into the risk model. The points against each factor were counted, and 1-, 3- and 5-year survival rates were also calculated. The nomogram was constructed using the rms package (https://CRAN.R-project.org/package=rms) (44). Additionally, the risk score was calculated as follows: Risk score=0.278 x PAXIP1 + 0.1718 x MYBL2 + 0.0175 x NRF1-0.2226 x FOXO1. Based on the risk score, patients with HCC were divided into the low-risk group and the high-risk group using the median risk score as the cutoff (45).

The human TF (hTF) target database (http://bioinfo.life.hust.edu.cn/hTFtarget#!/) represents an extensive resource dedicated to the regulation of hTFs and their respective targets (46). In the present study, this database was used to predict potential upstream TFs of PAXIP1. These TFs were then ranked according to their R- and P-value, providing a systematic evaluation of their potential regulatory roles. Find Individual Motif Occurrences (v4.10.0; https://meme-suite.org/meme/tools/fimo) was used to scan both the test and control sets, and then the numbers of recurrent motifs within the two sets were used to evaluate the significance of motifs for the TF (t-test with Bonferroni correction P<0.01) (46).

The TIMER database (https://cistrome.shinyapps.io/timer/) was used to analyze the expression profile of PAXIP1 and immune cell presence in HCC. For gene expression levels, log₂ transformed transcripts per million values were used.

UALCAN (http://ualcan.path.uab.edu) uses level 3 RNA-sequencing and clinical data from TCGA data of HCC. This platform facilitates a comprehensive analysis of gene expression, comparing tumor samples with healthy control tissues and examining variations across diverse tumor subgroups classified by cancer stage, tumor grade or other clinicopathological parameters. UALCAN was used to examine mRNA expression levels according to the online instructions. GEPIA (http://gepia.cancer-pku.cn/) was used to investigate MYBL2 and FOXO1 expression in LIHC.

Integrative Genomics Viewer (v2.17.0; http://software.broadinstitute.org/software/igv/home) was adopted to visualize ChIP-seq tracks, while ChIP-seq heat maps were generated using deepTools (47).

The Genomics of Drug Sensitivity in Cancer (GDSC; https://www.cancerrxgene.org/) database was utilized to evaluate the sensitivity of various chemotherapeutic agents. The pRRophetic package was employed to estimate the IC₅₀ of these drugs (48).

All experiments were performed in triplicate and repeated three times. All statistical analyses and subsequent visualization were implemented using R software (version 4.1.0). The data were assessed for normal distribution using the Shapiro-Wilk method and for homogeneity of variance using the Levene method. All two-group comparisons of normally distributed data were performed using unpaired Student's t test. For multigroup comparisons, one-way ANOVA with Tukey's post hoc test was used. Data that were not normally distributed or without homogeneity of variance were compared using Kruskal-Wallis with Dunn's post hoc test and Wilcoxon rank sum nonparametric tests. Data are presented as the mean ± SEM (error bars). All survival analyses were conducted using KM analysis, the log-rank test and the Cox proportional hazards model or the two-stage method (49). Pearson's test or Spearman's test was used to analyze the correlation of two variables. P<0.05 was considered to indicate a statistically significant difference.

To ascertain the potential involvement of PAXIP1 in HCC, the transcriptional levels of PAXIP1 across HCC studies were assessed using HCCDB. Analysis of 11 HCC cohorts from this database demonstrated that PAXIP1 mRNA was upregulated in HCC tissues compared with adjacent non-tumor tissues (Fig. 1A). A more granular examination of TCGA-LIHC samples using the UALCAN database further demonstrated a marked increase in PAXIP1 expression compared with that in healthy controls across all tumor grades (Fig. 1B and C). In addition, analysis of subgroups stratified by sex, age and ethnicity indicated that PAXIP1 expression was higher in patients with HCC compared with healthy controls (Fig. S1A-C). The methylation level of the PAXIP1 promoter region varied across groups based on sex, age and ethnicity (Fig. S1D-F). The mRNA expression matrix from the Cancer Cell Line Encyclopedia dataset corroborated these findings, showing that PAXIP1 was highly expressed in HCC cell lines (Fig. 1D). Furthermore, immunohistochemical analysis from the HPA database showed an absence of PAXIP1 protein in adjacent non-tumor liver tissues, while its expression was elevated in HCC tumor tissues (Fig. 1E and F). Therefore, the aforementioned results suggested that PAXIP1 expression was upregulated in HCC.

To elucidate the relationship between PAXIP1 expression and HCC prognosis, PAXIP1 genomic alterations were first analyzed in HCC using the cBioPortal website. The results showed that genomic alterations in PAXIP1 were present in 1.2% of patients (Fig. 2A). These alterations were diverse in nature (Fig. 2B). Among patients with HCC, amplification was one of the major types of PAXIP1 copy number variation (Fig. 2C). The prognostic significance of PAXIP1 expression was further analyzed in HCC. Analysis of the overall survival (OS) rate showed that patients with high PAXIP1 expression had a low survival rate (Fig. 3A). For the receiver operating characteristic (ROC) curves, the AUC value range was 0.66-0.59 for 1-, 3- and 5-year prognoses (Fig. 3B). To determine whether PAXIP1 could be used as an independent prognostic factor, univariate and multivariate Cox regression analyses were performed. Univariate Cox regression analysis indicated that PAXIP1 was a significant risk factor for OS in patients with HCC (Fig. 3C; P=0.03328).

However, further multivariate Cox regression analysis and nomogram results showed that PAXIP1 was not an independent risk factor for HCC (Fig. 3D-F). Collectively, these findings indicated that a high level of PAXIP1 may predict a low survival rate of patients with HCC.

To investigate the biological significance of PAXIP1 in HCC, PAXIP1 co-expression was examined in an HCC cohort using the LinkedOmics function module. It was shown that 5,710 genes were positively correlated with PAXIP1, while 3,317 genes were negatively correlated with PAXIP1 (false discovery rate <0.01; Fig. 4A). The heatmap further illustrates these correlations, highlighting the top 50 most significant genes, which were positively and negatively correlated with PAXIP1 expression (Fig. 4B and C). Subsequent GSEA was conducted to clarify the principal GO terms associated with PAXIP1 co-expressed genes. The analysis of GO biological process categories showed that the genes co-expressed with PAXIP1 were predominantly involved in processes such as ‘Protein alkylation’, ‘Covalent chromatin modification’, ‘Non-recombinational repair’ and ‘Maintenance of cell number’, whereas ‘Mitochondrial respiratory chain complex assembly’, ‘NADH dehydrogenase complex assembly’, ‘Translational elongation’, ‘mitochondrial gene expression’ and multiple metabolic processes were found to be downregulated (Fig. 4D; Table SI). KEGG pathway analysis indicated notable enrichment in pathways such as ‘MicroRNAs in cancer’, ‘Cell cycle’, ‘Complement and coagulation cascades’ and ‘Metabolism of xenobiotics by cytochrome P450’ (Fig. 4E; Table SII). These findings suggested the extensive impact of PAXIP1 on the global transcriptome, highlighting its potential role in HCC pathogenesis.

PAXIP1 was associated with the survival rate of patients with HCC. Therefore, a functional analysis of PAXIP1 target genes was performed. Using ChIP-seq data, 2,370 target genes were identified based on their binding scores (Table SIII). Using DAVID (26,50) and Metascape (27), functional annotation of the 2,370 PAXIP1 target genes was performed. The most enriched terms in the GO (Fig. 5A) and KEGG (Fig. 5B) analyses are shown. Notably, metabolism-related pathways or processes showed a high frequency of occurrence, including ‘Oxocarboxylic acid metabolism’, ‘Carbon metabolism’ and ‘Propanoate metabolism’, indicating a potential role of PAXIP1 in the metabolic mechanisms during HCC tumorigenesis.

The network of enriched terms elucidated the intricate interactions among the terms with considerable detail (Fig. 5C). To further understand the biological functions of the 2,370 PAXIP1 targets, the Metascape database (https://metascape.org/) was also employed to perform DisGeNET and PaGenBase enrichment analyses. The summary of enrichment analysis in DisGeNET and PaGenBase showed that the differentially expressed genes were mainly enriched in ‘Hepatomegaly’ (Fig. 5D) and ‘liver’ (Fig. 5E). These results suggested a tissue- or cell-specific role for PAXIP1 and its target genes in HCC. Therefore, PAXIP1 may be associated with HCC tumorigenesis via dysregulation of multiple pathways.

To understand the involvement of PAXIP1 in HCC and to expand the investigation to a genomic scale, the binding features of PAXIP1 in HCC cells were examined. By performing co-localization analysis of PAXIP1 ChIP-seq data (GSE104247) derived from HepG2 cells (22), 19 potential genes were identified and a multi-gene summary was performed using HCCDB. These genes were separated into two groups (upregulated and downregulated) based on their expression levels in HCC (Fig. 6A). Analysis using the GEPIA database revealed a marked increase in MYBL2 expression and a notable decrease in FOXO1 expression in HCC (Fig. 6B and C). Notably, FOXO1 expression was not significantly changed when combining TCGA data with Genotype-Tissue Expression data (Fig. 6C).

The genomic distribution of PAXIP1, MYBL2 (GSE32465) and FOXO1 (GSE104247) was then analyzed using data from previous studies (21,22). To visualize the overlap of PAXIP1, MYBL2 and FOXO1 binding sites, chromosomal folding was represented using a Hilbert curve, which preserves the spatial proximity of linearly adjacent regions (51). The distributions of PAXIP1, MYBL2 and FOXO1 exhibited similar patterns at the whole-genome scale (Fig. 6D). It was further revealed that PAXIP1 shares 14,324 peaks with MYBL2 and 6,991 peaks with FOXO1 (Fig. 6E and F). To further analyze the relationship between PAXIP1, MYBL2 and FOXO1, the signals of PAXIP1, MYBL2 and FOXO1 were plotted in descending order in the PAXIP1/MYBL2- and PAXIP1/FOXO1-cobound regions. PAXIP1 occupancy showed a similar pattern with MYBL2 and FOXO1 in their co-bound regions (Fig. 6G and H). Survival analysis revealed that in patients with HCC, elevated MYBL2 expression was associated with poorer OS, while reduced FOXO1 expression was also associated with worse OS outcomes. However, only the difference in OS related to MYBL2 expression was statistically significant, whereas the difference associated with FOXO1 expression was not significant (Fig. 6I and J). Taken together, these findings demonstrated that PAXIP1 functioned as a cofactor for MYBL2 or FOXO1 in HCC.

TFs are responsible for regulating gene expression in HCC development (52). To assess whether PAXIP1 was modulated by TFs, the upstream TFs of PAXIP1 were first predicted and the top 10 predicted TFs were identified using the hTF target database (Fig. 7A). Subsequently, survival analysis for all 10 TFs was performed (Figs. 7B, C and S2). Among these TFs, CTCF and NRF1 have been demonstrated to be prognostic markers in liver cancer (53,54), according to the HPA database. No statistically significant differences were found; however, there was a tendency that elevated CTCF and NRF1 expression levels were associated with poorer patient outcomes in HCC (Fig. 7B and C). To confirm the role of NRF1 and CTCF in the regulation of PAXIP1 in HCC, NRF1 and CTCF were knocked down with targeted siRNA. Decreased expression of NRF1 resulted in decreased PAXIP1 expression in HuH-7 and PLC-PRF-5 cells; however, CTCF knockdown did not affect PAXIP1 expression (Fig. 7D). The aforementioned findings suggested that NRF1 may regulate PAXIP1 expression during the development of liver cancer.

Since potential cofactors and regulators of PAXIP1 were identified, the prognostic significance of these genes in HCC was assessed. Patients with HCC were redivided into the high-risk (n=185) and low-risk groups (n=185) according to the median risk score (Fig. 7E-G). Those in the high-risk group exhibited worse outcomes than their low-risk counterparts (Fig. 7H). For the ROC curve, the AUC range was 0.735-0.649 for the 1-, 3- and 5-year prognoses (Fig. 7I). The results indicated that the combined expression of four genes, PAXIP1, MYBL2, FOXO1 and NRF1, could serve as an effective prognostic marker for HCC.

To elucidate the molecular mechanisms underlying the role of PAXIP1 in HCC, patients were stratified into two subgroups based on PAXIP1 expression levels: i) The PAXIP1-high group (PAXIP1-High; n=186); and ii) the PAXIP1-low group (PAXIP1-Low; n=185), determined by the median expression value. Differential expression analysis between these two subgroups was conducted (Fig. 8A and B). The subsequent GO and KEGG results revealed that the upregulated genes in the PAXIP1-High group were predominantly associated with ‘cell cycle’ pathways (Fig. 8C), as well as with processes related to ‘organelle fission’, ‘nuclear division’ and ‘chromosome segregation’ (Fig. 8D). Conversely, the downregulated genes in the PAXIP1-High group were involved in ‘PPAR signaling pathway’, ‘metabolism of xenobiotics by cytochrome P450’ and ‘complement and coagulation cascades’ (Fig. 8E), and in ‘lipid localization’, ‘glycerolipid metabolic process’ and ‘acute inflammatory response’ (Fig. 8F). These findings indicated that the differential expression of PAXIP1 target genes may result in the dysregulation of cell division, immune responses and metabolic processes.

The chemotherapeutic response of each sample was assessed using the GDSC database. The 50% maximal inhibitory concentrations for the samples were estimated via ridge regression, and all the involved parameters were set to their default values. To eliminate batch effects, combat normalization was applied, and duplicate gene expression values were summarized by taking their mean. The results showed that the sensitivities to the chemotherapeutic drugs axitinib, lenalidomide, pazopanib, sorafenib and XL-184 (Fig. 9A-E; P<0.001) were significantly negatively correlated with PAXIP1 expression in HCC. These analyses indicated that upregulation of PAXIP1 expression may increase sensitivity to the aforementioned chemotherapy drugs in HCC treatment.

Given the possible oncogenic function of PAXIP1 in HCC, it is imperative to explore the relationship between PAXIP1 and immune events in HCC. Utilizing the TIMER database, the correlation between PAXIP1 expression and immune cell infiltration was evaluated. The analysis revealed a positive association between PAXIP1 expression and the presence of CD4⁺ T cells, neutrophils, macrophages, B cells and myeloid dendritic cells within HCC tissues (Fig. 10A and B). Immune checkpoints such as programmed cell death protein 1 (PD1/PDCD1)/programmed death-ligand 1 (PD-L1/CD274) and cytotoxic T-lymphocyte associated protein 4 (CTLA4) are crucial for the modulation of immune responses and tumor immune evasion (55). The current study demonstrated a significant positive correlation between PAXIP1 expression and the levels of PDCD1, CD274 and CTLA4 (Fig. S3A-C; P<0.05). Further analysis using the GEPIA database corroborated these findings, revealing a weak positive correlation between PAXIP1 and PDCD1, CD274 and CTLA4 expression in HCC (Fig. S3D-F). Therefore, PAXIP1 may contribute to carcinogenesis in HCC through mechanisms involving immune cell infiltration and tumor immune escape.