The Cancer Genome Atlas Liver Hepatocellular Carcinoma (TCGA-LIHC) and Gene set enrichment 14520 (GSE14520) datasets were obtained from TCGA (https://tcga-data.nci.nih.gov/tcga/) and Gene Expression Omnibus (GEO, https://www.ncbi.nlm.nih.gov/geo/) datasets. APPGs were obtained through ImmPort (https://www.immport.org/home).

The R package limma (version 3.6.1) was utilized to process the processing of gene expression matrix data. With the addition of APP.Score, immune scores and immune cell infiltration were assessed by the ESTIMATE algorithm and TIMER database, respectively [24, 25]. In addition, we analyzed the expression of immune cell marker genes with different APP.Scores of LIHC samples using the TIMER database. The Gene Set Variation Analysis (GSVA) package was used for analyzing the Kyoto Encyclopedia of Genes and Genomes (KEGG, https://www.genome.jp/kegg/) pathway via GSVA [26]. Additionally, correlation analyses were conducted among prognostic scoring, gene scoring, and functional enrichment pathways. Subsequently, the ClusterProfiler R package was utilized for performing GSE analysis (GSEA) based on the APP.Score.

The mutated genes in the high and low APP.Score groups were analyzed using the “maftools” R package [27]. The copy number variation (CNV) landscape and the assessment of copy number gains or losses at amplified or deleted peaks, was evaluated using GISTIC 2.0 analysis (https://gatk.broadinstitute.org) [28].

The somatic mutations and CNV profiles were collected from the TCGA-LIHC dataset. The drug sensitivity of the high APP.Score and low APP.Score groups was analyzed using the oncoPredict package (version 0.2) [29].

Human normal hepatic stellate cells (LX-2, AW-CNH008; Abiowell Biotechnology Co., Ltd., Changsha, China) and HCC lines, including HepG2 (AW-CCH024; Abiowell) and Huh-7 (AW-CCH089; Abiowell), were cultured in a humidified incubator at 37 °C and 5% CO${}_2$. Phospholipase A2, group 7 (PLA2G7) silencing (si-PLA2G7) and its negative control (si-NC) were transfected into HepG2 and Huh-7 cells using Lipofectamine 2000 (2028090; Invitrogen, Carlsbad, CA, USA) for a duration of 48 h. All cell lines were validated by short tandem repeat profiling and tested negative for mycoplasma.

The six selected genes were validated by quantitative real-time PCR (qRT-PCR) using a fluorescence qRT-PCR instrument (PIKOREAL96, Thermo, Waltham, MA, USA). Primer sequences were designed utilizing primer 5 software, with $\beta$-actin serving as an internal mRNA reference and listed in Supplementary Table 1. Relative gene expression was analyzed employing the 2${}^{-\mathrm{\Delta }\mathrm{\Delta }\text{Ct}}$ method.

The cells were lysed, centrifuged, and concentrated. Then the proteins were separated by electrophoresis and electrotransferred to a nitrocellulose membrane. The membrane was incubated overnight at 4 ℃ with the following primary antibodies: killer cell lectin like receptor B1 (KLRB1), PLA2G7, musculin (MSC), heat shock protein family A (HSPA6), secreted phosphoprotein 1 (SPP1), neutrophil cytosolic factor 2 (NCF2), programmed death-ligand 1 (PD-L1), phosphoinositide 3-kinase (PI3K), AKT, phosphorylated signal transducer and activator of transcription 1 (p-STAT1), STAT1, phosphorylated nuclear factor kappa B (p-NF-$\kappa$B), and NF-$\kappa$B with $\beta$-actin as the internal reference. The membrane was incubated with horseradish peroxidase-conjugated goat anti-mouse/rabbit immunoglobulin G (IgG) secondary antibody at 37 ℃ with shaking for 90 min, followed by the development and visualization of the proteins. After washing with phosphate buffer, the membrane was immersed in ECL reaction solution and incubated at room temperature for 1 minute. The protein was exposed and imaged using a chemiluminescence imaging system (ChemiScope 6100, CLINX, Shanghai, China). Detailed antibody information was provided in Supplementary Table 2.

The Cell Counting Kit-8 (CCK-8) assay was used to measure cell proliferation. Cells were separately cultured for different time points (12, 24, 48 h) and then added to medium containing 10% CCK-8 (NU679, DOJINDO Laboratories, Kumamoto, Japan). The cells were incubated at 37 ℃ for an additional 4 h, and their optical density values were analyzed at 450 nm.

The cells were digested into single cells, and 200 cells/per well were seeded and maintained for 10 days. After the cells were fixed in 4% paraformaldehyde and stained with crystalline violet, photographs were taken for observation.

As previously described [30], a monolayer of cells was scratched. The cells were rinsed three times with phosphate-buffered saline to remove the scratched cells, and serum-free medium was added. Cell healing was observed under an optical microscope (DSZ2000X, Beijing Zhongyantaihe Medical Instrument Co., Ltd., Beijing, China) at 0 and 48 h.

Matrigel (354262, Becton Dickinson and Co., Franklin, NJ, USA) that had been prediluted was added to the upper chamber, followed by the addition of a cell suspension of 2 $\times$ 10${}^4$ to the same chamber. The lower chamber was filled with medium containing 20% fetal bovine serum. After incubating and wiping down the upper chamber, the cells were fixed, stained, and observed under an optical microscope. The results were analyzed using ImageJ software (version 1.49, National Institutes of Health, Bethesda, MD, USA).

Cell apoptosis was detected using an apoptosis detection kit (KGA1030, KeyGEN BioTECH, Jiangsu, China) through the instructions. After incubation, the results were immediately evaluated using FCM (A00-1-1102, Beckman, Brea, CA, USA).

Cells were fixed, permeabilized, and closed. Then cells were incubated overnight at 4 °C with primary antibody against PD-L1 (1:50, 28076-1-AP, Proteintech, Chicago, IL, USA), followed by incubation for 90 min at 37 °C with anti-rabbit IgG secondary antibody (1:200, SA00013-4, Proteintech). Finally, the cell nuclei were stained and visualized under a fluorescence microscope (BA210T, Motic Microscopes, Chautauqua County, NY, USA).

We conducted our statistical analysis and data visualization by utilizing R version 3.6.1 along with the ggplot2 package. Additionally, for variables that were not normally distributed, we employed both Wilcoxon and Kruskal-Wallis tests as part of our analytical approach. To calculate correlation coefficients, we used Pearson correlation. To generate and display survival curves for our subgroups, we used the Kaplan-Meier method. All of our tests were two-sided. To determine statistical significance, we used GraphPad Prism software package version 8.0 to perform analysis of variance (ANOVA) or unpaired Student’s t tests. Statistical significance was determined by a p-value less than 0.05.

Based on the 23 APP genes, we obtained two genotypes through cluster analysis: APP cluster 1 and APP cluster 2 (Fig. 1A). The principal component analysis plot showed that APP cluster 1 and APP cluster 2 had good dispersion (Fig. 1B). We conducted survival analysis on the two clusters. When the survival rate was 0.5, the survival time of cluster 2 was longer than that of cluster 1 (Fig. 1C). Subsequently, the differential expression of the 23 APP genes in cluster 1 and cluster 2 was analyzed (Fig. 1D,E). The results showed that compared with cluster 1, 17 genes including adaptor-related protein complex 3 subunit beta 1, Fc Fragment of IgE Receptor Ig (FCER1G), HSP90AA1, HSP90AB1, HSPA4, HSPA6, HSPA8, nuclear transcription factor Y subunit gamma, proteasome 26S subunit ATPase 4, proteasome 26S subunit, non-ATPase 11 (PSMD11), PSMD14, PSMD2, PSMD7, proteasome activator complex subunit 3, retinoic acid early transcript 1G (RAET1G), RAET1L, and transient receptor potential cation channel, subfamily C, member 4 associated protein were downregulated in cluster 2. The above results reveal that the APP-related genotyping analysis is reasonable.

Fig. 1.

Recognition of antigen-presenting prognostic (APP)-related genotyping. (A) Cluster analysis. (B) Principal component analysis (PCA) analysis. (C) Survival analysis. (D) Heatmap. (E) Bar plot. *p 0.05; **p 0.01; ***p 0.001; ****p 0.0001; ns means no statistical significance; vs. the APP. cluster 1.

The minimum soft threshold value for building a scale-free network was 5 (scale-free fit index = 0.9). Therefore, 5 was chosen as the optimal soft threshold value for subsequent analyses. Another Figure showed the network connectivity under soft thresholding power (Fig. 2A). Subsequently, a gene clustering tree was constructed (Fig. 2B). The correlation and significance of different gene modules with cluster 1 and cluster 2 are shown in Fig. 2C. The green gene module exhibited the highest correlation with both cluster 1 and cluster 2. Further analyses revealed a positive correlation between module membership and gene significance (Fig. 2D), suggesting that genes highly related to cluster 1 and cluster 2 were also crucial in the green gene module.

Fig. 2.

Weighted gene co-expression network analysis (WGCNA). (A) Network topology analysis. (B) Module clustering tree. (C) Phenotypic and module correlation plot. (D) Correlation analysis between genes and phenotypes in modules.

Subsequently, 39 genes were obtained by performing univariate analysis on the genes in the green module (Fig. 3A). Nine significant genes, including KLRB1, NCF2, PLA2G7, MSC, HSPA6, NPL, FCER1G, SPP1, and FABP5P7, were obtained through further RF analysis (Fig. 3B). Then a LASSO model was constructed. Seven genes, including KLRB1, NCF2, PLA2G7, MSC, HSPA6, SPP1, and FABP5P7, were selected. Finally, the APP.Score was obtained by calculating the formula (APP.Score = –2.2047 $\times$ KLRB1 + 0.0197 $\times$ NCF2 + 1.2247 $\times$ PLA2G7 + 0.7626 $\times$ MSC + 0.78 $\times$ HSPA6 + 0.7743 $\times$ SPP1 + 0.8317 $\times$ FABP5P7) (Fig. 3C). Upon further investigation, it was discovered that patients with a high APP.Score experienced a greater number of deaths than those with a low APP.Score (Fig. 3D). The area under the curve (AUC) values were 0.734, 0.747, and 0.679 for survival periods of 1, 2, and 3 years, respectively (Fig. 3E–G). These results suggest that APP.Score is a risk factor for HCC.

Fig. 3.

APP.Score was a risk factor for Liver Hepatocellular Carcinoma (LIHC). (A) Univariate analysis. (B) Random forest (RF). (C) Lasso analysis. (D) Risk evaluation of APP.Score. (E) The Cancer Genome Atlas (TCGA) survival analysis. (F) Receiver operator characteristic (ROC) curve. (G) GSE14520 survival analysis.

Supplementary Fig. 2A shows that with the increase in APP.Score, the expression of KLRB1, NCF2, PLA2G7, MSC, HSPA6, SPP1, and FABP5P7 was significantly increased. With the increase in APP.Score, stromal score, immune score, estimate score, B cells, CD4${}^{+}$ T cells, CD8${}^{+}$ T cells, neutrophil, MAC, and DC enrichment were significantly increased, whereas tumor purity enrichment was significantly decreased (Supplementary Fig. 2B). With the increase in APP.Score, the related genes were significantly enriched including AP-related (e.g., HLA-DRB5, HLA-DQB1, and HLA-DPB1), cell adhesion-related (integrin subunit beta 2 and intercellular adhesion molecule 1), co-inhibitor-related (e.g., cluster of differentiation 274 [CD274], butyrophilin subfamily 3 member A2 [BTN3A2], and BTN3A1), co-stimulator-related (CD80 and CD28), ligand-related (e.g., tumor necrosis factor [TNF], TNFSF4, CD70, and interferon alpha 1), receptor-related (e.g., T cell immunoreceptor with Ig and immunoreceptor tyrosine-based inhibitory motif domains, interleukin 2 receptor subunit alpha, and TNF receptor superfamily member 14), and others (e.g., indoleamine 2, 3-dioxygenase 1, granzyme A, perforin 1, high mobility group box 1 protein, and ectonucleoside triphosphate diphosphohydrolase 1), whereas ARG1 enrichment was significantly decreased (Supplementary Fig. 2C).

We further analyzed the correlation of mutated genes in the high and low APP.Score groups. In the high APP.Score group, tumor protein p53 (TP53) and retinoblastoma 1 (RB1) might have co-occurring mutations; $\beta$-catenin (CTNNB1) might have co-occurring mutations with WD repeat domain 87 (WDR87), obscurin (OBSCN), and mucin 16 (MUC16); TTN might have co-occurring mutations with low-density lipoprotein receptor-related protein 1B (LRP1B) and dedicator of cytokinesis protein 2 (DOCK2); MUC16 might have co-occurring mutations with ATP-binding cassette sub-family A member 13 (ABCA13); albumin (ALB) might have co-occurring mutations with Piccolo presynaptic cytomatrix protein (PCLO); ryanodine receptor 2 (RYR2) might have co-occurring mutations with WDR87, adhesion G protein-coupled receptor V1 (ADGRV1), and DOCK2; OBSCN might have co-occurring mutations with filaggrin (FLG); PCLO might have co-occurring mutations with WDR87 and ABCA13; ABCA13might have co-occurring mutations with CUB and Sushi domain-containing protein 3; DOCK2 might have co-occurring mutations with FLG; and ADGRV1 might have co-occurring mutations with WDR87. In the low APP.Score group, TP53 might have co-occurring mutations with KMT2D; MUC16 might have co-occurring mutations with RYR2; ALB might have co-occurring mutations with XIRP2; HMCN1 might have co-occurring mutations with RYR2; LRP1B might have co-occurring mutations with FBN2; and CTNNB1 might not have co-occurring mutations with AXIN1 and TP53. These results indicate that the number of gene pairs with simultaneous mutations would be greater in the high APP.Score group than in the low APP.Score group (Fig. 4A). In addition, there were differences in CNV frequency between the high and low APP.Score groups (Fig. 4B). These results indicate the significant involvement of APP.Score in both mutations and CNVs.

Fig. 4.

Mutation and copy number variation (CNV) analysis. (A) Correlation analysis between mutated genes. (B) CNV frequency.

We further analyzed the difference in drug sensitivity between the high and low APP.Score groups. In all of the drug sensitivity analyses, including temozolomide 1375, otx015 1626, leflunomide 1578, i-bet-762 1624, gdc0810 1925, dasatinib 1079, bpd-00008900 1998, bdp-00009066 1866, and azd5153 1706, the drug sensitivity of the high APP.Score group was higher than that of the low APP.Score group (Fig. 5).

Fig. 5.

Prediction of drug sensitivity.

We further analyzed the expression differences in KLRB1, PLA2G7, MSC, HSPA6, SPP1, and NCF2 between the tumor and normal groups using the TCGA-LIHC dataset. Compared with the LX-2 group, the expression of KLRB1, PLA2G7, MSC, HSPA6, and SPP1 but not NCF2 was significantly upregulated in the HepG2 and Huh-7 groups. The difference in KLRB1 and PLA2G7 expression at the gene and protein levels was the most significant (Fig. 6A,B). The OS rate was 0.5, and the survival time of patients with high expression of PLA2G7, MSC, HSPA6, and SPP1 was shorter, whereas that with high KLRB1 expression was longer (Fig. 6C). Therefore, PLA2G7 was selected for further study.

Fig. 6.

In vitro validation of gene expression of APP.Score. (A,B) The expression of genes (KLRB1, PLA2G7, MSC, HSPA6, SPP1, and NCF2) of APP.Score. (C) Survival curves of the six genes (KLRB1, PLA2G7, MSC, HSPA6, SPP1, and NCF2) used to construct the APP.Score were analyzed. Analysis of variance (ANOVA). *p 0.05, vs. LX-2.

Firstly, PLA2G7 was silenced in HepG2 and Huh-7 HCC cell lines (Fig. 7A,B). Compared to the si-NC group, the cell proliferation, migration, and invasion abilities of the si-PLA2G7 group were suppressed, while the level of cell apoptosis was significantly increased (Fig. 7C–G).

Fig. 7.

Inhibiting PLA2G7 affects the functions of hepatocellular carcinoma (HCC) cells HepG2 and Huh-7. (A,B) The expression of PLA2G7 was detected. (C) The cell proliferation (0, 12, 24, 48 h) was measured. (D) The cell proliferation ability was analyzed. (E) The cell migration was measured. (F) The cell invasion was analyzed. (G) The cell apoptosis level was assessed. ANOVA. *p 0.05, vs. the si-PLA2G7 group (HepG2); #p 0.05, vs. the si-PLA2G7 group (Huh-7).

We performed GSEA for PLA2G7, which was enriched for PD-L1 expression and the PD-1 checkpoint pathway in cancer, and the PI3K/AKT, and Janus kinase/STAT pathways; si-PLA2G7 inhibited PD-L1 expression (Supplementary Fig. 3 and Fig. 8A–C). The levels of PI3K, AKT, and p-NF-$\kappa$B/NF-$\kappa$B were elevated in the si-PLA2G7 group, while the level of p-STAT1/STAT1 was clearly diminished (Fig. 8D). The STAT1 activator, 2-NP, reversed the decreased PD-L1 expression caused by si-PLA2G7 (Fig. 8E), showing that PLA2G7 can affect PD-L1 expression via STAT1.

Fig. 8.

PLA2G7 affects PD-L1 expression via STAT1. (A) PLA2G7 was enriched in the PD-L1 pathway. (B) The expression of PD-L1 was detected. (C) PLA2G7 was enriched in the PI3K-Akt, JAK-STAT, and nuclear factor kappa B (NF-$\kappa$B) pathway. (D) The expression of PI3K, Akt, phosphorylated signal transducer and activator of transcription 1 (p-STAT1), STAT1, p-NF-$\kappa$B, and NF-$\kappa$B was measured. (E) The cells were treated with STAT1 activator (2-NP, 45 µmol/L) for 1 h. The fluorescence intensity of PD-L1 was observed. ANOVA. *p 0.05, vs. the si-NC group; #p 0.05, vs. the si-PLA2G7 group.