The Liver HCC dataset was downloaded from the UCSC Xena platform (https://xenabrowser.net). Sample selection from TCGA database employed these inclusion criteria: (1) histologically confirmed primary solid tumor specimens with matched adjacent non-tumor tissues; (2) transcriptomic profiling data derived from frozen patient samples. The final cohort comprised 369 tumor specimens and 50 matched normal controls. Additionally, a reference set of 79 peroxisome-associated genes was obtained from the KEGG database (https://www.kegg.jp). All datasets were processed under uniform inclusion criteria to ensure comparability across samples.

The bioinformatic analyses for biomarker discovery and validation were performed in accordance with the REMARK guidelines. DEGs between hepatocellular carcinoma and adjacent normal tissues were identified using the edgeR package in R. Intersection of these DEGs with the peroxisome-related gene set obtained from KEGG yielded a subset of pDEGs.

Prognostically significant pDEGs were identified through univariate and multivariate Cox proportional hazards regression implemented with R survival package. Resultant hazard ratios were visualized using forest plots. For refined variable selection, LASSO regression was applied with optimal penalty parameters determined via 10-fold cross-validation. To evaluate the predictive efficacy of selected pDEGs, TCGA survival records were curated by removing duplicate entries and randomly partitioning patients into training and validation sets at a 1:1 ratio. Time-dependent ROC analysis was conducted using survival ROC software, where larger AUC values were interpreted as enhanced predictive accuracy.

PPI networks were generated through the GeneMANIA web resource (http://genemania.org/). Functional annotation of gene sets was performed via Gene Ontology (GO) and KEGG pathway enrichment analyses using the ClusterProfiler package in R, with a statistical significance threshold set at P < 0.05.

Samples were stratified into high- and low-expression subgroups according to median Glyceronephosphate O-acyltransferase (GNPAT) mRNA levels. Associations between GNPAT expression and clinicopathological staging were visualized through boxplots generated with the ggpubr package in R.

Applying the ESTIMATE algorithm that provides the abundance of stromal and immune cells using transcriptomic profiling was used to characterizeTME composition in HCC specimens. The identification of key genes made the threads split into high- and low-expression groups using median values of expression to compare the stromal and immune scores. Gene expression signature-immune infiltration level correlation patterns were evaluated by applying the ESTIMATE pack in R. Moreover, relative fractions of 22 types of immune cells were deconvoluted using CIBERSORT computational method (https://cibersortx.stanford.edu/). Once the samples of HCC were partitioned into two expression-based clusters, the comparative distributions of populations of immune cells were represented in the form of violin plots created in the ggplot2 package in R.Together, these analyses enabled integrated quantification of stromal/immune content and immune cell composition across GNPAT-defined subgroups.

Tumour Immune Dysfunction and Exclusion(TIDE) computational model was used to determine possible reaction of HCCs to immune checkpoint blockade therapy. This is a methodology that measures tumor immune evasion based on a dysfunctional or exclusion mechanism of T cells which presents an estimate of likely sensitivity to checkpoint blockade. Further, Subclass Mapping was done to compare the progress of transcriptomics in the high and low-expression subgroups of GNPAT to reference samples of reported immunotherapy responders. This comparative method evaluated the presence of expression signatures in the individual subgroups that were found to be molecularly homologous to the signatures of the presence of the responses to either CTLA-4 [anti-cytotoxic T-lymphocyte-associated antigen 4] or PD-1 [anti-programmed cell death protein 1] therapeutic antibodies.

The data available in the Genomics of Drug Sensitivity in Cancer resource was subjected to systematic analysis of chemosensitivity relationships by comparing the expressions of GNPAT and half-maximal inhibitory concentrations (IC50) of a variety of conventional chemotherapeutic agents. R pRRophetic algorithm was used to make drug sensitivity predictions on all specimens that were using transcriptional profiles to predict therapeutic response. There was a comparative evaluation of the estimated IC50 values between subgroups of GNPAT high- and low-expression groups. Differential drug sensitivity was statistically evaluated by the Wilcoxon rank-sum test, providing a non-parametric assessment of differences between expression-based cohorts.

The analysis of lipid in biological samples was done using an optimized version of the Bligh-Dyer method, which enables the recovery of the largest amount of lipid in the organic layer and conservation of the integrity of volatile metabolic compounds. After extraction, the dried lipids film was resuspended in the isotope-labeled internal standard solution according to the next quantification and normalization required in the consequent analysis. A comprehensive screening of lipidomics in the Exion UPLC system was improved with a QTRAP 6500Plus mass spectrometer (Sciex, USA) together with electrospray ionization to provide a high sensitivity level and specific detection of the molecular forms. The polar lipids were chromatographically separated using the Phenomenex Luna silica column (150 x 2.0 mm, 3 m) with binary mobile phase system, mobile phase A (comprised of chloroform/methanol/ammonia 89.5: 10:0.5) and mobile phase B (comprised of chloroform/methanol/ammonia/water 55:39: 0.5: 5.5). Multiple reaction monitoring (MRM) was used to identify and determine lipid species with good accuracy in relative quantification and unequivocal ion discrimination.

In order to perform ultrastructural analysis, the 5×10⁷ cells of hepatocellular carcinoma were fixed, dehydrated, and embedded following a 3-step procedure involving an acetone gradient at ambient temperature. Thin slices (around 1μm) were stained with uranyl acetate then stained with lead citrate to increase the visibility of peroxisomes. TEM was used to perform imaging and random microscopic fields were chosen to be systematically evaluated with regard to peroxisomal morphology.

Paired tumor and adjacent non-tumorous liver tissues were obtained from HCC patients. Total RNA was isolated from all specimens following a standardized extraction protocol to preserve RNA integrity for downstream applications. Purified RNA was subsequently converted to complementary DNA using a commercial reverse transcription system (Toyobo, Japan). Gene expression quantification was performed via RT-qPCR, with relative transcript levels determined by the 2^−ΔΔCt method using endogenous reference genes for normalization. Primer sequences validated for target gene amplification are comprehensively detailed in Table 1.

Tissue sections were blocked with goat serum prior to incubation with a primary antibody targeting GNPAT. Following extensive washing, specimens were sequentially treated with a biotin-conjugated secondary antibody and horseradish peroxidase-streptavidin complex. Diaminobenzidine (DAB) served as the chromogenic substrate for visualization, with hematoxylin providing nuclear counterstaining. Stained sections were examined using a Leica microscope, and relative GNPAT protein abundance was quantified through staining intensity analysis performed with HistoQuest software.

The human HCC cell line HepG2 (ATCC, USA) was maintained in DMEM (HyClone, USA) supplemented with 10% FBS (HyClone, USA). To investigate GNPAT’s functional role in plasmalogen biosynthesis, cells were allocated into three experimental conditions: control, GNPAT-knockdown (si-GNPAT), and GNPAT-overexpression (OE-GNPAT) groups. For mechanistic studies examining plasmalogen-mediated tumor proliferation via PPARγ signaling, four treatment conditions were implemented: control, OE-GNPAT, OE-GNPAT with PPARγ antagonist T0070907 (OE-GNPAT+T007, 10 µM, MedChemExpress, USA), and si-GNPAT with hexadecyl plasmalogen supplementation (si-GNPAT+C16 P, using 10 µM C16:0 plasmalogen, Avanti Polar Lipids, USA).

Macrophage polarization assays utilized a Transwell co-culture platform, where M0 macrophages were co-cultured with HepG2 cells transfected with either GNPAT-overexpression constructs (OE-HepG2) or control vectors (NC-HepG2). Experimental configurations included: OE-HepG2+M0 (co-culture with GNPAT-overexpressing cells), OE-HepG2+M0+T007 (co-culture under PPARγ inhibition), and M0+C16 P (macrophages exposed to C16:0 plasmalogen alone). All cell lines underwent authentication through short tandem repeat profiling before experimental use.

Plasmenylethanolamine (pPE) isolation and structural characterization followed previously established protocols (28). Total lipid extracts were prepared, followed by phospholipid enrichment through cold acetone precipitation and subsequent silica gel column chromatography. Further purification of plasmalogen and ether phospholipid species was achieved using large-scale silica gel chromatography (280 × 800 mm). Structural identification of purified fractions involved LC-MS analysis, with compound verification based on retention behavior, accurate mass measurements, and characteristic tandem mass spectral fragmentation. Sample purity was assessed by comparing chromatographic profiles of isolated plasmalogens against a phosphatidylethanolamine (PE 14:0/14:0) reference standard using liquid chromatography-mass spectrometry.

For immunofluorescence studies, cells were plated on coverslips or in 24-well plates and exposed to various treatments. Fixation and permeabilization conditions were optimized according to subcellular localization of target antigens. Detection of Arg1 and ABCD3 employed fixation with 2% paraformaldehyde in PBS following the membrane permeabilization in 0.1% Triton X-100 (20 min at room temperature). The two-step protocol also facilitated antibody penetration. Following aldehyde quenching with 100 mM ammonium chloride for 10 min, the cells were labeled for PPARγ with pre-chilled acetone/methanol (1:1) for 20 min.

Following permeabilization, nonspecific binding sites were blocked by incubating the sections with 10% normal goat serum for 1 h. Primary antibody incubation proceeded overnight at 4 °C to ensure specific antigen recognition. After thorough washing, specimens were exposed to fluorophore-conjugated secondary antibodies (FITC or Alexa Fluor derivatives) for one hour under dark conditions. Nuclear counterstaining utilized DAPI (30 minutes) when required. Mounting was performed with 20% glycerol or Vectashield H1000 medium to preserve fluorescence integrity. Specificity controls included PBS-treated specimens, with final imaging conducted using fluorescence or confocal microscopy.

HepG2 cells were plated in 96-well plates (1×10⁴ cells per well) to ensure uniform distribution and attachment. Following a 24-hour incubation under standard culture conditions to establish monolayer integrity, cells were pulsed with 100 µL of 50 µM EdU for 2 hours to label replicating DNA. Cellular architecture and maintenance of nucleic acid integrity was ensured by subsequent fixation in 4% paraformaldehyde (20 minutes room temperature). The procedure used to detect incorporated EdU was as follows: the Cell-Light™ EdU Apollo^® 488 Imaging Kit (RiboBio, China) using DAPI, to visualize nucleus, to perform incorporation of EdU detecting plastic flow into nucleus. Fluorescence microscopic analysis was done after final washing procedures and indicated the successful detection of EdU positive nuclei, which indicated proliferating cell populations.

The assays were Transwell chamber assays which were used to measure the cell migration and invasion ability of the cells through porous membranes. For migration assays, serum-free medium containing 5×10⁴ cells was added to the upper chamber, and the bottom chamber was loaded with 20% FBS medium, which served as a chemoattractant. Invasion assays were conducted using Matrigel-coated membranes (BD Biosciences) to replicate physiological barriers, where these coated membranes were incubated 4for 4 h at 37°C before cell seeding. After a 24-hour incubation, the cells that had transmigrated to the lower surface of the membrane were fixed with 4% paraformaldehyde and then stained with crystal violet. Microscopic quantification of stained cells was used as a quantitative analysis of migratory and invasive potential.

Protein extraction from cultured cells utilized RIPA buffer (Thermo Fisher Scientific, China) following established protocols to ensure complete solubilization of cellular proteins. Protein concentrations were detected colorimetrically by a bicinchoninic acid (BCA) assay kit. For immunoblotting, equal protein aliquots (50 μg) underwent electrophoretic separation on 10% SDS-polyacrylamide gels, followed by transfer to PVDF membranes. Membranes were blocked with 5% skim milk in TBST for 90 minutes at 25°C before overnight incubation at 4°C with primary antibodies from Proteintech: GNPAT (14931-1-AP; 1:1,000), E-cadherin (20874-1-AP; 1:10,000), N-cadherin (22018-1-AP; 1:2,000), Vimentin (10366-1-AP; 1:20,000), and PPARγ (16643-1-AP; 1:1,000). After washing, membranes were probed with HRP-conjugated secondary antibodies (CW0103/CW0110S; CWBio; 1:1,000) for 2 hours. Protein bands were visualized via enhanced chemiluminescence (Thermo Fisher, USA) and quantified using ImageJ software, with β-actin serving as loading control for normalization.

Cytokine concentrations (TGF-β, IL-10, VEGF) in culture supernatants were quantified with commercial ELISA kits (JRDUN Biotechnology, China). Absorbance readings were normalized to total protein content determined by the Pierce BCA Protein Assay (Thermo Scientific, USA), with final results expressed as nanograms per gram of total protein.

CD206 surface expression on co-cultured macrophages was quantified through flow cytometric analysis, with untreated M0 macrophages serving as negative controls. Following detachment using cell scrapers and resuspension in chilled PBS with 2% FBS, cells were pretreated with human Fc receptor blocking reagent to minimize nonspecific antibody binding. For specific detection of the M2 polarization marker CD206, cells were stained with PE-conjugated anti-human CD206 monoclonal antibody for 30 minutes at 4°C under light-protected conditions. After thorough washing to remove unbound antibodies, samples were analyzed on a BD FACSCanto II system (BD Biosciences, USA), enabling precise quantification of CD206-positive populations based on fluorescence emission profiles.

R (v4.3.1) was used to carry out all statistical calculations and data representation, which is a complete environment in the application of analytical processes. Enrichment analysis p-values were adjusted using BenjaminiHochberg adjustment in order to control the detection of false findings. Wilcoxon rank-sum test was used as a non-parametric distribution comparison tests as a group comparison. Further statistical analysis and visualization were done using GraphPad Prism 9 (GraphPad Software). To demonstrate the quantitative data of 3 or more independent replicates, the result is represented as mean with standard deviation. Two-group comparisons were using Student t-tests (independent or paired design), and multi-group comparisons were performed using the one-way ANOVA with post-hoc tests of Dunnett. The analysis of the flow cytometry data and the population quantification was done using the FlowJo software (Tree Star). Statistical significance was defined as p < 0.05.

In-depth review of the TCGA-LIHC cohort showed that there were 12, 718 differentially expressed genes in hepatocellular carcinoma which included 9, 746 up-regulated and 2, 972 down-regulated transcripts (Figure 1A). Intersection of the peroxisome-associated gene set with the DEGs revealed 38 pDEGs, the expression patterns of which were decomposed using hierarchical clustering (Figures 1B, C). GNPAT was observed to be the most important prognostic factor for overall survival by univariate Cox and LASSO regression analyses (Figures 1D, F, H). The predictive model was shown to be very valid with values above 0.6 on the ROC analysis (Figures 1E, G). The network analysis of protein-protein interaction showed that GNPAT is associated with the key lipid-metabolic enzymes, such as AGPS, FAR1, and PTS (Figure 1I). GNPAT was functionally enriched and associated with peroxisomal pathways, glycerophospholipid metabolism and ether lipid biosynthesis (Figures 1J, K). Clinical validation demonstrated that GNPAT increased progressively as the TNM stage advanced with a peak of TNM stage III (Figures 1L, M).

Stratification based on median GNPAT expression revealed significant alterations in the TME. While tumors with high GNPAT expression exhibited elevated immune scores (Figure 2A), further analysis indicated that this was driven by a distinct leukocyte composition, characterized by significant enrichment of immunosuppressive cell types, particularly M0 macrophages and neutrophils (Figure 2B). This pattern of infiltration is associated with impaired anti-tumor immunity, which aligns with the significantly higher TIDE scores observed in the high-GNPAT group (Figure 2E), predicting greater immune evasion and poorer response to immunotherapy. Analysis of immunophenoscore suggested a possibility of response to combined PD-1/CTLA-4 blockade in specific subsets (Figure 2C), although no significant difference was found in the expression of conventional immune checkpoint molecules (Figure 2D). Consistent with an altered metabolic state, pharmacogenomic analysis predicted increased sensitivity to several therapeutic agents, including Sorafenib, Cyclopamine, Embelin, PAC-1, and AKT inhibitor VIII, in high-GNPAT tumors (Figure 2F).

Clinical specimen multimodal validation showed there were significant peroxisomal changes in HCC. The transmission electron microscopy showed a higher density of the peroxisomes with fewer organellar diameters in tumor tissues (Figures 3A–C). Lipidomic profiling of 20 matched pairs revealed a complex alteration in plasmalogen species. Most plasmalogen species, including both plasmenylethanolamine (pPE) and plasmenylcholine (pPC) subtypes, showed an overall decrease in tumor tissues (Figures 3D, E). This suggests a general reduction rather than a selective remodeling of the plasmalogen pool in HCC. The expression of GNPAT in malignant tissues was always corroborated by transcriptional and protein analysis (Figures 3F, G).

Efficient GNPAT modulation in HepG2 cells was confirmed through transcriptional and protein assessment (Figures 4A–C). Immunofluorescence analysis demonstrated that GNPAT overexpression enhanced peroxisomal abundance while knockdown produced opposite effects (Figures 4D, E). LC-MS/MS quantification revealed that GNPAT elevation significantly increased pPE levels (Figure 4F), supporting a central role for GNPAT in regulating plasmalogen biosynthesis.

Mechanistic investigation revealed PPARγ pathway activation upon GNPAT overexpression, characterized by enhanced protein expression and increased nuclear accumulation of PPARγ (Figures 5A–C). This observed nuclear accumulation is consistent with pathway activation and is further corroborated by the functional rescue upon pharmacological inhibition of PPARγ. Functional assays demonstrated that GNPAT potentiated proliferation, migration, epithelial-mesenchymal transition, and apoptosis resistance (Figures 5D–M). Specifically, GNPAT overexpression significantly enhanced cell proliferation (Figure 5D), promoted migratory and invasive capacities (Figures 5E–G), and induced a shift towards a mesenchymal phenotype, as evidenced by decreased E-cadherin and increased N-cadherin/Vimentin expression (Figures 5J–M). Conversely, GNPAT knockdown or PPARγ inhibition with T0070907 consistently reversed these pro-tumorigenic phenotypes, confirming the dependency of these processes on the GNPAT-PPARγ axis.

To investigate the immunomodulatory effects of GNPAT in the tumor immune microenvironment, a co-culture system was established that would include cells of HCC and THP-1-differentiated macrophages. The immunofluorescence analysis revealed substantial upregulation of arginase-1 (Arg1) which is one of the canonical M2 polarization markers in macrophages in co-culture with GNPAT-overexpressing HepG2 cells relative to the control conditions (Figure 6A). Analysis of PPARγ signaling revealed that macrophages upon observation of high-GNPAT HCC cells exhibited increased transcriptional levels of PPARγ and its downstream metabolic targets ACOX1 and CD36 (Figures 6B–D). Direct administration of C16:0 plasmalogen to macrophages was sufficient to induce PPARγ pathway activation and M2 polarization, supporting the role of plasmalogens as key signaling mediators in this cross-talk. Cytokine profiling of the co-culture supernatants revealed increased release of immunosuppressive factors (TGF-β, IL-10) and the angiogenic mediator VEGF (Figure 6E). Stunningly, all the above effects were inhibited by the PPARγ antagonist T0070907 and demonstrated that PPARγ is the primary regulator of this immunometabolic axis.