Air pollution-related genes were obtained from Pu et al (11), who investigated seven air pollutants in HCC. To expand the pollutant list, three additional pollutants, namely 1,3-butadiene, sulfur trioxide and vinyl chloride, were included in the present study. These three pollutants were searched in PubChem (https://pubchem.ncbi.nlm.nih.gov/) using their exact chemical names as search terms, namely ‘1,3-butadiene’, ‘sulfur trioxide’ and vinyl chloride. The corresponding pollutant-related genes were then retrieved and integrated with the genes associated with the seven pollutants reported by Pu et al (11), resulting in a ten-pollutant air pollution-related gene set. Pu et al (11) mainly focused on air pollutant-related immune genes in HCC based on seven pollutants. By contrast, the present study was based on a broader set of genes associated with ten air pollutants and was not limited to air pollution-related immune genes. This approach may provide a more comprehensive view of the potential molecular association between air pollution exposure and HCC. The ADMETLAB 3.0 platform (https://admetlab3.scbdd.com) and the ProTox3 database (https://tox.charite.de/protox3) were used to assess carcinogenicity with their default prediction settings. For each AP, the canonical chemical structure was submitted directly to both platforms without additional parameter modification. The carcinogenicity-related probability output generated by each platform was recorded and a value >0.5 was considered indicative of carcinogenic potential. Pollutants were defined as carcinogenic when either ADMETLAB 3.0 or ProTox3 returned a carcinogenicity probability >0.5.

The human target genes for the 10 APs were obtained from the SuperPred database (https://prediction.charite.de/), Swiss Target Prediction database (http://www.swisstargetprediction.ch/), STITCH database (http://stitch.embl.de/) and Similarity Ensemble Approach database (https://sea.bkslab.org/). The genes from all databases were merged and duplicates were removed.

RNA-sequencing data for HCC were derived from The Cancer Genome Atlas (TCGA; 50 normal samples; 374 HCC samples; http://www.cancer.gov/ccg/research/genome-sequencing/tcga), the International Cancer Genome Consortium (ICGC; 232 HCC samples; http://dcc.icgc.org/) and GSE14520 (221 HCC samples; http://www.ncbi.nlm.nih.gov/geo/). Within each dataset, expression matrices were preprocessed and normalized according to the original data format, and cross-dataset batch effects were corrected using the ComBat function in the ‘sva’ R package (version 4.3.2; http://www.r-project.org/) before downstream integrative analyses. In the weighted gene co-expression network analysis (WGCNA), genes with a standard deviation of expression level <0.5 were removed. The optimal soft threshold was determined to be 18. Dynamic tree cutting and hierarchical clustering were then used to identify modules. The cutting height was 0.3 and the minimum module size was 100 genes.

Differentially expressed genes (DEGs) between HCC and normal liver samples with |log₂ fold change|>1 and false discovery rate (FDR) q<0.05 were subsequently identified.

HCC-related genes were retrieved from the GeneCards database (https://www.genecards.org/) with a threshold of >5 score.

The ‘Org.Hs.eg.db’, ‘enrichplot’ and ‘clusterProfiler’ packages in R software were used to perform Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) enrichment analyses of AP-related and HCC-related genes (AP-HCC-Gs). GeneMANIA (https://genemania.org/) was used to evaluate the biological functions.

TargetScan Human (https://www.targetscan.org), miRDB (https://mirdb.org/), miRanda (https://www.microrna.org/) and miRWalk (https://mirwalk.umm.uni-heidelberg.de/) were utilized to predict the binding of miRNAs and mRNAs. Cytoscape (v3.8.2; http://cytoscape.org/) was used for visualization.

Batch effects between the TCGA cohort (training set) and ICGC cohort (external validation set) were corrected using the ‘sva’ R package. After univariate Cox regression analysis was performed in the training set to identify AP-HCC-Gs associated with overall survival, 10 machine-learning algorithms, including least absolute shrinkage and selection operator (LASSO), Ridge, StepCox, CoxBoost, random survival forest, generalized boosted regression modeling, survival support vector machine, supervised principal components (SuperPC), elastic net and partial least squares regression for Cox, were used to develop 101 prognostic models with 10-fold cross-validation. These models included both single-algorithm models and sequential combinations, in which the first algorithm was applied for variable selection and the second algorithm was used for model fitting. For single-algorithm models, all candidate variables were entered directly into model construction. Only genes shared by the TCGA and ICGC datasets were retained, and the expression matrices were standardized before analysis. Model combinations yielding ≤3 variables at the variable-selection step or retaining ≤3 variables in the final fitted model were excluded (seed, 1,234). Model performance was evaluated using the concordance index (C-index) in both the training and external validation cohorts, and models were ranked according to the average C-index across cohorts. The top-ranked model was retained for downstream analyses.

Principal component analysis (PCA), the area under the curve (AUC), Kaplan-Meier (K-M) curves, risk curves and univariate and multivariate Cox analyses were used to further validate the robustness.

Protein structures were obtained from the Alphafold (https://alphafold.ebi.ac.uk/) and Protein Data Bank (https://www.rcsb.org/) databases. CB-Dock2 (https://cadd.labshare.cn/) was used for molecular docking analysis.

For diagnostic marker screening, LASSO, support vector machine recursive feature elimination (SVM-RFE) and random forest algorithms were applied in the discovery dataset using 10-fold cross-validation or the corresponding internal optimization procedure as implemented in each algorithm. Genes consistently retained across the three methods were defined as candidate diagnostic markers and were subsequently evaluated by receiver operating characteristic (ROC) analysis (‘timeROC’ R package). Because these markers were derived from the discovery pipeline, their diagnostic performance should be interpreted as preliminary until further independent validation is available. Subsequently, AP-HCC-DMs were further intersected with AP-HCC-PS genes to obtain hub genes. The AUC was used to evaluate the diagnostic value.

Using the ‘rms’ R package, prognostic and diagnostic nomograms were constructed. The C-index was used to evaluate the consistency of the predicted values and the observed values.

Single-cell analysis was used to localize the hub gene within immune cell subpopulations and, through cell-cell communication analysis, to predict the communication frequency between the hub gene-enriched immune subset and other immune cells, thereby indicating potential AP-gene-immune cell pathways through which AP-induced gene expression may promote an immunosuppressive tumor microenvironment in HCC.

THLE-2 and Huh-7 cells were obtained from the Institute of Cell Research, Chinese Academy of Sciences. Huh-7 cells were cultured in DMEM (Gibco; Thermo Fisher Scientific, Inc.), whereas THLE-2 cells were cultured in bronchial epithelial growth medium (BEGM; Lonza, Switzerland) with the recommended growth supplements added. Both media were supplemented with 10% fetal bovine serum (Gibco; Thermo Fisher Scientific, Inc.) and 1% penicillin-streptomycin (Gibco; Thermo Fisher Scientific, Inc.). All cells were maintained at 37°C in a humidified incubator with 5% CO₂. Total RNA was isolated using TRIzol reagent (cat. no. 15596-026; Invitrogen; Thermo Fisher Scientific, Inc.) according to the manufacturer's instructions. cDNA was synthesized using a reverse transcription kit (K1622; Fermentas; Thermo Fisher Scientific, Inc.) according to the manufacturer's instructions. RT-qPCR was performed using a SYBR Green PCR kit (K0223; Thermo Fisher Scientific, Inc.) on an ABI 7300 Real-Time PCR System (Applied Biosystems; Thermo Fisher Scientific, Inc.). The PCR reaction mixture (25 µl) contained 10 µl SYBR Green Mix, 1 µl forward primer, 1 µl reverse primer, 2 µl cDNA template and 11 µl nuclease-free water. The thermocycling conditions were as follows: Initial denaturation at 95°C for 10 min, followed by 40 cycles of denaturation at 95°C for 15 sec and annealing/extension at 55°C for 45 sec. Melting curve analysis was subsequently performed at 95°C for 15 sec, 60°C for 1 min, 95°C for 15 sec and 60°C for 15 sec. Relative mRNA expression levels were quantified using the 2^−ΔΔCq method with GAPDH as the internal control (12). The primer sequences used for RT-qPCR were as follows: Proteasome 20S subunit β5 (PSMB5) forward, 5′-GCTACAGCGGGTGCTTACAT-3′ and reverse, 5′-TTCCCAGAAGCTGCAATCCG-3′; GAPDH forward, 5′-CACCCACTCCTCCACCTTTG-3′ and reverse, 5′-CCACCACCCTGTTGCTGTAG-3′.

Cells were transfected with PSMB5-targeting siRNA or the corresponding negative control (NC) siRNA using Lipofectamine 2000 (cat. no. 11668-019; Invitrogen; Thermo Fisher Scientific, Inc.) according to the manufacturer's instructions. At 24 h post-transfection, cells were harvested for total RNA extraction, and RT-qPCR was performed to determine the knockdown efficiency of PSMB5 at the transcript level. The siRNA sequences used for transfection were as follows: siNC, 5′-UUCUCCGAACGUGUCACGUTT-3′; and siPSMB5, 5′-CGAAAUGCUUCAUGGAACA-3′.

For the wound-healing assay, cells in the control and experimental groups were seeded in 35-mm dishes at a density of 8×10⁵ cells/dish and cultured until reaching a confluent monolayer. A linear scratch was created perpendicular to the reference line marked on the bottom of the dish using a 10-µl pipette tip. Detached cells were removed by PBS washing and cells were then incubated in serum-free medium. Images were captured at 0, 24 and 48 h at the same marked position under a microscope to evaluate cell migration. Wound closure was quantified using ImageJ software (version 1.53; National Institutes of Health).

Single-cell analysis was used to localize the hub gene within immune cell subpopulations and, through cell-cell communication analysis, to predict the communication frequency between the hub gene-enriched immune subset and other immune cells, thereby indicating potential AP-gene-immune cell pathways through which AP-induced gene expression may promote an immunosuppressive tumor microenvironment in HCC.

In the present study, scRNA-seq experiments were not performed on newly collected samples. Instead, publicly available scRNA-seq data from 12 HCC samples, consisting of treatment-naïve tumor tissues from patients with pathologically confirmed primary HCC, were obtained from the China National Gene Bank Nucleotide Sequence Archive (CNP0000650; http://db.cngb.org/cnsa).

After excluding genes expressed in <3 cells and retaining cells with an nFeature of 6,000-10,000, nCount of 4,000-15,000 and percentage of mitochondrial gene counts (percent.mt) <5, 19,516 genes and 8,680 cells were obtained. For scRNA-seq data, LogNormalize was applied for data normalization, the ‘Harmony’ package was applied for correction of batch effects and PCA was applied for dimensionality reduction. The ‘FindNeighbors’ and ‘FindClusters’ functions were used for cell clustering and resolution determination. Uniform manifold approximation and projection were used for visualization. Cell annotation was performed using marker genes from previous studies (13–15). Cell-cell communication was analyzed using the CellChat R package.

The infiltration abundance of 22 immune cell types was calculated using CIBERSORT (https://cibersortx.stanford.edu/). Pearson analysis was further utilized to evaluate the correlation between the hub genes and immune cell types in HCC.

Statistical analyses were performed using R software (version 4.3.2; http://www.r-project.org/). Unless otherwise specified, two-group comparisons of continuous variables were performed using Student's t-test or the Wilcoxon rank-sum test, depending on data distribution. Categorical variables were compared using the χ² test. Survival analyses were conducted using K-M curves (‘survival’ R package) and Cox regression models. P-values were two-sided, and P<0.05 was considered to indicate a statistically significant difference. Where appropriate, multiple testing was adjusted for using the FDR method.

Fig. S1 shows the research flowchart of the present study. All 10 APs exhibited carcinogenicity based on both toxicity prediction platforms (Table SI; Fig. 1A). Target genes for the 10 APs were obtained: 133 for 1,3-butadiene, 82 for benzene, 70 for CO, 92 for NO, 87 for NO₂, 88 for O₃, 92 for SO₂, 83 for SO₃, 102 for toluene and 173 for vinyl chloride. After removing duplicates, a total of 354 AP-related target genes were identified (Table SII).

WGCNA identified six co-expression modules (Fig. 1B). Among them, the blue module (3,344 genes) had the highest correlation with HCC (r=0.7; P=3×10⁻¹²; Fig. 1D and E). Module membership and gene significance showed a significant correlation (r=0.64; P<1×10⁻²⁰⁰; Fig. 1C).

Differential gene expression analysis of mRNAs in TCGA-liver HCC dataset identified 6,479 DEGs (Fig. 1F).

A total of 4,423 HCC-related genes were obtained from GeneCards. After taking intersection, 655 HCC-related genes were obtained (Fig. 1G; Table SIII).

The intersection of 354 AP-related genes and 655 HCC-related genes was determined, revealing 43 AP-HCC-Gs (Fig. 1H; Table SIII).

GO analysis revealed that the AP-HCC-Gs were primarily involved in ‘mitotic cell cycle phase transition’ (GO:0044772), ‘regulation of protein serine/threonine kinase activity’ (GO:0071900), ‘protein autophosphorylation’ (GO:0046777), ‘positive regulation of transferase activity’ (GO:0051347), ‘nuclear division’ (GO:0000280) and ‘G₂/M transition of mitotic cell cycle’ (GO:0000086) (Fig. S2A; Table SIV).

KEGG pathway analysis suggested that the top three most enriched pathways were ‘microRNAs in cancer’, ‘cell cycle’ and ‘MAPK signaling pathway’ (Fig. S2B; Table SV). Furthermore, a number of immune and metabolic pathways were enriched, including the ‘Fc epsilon RI signaling pathway’, ‘neutrophil extracellular trap formation’, ‘growth hormone synthesis, secretion and action’, ‘chemical carcinogenesis-reactive oxygen species’ and ‘thyroid hormone signaling pathway’.

In addition, due to the highly enriched ‘miRNAs in cancer’ pathway, a miRNA-mRNA network was constructed, providing a reference for relevant studies to understand core regulatory axes driven by miRNAs (Fig. S3; Table SVI).

The univariate Cox analysis identified 36 AP-HCC-Gs significantly associated with HCC prognosis (Fig. 2G; Table SVII). Among the 101 models, the StepCox (both) + SuperPC model and the StepCox (backward) + SuperPC model achieved the highest average C-index (Fig. 2C). As both models showed equivalent predictive performance and retained the same seven genes, namely FKBP prolyl isomerase 1A (FKBP1A), NAD(P)H quinone dehydrogenase 1 (NQO1), aurora kinase B (AURKB), hydroxymethylbilane synthase (HMBS), enhancer of zeste 2 polycomb repressive complex 2 subunit (EZH2), PSMB5 and X-ray repair cross complementing 1 (XRCC1), the StepCox (both) + SuperPC model was retained as the final model because it preserved variables selected in both directions and was therefore considered to be more comprehensive. In this signature, lower expression of these genes was associated with improved prognosis (Fig. S4A-F).

PCA and K-M curve analysis indicated that the AP-HCC-PS could distinguish between different prognoses and had good predictive performance (Fig. 2E and F). The risk curve and heat map indicated that the higher the AP-HCC-PS gene expression level, the worse the prognosis (Fig. 2D and H). The univariate Cox analysis [hazard ratio (HR), 2.224; 95% CI, 1.701-2.908; P<0.001; Fig. S4G] and multivariate Cox analysis (HR, 2.076; 95% CI, 1.548-2.783; P<0.001; Fig. S4H) showed that the AP-HCC-PS exhibited predictive performance independent of clinical features. Decision curve analysis (DCA) revealed that the AP-HCC-PS yielded a higher net benefit than other clinical features across a range of thresholds (Fig. S4I). The ROC curve suggested that the AP-HCC-PS had better predictive performance (AUC, 0.768; Fig. 2B) than clinical features, and had good predictive performance in predicting 1-year (AUC, 0.768), 3-year (AUC, 0.679) and 5-year (AUC, 0.657) survival (Fig. 2A). In the training set, the model sensitivity decreased with longer prediction horizons, which may be partly attributable to the limited sample size and the introduction of additional noise over time. The ICGC cohort and the additional validation cohort GSE14520 further validated the generalizability and robustness of the model (Figs. S5A-G and S6A-G).

LASSO (Fig. 3A and B), Random-forest algorithm (Fig. 3C and D) and SVM-RFE (Fig. 3E and F) methods identified 16, 7 and 15 genes, respectively. After taking the intersection, PKMYT1 and PSMB5 were identified as AP-HCC-DMs (Fig. 3G; Table SVIII). The AUCs indicated that PKMYT1 (AUC, 0.971) and PSMB5 (AUC, 0.965) showed high diagnostic value (Fig. 3H).

The prognostic nomogram (C-index, 0.722; 95% CI, 0.674-0.769; Fig. S7B and D) and the diagnostic nomogram (C-index, 0.984; 95% CI, 0.972-0.995; Fig. S7A and C) showed good accuracy in predicting HCC incidence and prognosis. DCA revealed that both prognostic and diagnostic nomograms yielded a superior net benefit across a range of thresholds, indicating favorable clinical accuracy (Fig. S7E and F). The DCA demonstrated that both the AP-HCC-PS risk score and the integrated clinical nomogram achieved a higher net benefit than individual clinical features over a wide spectrum of threshold probabilities, with the nomogram outperforming the risk score alone (Fig. S7G).

The molecular docking analysis showed that the predicted docking scores between the AP-HCC-PS/AP-HCC-DM targets and the 10 APs were generally low, with average binding energies ranging from −2 to −3 kcal/mol. These findings provided only exploratory support for potential target-pollutant interactions (Fig. 4A; Table SIX). Among these exploratory results, PSMB5 showed relatively lower docking scores with toluene (−4.2; Fig. 4B), benzene (−3.8; Fig. 4C) and 1,3-butadiene (−2.2; Fig. 4D), suggesting that it may warrant further investigation. Furthermore, molecular docking of known inhibitors for the seven targets was performed to aid the interpretation of the AP-gene relationships. For all seven targets, the inhibitor-target binding energies were ≤-6 kcal/mol, whereas the pollutant-target docking scores were substantially weaker (Table SX). This result confirmed the robustness of the docking procedure and suggested that AP-gene interactions differ from drug-target inhibition: While drugs achieve rapid inhibition through strong binding, APs may modulate target traits through long-term and cumulative effects on gene expression.

After intersecting the AP-HCC-PS with the two AP-HCC-DMs (Fig. S8A), PSMB5 was identified as the hub gene, and its high expression was found to be associated with poor prognosis in HCC (Fig. S8B).

RT-qPCR results showed that PSMB5 was highly expressed in HCC cells (Fig. 5A; Table SXI). PSMB5 expression was significantly decreased following siRNA-mediated knockdown (Fig. 5B; Table SXII). The wound-healing assay showed no significant difference among the three groups at 0 h (Fig. 5C). At 24 h, the PSMB5 siRNA group differed significantly from both the blank group (P=0.0082) and the NC group (P=0.0003), while no significant difference was found between the blank and NC groups. At 48 h, the wound-healing rate in the PSMB5 siRNA group was significantly lower than that in the blank group (P<0.0001) and the NC group (P=0.0002) (Fig. 5C-I; Table SXIII).

PSMB5 mainly interacts with its family members and is involved in biological functions such as immune responses, oxidative stress and classical cancer pathways (Fig. S8C). scRNA-seq analysis was strategically employed to resolve the cellular context of PSMB5. scRNA-seq results showed that nCount RNA was negatively correlated with percent.mt (r=−0.72) and positively correlated with nFeature RNA (r=0.54) (Fig. 6A). Based on the elbow plot, the top 10 principal components were selected for subsequent analysis (Fig. S9C). Six clusters were identified (resolution, 0.2; Figs. 6B, and S9A and B), which were subsequently annotated into five cell types (Fig. 6C and D): B cells, stromal cells, regulatory T cells (Tregs), hepatocytes and cytotoxic T cells (CTLs). Among the immune cell populations, PSMB5 was mainly expressed in Tregs (Figs. 6E and S9E). CIBERSORT revealed a positive correlation between PSMB5 and Tregs (Fig. 6F; Table SXIV). Tregs (sender) and CTLs (receiver) exhibited the highest interaction number, which was mainly mediated by the MIF-(CD74+CXCR4) axis (Figs. 6G and H, and S9D). These findings suggested that PSMB5 may be related to Treg activity and may participate in Treg-CTL communication via the MIF-(CD74+CXCR4) axis, potentially contributing to an immunosuppressive tumor microenvironment in HCC.