Three GEO datasets, including GSE133039, GSE131329, and GSE81928, were applied for data mining. A brief description of each and their uses are shown in Supplementary Table 1. The clinical characteristics of these datasets are presented in Supplementary Tables 2-3. All transcriptome data was standardized by log2 (FPKM+1) transformation.

An FR gene set was constructed utilizing 253 homo-sapiens ferroptosis genes from FerrDb (http://www.zhounan.org/ferrdb/current/), recognized as the most comprehensive database dedicated to ferroptosis regulators and ferroptosis-disease associations (14). The protein-protein interaction (PPI) network of these FR genes was depicted through the STRING database (https://string-db.org/) (15) and Cytoscape v3.71 software (16). Their biological functions were analyzed using the online Metascape tool (http://metascape.org/) (17).

Consensus clustering analysis was employed for identifying HB diagnostic patterns based on the expressions of FR genes. This analysis was accomplished by using the ‘ConsensusClusterPlus’ package in R v4.1.2 (18). The optimal number of clusters could be determined according to the K-means algorithm and cumulative distribution function (CDF) (18). A heatmap of the consistency matrix was applied to visualise the intra- and intergroup variabilities within each cluster.

In the present study, two approaches were combined to identify the core FR genes for characterizing HB, including Lasso regression analysis and the support vector machine recursive feature elimination (SVM-RFE) algorithm. The Lasso regression analysis was conducted using the ‘glmnet’ R package based on the 7-fold cross-validation scheme. SVM-RFE, a common machine-learning method, was used to reduce the number of parameters by creating weight vectors for each parameter. These weights determine the importance of the parameters, thereby maximizing data separation via an arithmetical hyperplane (19). SVM-RFE was performed using the ‘KeBABS’ R package (20).

Considering that only the GSE131329 dataset harbored HB survival outcomes, yet without survival time, the differences in expression of three candidate genes between deceased and surviving cases were compared for this dataset. Their predicted accuracy was also assessed using the receiver operating characteristic (ROC) curve. The relationships between the expression of the three core genes and HB clinicopathological features were analyzed based on the analysis of variance test, including histological type, metastatic event, and clinical stage.

The changes in infiltration levels of 21 immune cells were evaluated through the CIBERSORT algorithm, which quantifies the cell composition of sample tissue based on the gene expression profiles (21).

To verify the immunogenic effects of NOS2, we performed immunofluorescence detections on two pairs of clinical samples with different NOS2 protein levels. Immunofluorescence was conducted as previously described (22). Briefly, 4-mm tissue sections of HB clinical samples were prepared, and dewaxed with xylene and anhydrous ethanol. Then, the tissue sections underwent antigen repair, fixation, and nonspecific blocking. Next, the sections were incubated with first and secondary antibodies (NOS2, ab178945; CD8, ab33786; CD163, ab87099; Abcam, Cambridge, UK) in conditions protected from light. The slides were then sealed with anti-fluorescence quenching sealing tablets. CD8 and CD163 were stained with the Cy3-labeled goat anti-mice IgG (H&L) antibody (P0193; Beyotime, Shanghai, China). NOS2 was stained with the FITC-labeled goat anti-rabbit IgG (H&L) antibody (P0186; Beyotime). Nuclei were stained with Hoechst reagent (C0003; Beyotime). Subsequently, the slides were analyzed with an automatic fluorescent microscope using a 40× objective lens (BX53; Olympus, Tokyo, Japan).

GSEA was used to probe into the effects of NOS2 on multiple metabolic processes, including glycolysis, fatty acid metabolism, and cholesterol biosynthesis. The corresponding gene sets were obtained from the MsigDB database (https://www.gsea-msigdb.org/) and their detailed descriptions are presented in Supplementary Table 4. The parameters for GSEA analyses were set as follows: phenotype labels set as high- versus low-NOS2 expression samples; the number of permutations 1000; and no collapse in gene symbols. Transcriptome data were derived from the GSE133039 dataset.

The Genomics of Drug Sensitivity in Cancer (GDSC) database (www.cancerRxgene.org) was utilized for the connections between NOS2 and the sensitivity of both molecular targeted therapy (MTT) and chemotherapy drugs (23). Drug sensitivity was based on the half-maximal inhibitory concentration (IC₅₀). Correlations were assessed based on Pearson’s coefficient. Mutational features of the three candidate genes in the two pediatric projects (n=1,064 samples) were displayed using the cBioPortal database (http://cbioportal.org) (24).

The TRRUST database (https://www.grnpedia.org/trrust/) (25), hTFtarget database (http://bioinfo.life.hust.edu.cn/) (26), and a subsection of the GeneCards database termed ‘Promoters/Enhancers’ (https://www.genecards.org/) (27) were employed for predicting the potential upstream transcription factors (TFs) of NOS2. The promoter sequence of NOS2 was obtained from the UCSC database (http://genome.ucsc.edu/). The possible binding sites between the predicted TFs and NOS2 were achieved through the JASPAR database (https://jaspar.genereg.net/) (28). The motif sequence of the candidate TF was also exhibited through the JASPAR database. Potential enhancers of NOS2 were predicted using the EnhancerDB database (http://lcbb.swjtu.edu.cn/EnhancerDB/) (29).

RT-qPCR analyses were conducted on 25 pairs of HB tissues and adjacent normal tissues from the Second Affiliated Hospital of Xi’an Jiaotong University to confirm differential expression of NOS2. Prior to the clinical sample analysis, informed consent was obtained from the parents of the children involved. The study protocol was approved by the Ethics Committees of the Second Affiliated Hospital of Xi’an Jiaotong University. All methods in this study were performed in accordance with relevant guidelines and regulations.

For RNA extraction, TRIzol reagent (Takara, Dalian, China) was employed. The purity of each RNA sample was assessed using the optical density ratio at 260nm and 280nm (A260/A280). Reverse transcription was performed using the PrimeScript RT reagent Kit (Takara). The RT-qPCR reaction was tracked using the SYBR-Green PCR Reagent (Takara) and subsequently detected on an ABI Prism 7900 sequence system. GAPDH was used as the internal reference gene. The relative gene expression was calculated according to the 2^-ΔΔCT method. The list of primer sequences is shown in Supplementary Table 5.

Western blot assays were performed as previously described (10). Transfected cells were lysed using RIPA buffer (Beyotime) while kept on ice. The protein concentration of lysastes was measured using a BCA kit (Phygene Life Sciences Company, Fuzhou, China). Protein samples were separated using 10% SDS-PAGE. After electrophoresis, protein samples were transferred to PVDF membranes (BestBio, Shanghai, China). The PVDF membranes were blocked using 5% skim milk at 37°C for 2 h. The membranes were subsequently incubated with the primary and secondary antibodies, respectively (all purchased from Abcam). The primary antibodies were as follows: rabbit monoclonal anti-NOS2 (1:1,000; ab178945) and rabbit polyclonal anti-GAPDH antibody (1:1,000; ab22555). The secondary antibody was goat-anti rabbit IgG-HRP secondary antibody (1:2,000; ab205718).

Two HB cell lines (HepG2 and HuH-6) were purchased from Procell Life Company (Wuhan, China). These cells were maintained in distinct culture mediums: HepG2 cells were cultured in MEM medium, while HuH-6 cells were cultured in DMEM medium. Both mediums were supplemented with 10% fetal bovine serum (FBS) and 1% penicillin/streptomycin. The cells were cultured at 37 °C and 5% CO₂.

For gene manipulation, short hairpin RNA (sh-NOS2) for NOS2 knockdown and overexpression plasmids (OE-NOS2) were designed and synthesized by HanHeng Biotechnology (Shanghai, China). The transfection process followed the manufacturer’s instructions.

The colony formation assay was used to evaluate the proliferative abilities of HB cells. Transfected cells in the logarithmic growth phase were seeded into 6-well plates (1×10³ cells per well). After 1-week incubation, visible colonies formed. Then, the cells were fixed and stained using methanol and Giemsa, respectively. Colonies were counted under a microscope by examining five random fields.

The workflow was consistent with that described in previous research (10, 30). For migration assays, mediums containing 0.1% FBS were added to the upper chambers, while medium containing 10% FBS was added to the lower chambers. Cells were seeded in the upper chamber at a density of 1×10⁴ cells per well. After 24 h, non-migrated cells were removed washing with PBS and gently scraping with a cotton swab. The cells adhering to the lower membrane surface represented migrated cells. These cells were fixed by paraformaldehyde for 20 min and stained with 0.1% crystal violet for 20 min. Cell counting was performed under a high magnification microscope (100-fold) from five randomly-selected visual fields.

For invasive assays, the upper chambers were precoated with Matrigel (CORNING, Corning, NY, USA). The subsequent steps were similar to the migration experiment.

The flowchart of this study is presented in Figure 1. First, we selected three potential core ferroptosis genes in HB using Lasso regression analysis and the SVM-RFE algorithm based on a comprehensive and reliable FR gene set (n=253). Second, their clinical relevance was explored through three cohorts. Third, NOS2 was selected for further multi-omics investigation, including its impact on immune responses, its associations with HB metabolism, and its role in determining the therapeutic efficacy of both MTT and chemotherapy agents. Moreover, we explored the functional implications of NOS2 in driving the malignant progression of HB.

A reliable gene set is a prerequisite for finding key genes regulating cancer development (31). With the aid of the FerrDb database, we included 253 human ferroptosis genes in the FR gene set. The PPI network of these genes is exhibited in Figure 2A. Functional analysis revealed that these regulators were closely involved in multiple pivotal sessions of ferroptosis, such as the reactive oxygen species metabolism process, cellular response to oxidative stress, and antioxidant activity (Figure 2B). KEGG analysis also showed that the selected genes were highly enriched in the ferroptosis process (Figure 2C). These findings confirmed the reliability of the FR gene set. Through module analysis, 36 core ferroptosis genes were obtained for further analyses (Figure 2D).

Given that only the GSE131329 dataset provided the survival status of patients with HB among the three GEO datasets, we used it as the training cohort to identify the critical regulators involved in HB progression and clinical outcome. Before further bioinformatic analyses, the difference in transcriptome data of FR genes between normal and tumor samples was first determined (Figure 3A). Notably, 15.8% of FR genes (40/253) were differentially expressed in the HB samples. Differentially expressed genes (DEGs) in ferroptosis revealed two cluster patterns for characterizing HB (Figures 3B–D). The ferroptosis clusters not only demonstrated an excellent diagnostic capacity for distinguishing between tumor and normal cases (AUC=0.953, Figure 3E), but also had substantial potential for predicting the survival outcomes of patients with HB (AUC=0.854, Figure 3F).

Figure 3G displays the heatmap of FR DEGs based on distinct survival statuses. To identify key FR genes in HB pathogenesis, we performed Lasso regression analysis and applied the SVM-RFE algorithm using survival status as the outcome variable (Figures 3H, I). The genes that emerged as common outcomes from both machine-learning methods were identified as the intersection (Figure 3J). From this selection, NQO1, NOS2, and SLC1A4 were screened for further clinical correlation analyses.

The basis for functional gene activity lies in their expressive differences. Hence, we confirmed the expressions of NQO1, NOS2 and SLC1A4 in HB. All three candidate genes were overexpressed in tumor samples in the GSE133039 and GSE131329 datasets (Figures 4A, B). However, there was no difference in SLC1A4 expression between normal and tumor samples in the GSE81928 dataset (Figure 4C). ROC curve analyses indicated that the expression of these genes could differentiate between normal and HB samples in most datasets (Figures 4D-F), with an average diagnostic accuracy of ~0.80 (Figure 4G). Specifically, NOS2 emerged as the best diagnostic biomarker with the highest AUC value of 0.823 and the lowest standard difference of 0.007 (Figure 4G).

Further, we investigated the relationships between three candidate genes and HB clinical status. In the GSE131329 dataset, NOS2 had moderate competence for predicting survival outcomes of patients with HB, which was slightly higher than NQO1 and SLC1A4 (AUC=0.649, Figure 5A). There were no differences in the expressions of NQO1 and SLC1A4 between different survival outcomes, whereas higher NOS2 expression was observed in surviving cases compared to deceased cases (Figure 5B). As for histological types, no differences in the expressions of the three genes were found between poorly- and well-differentiated patterns (Figure 5C). Moreover, higher NOS2 expression was concordant with metastatic events and late clinical stage in the GSE131329 dataset, while no correlations were found for NQO1 (Figures 5D, E). Collectively, although NQO1, NOS2, and SLC1A4 were all aberrantly expressed in HB, only NOS2 was closely associated with clinicopathological features. Therefore, we further investigated the effects of NOS2 on the tumor immune microenvironment (TIME) and metabolism in HB.

Through using the CIBERSORT algorithm, we found that different NOS2 expression resulted in alteration of immune cell abundance (Figures 6A–C). For instance, in the GSE133039 dataset, the infiltration levels of CD8+ T cells were significantly increased in the high-NOS2 expression group compared with the low expression group, whereas that of macrophages and naïve B cells was markedly increased (Figure 6A). An overview of the immune analyses findings is summarized in Figure 6D. Notably, CD8+ T cells and macrophages were the dominant cell types affected by NOS2 expression. Given the functions of these cells in anti-tumor immune response and their variation in HB samples with high-NOS2 expression (Table 1), it became apparent that high-NOS2 expression indicated an unfavorable anti-tumor immune process.

To verify the effects of NOS2 on the TIME in HB, we performed immunofluorescent staining on HB clinical samples with different NOS2 expression levels. As shown in Figure 7, the green fluorescence intensity of labelled NOS2 was consistent with NOS2 expression. However, the red fluorescence intensities of labelled CD8+T cells and macrophages exhibited vastly different levels. In samples with high-NOS2 expression, the red fluorescence of CD8+T cells was rarely visible, whereas that of macrophages filled the field of view. Inversely, in samples with low-NOS2 expression, the fluorescence intensity of CD8+ T cells was strong, whereas that of macrophages was extremely weak. These findings confirmed that NOS2 expression acted as a biomarker of immune status in HB.

Metabolic reprogramming, particularly characterized by enhanced glycolysis, is a critical hallmark of cancer cells (32, 33). Hence, we investigated the associations of NOS2 with three metabolic pathways, including glycolysis, fatty acid metabolism, and cholesterol biosynthesis. Interestingly, using GSEA analysis, there were no enrichment differences in these metabolic processes between samples with high- and low-NOS2 expression (Figures 8A–C). In light of the above, it was speculated that alterations in NOS2 expression might not significantly impact the metabolic processes within HB metabolic.

Using the GDSC database, the potential association between NOS2 expression and the therapeutic effects of multiple drugs were explored. For MTT, NOS2 expression was not correlated with the sensitivity of most MTT drugs; however, it was positively correlated with that of trametinib and lapatinib (Figure 9A). Similarly, no correlations were observed between NOS2 expression and commonly used chemotherapeutic drugs, such as 5-fluorouracil, doxorubicin, and vinblastine (Figure 9B). Only the IC₅₀ value of cisplatin was positively correlated with NOS2 expression (R=0.245, P=0.037). Therefore, high expression of NOS2 was indicative of resistance to trametinib, lapatinib, and cisplatin treatments.

Moreover, the mutational characteristics of several crucial genes (NQO1, NOS2, SLC1A4, NFE2L2, KEAP1, and BACH1) were briefly elaborated through the cBioportal database. The somatic mutations of these genes rarely appeared in pediatric cancer samples with a mutational frequency of ~0.3% (Figure 9C). More specifically, we also focused on their mutational status in HB samples (16 cases from the ‘Columbia U’ project and 0 cases from the ‘DKFZ’ project). As shown in Supplementary Table 6, only 1 HB sample was concomitant with the NFE2L2 missense mutation, whereas any mutations of the other genes were not found in all 16 HB samples.

The three-dimensional protein structure of NOS2 was obtained through the AlphaFold database (Figure 9D), which provided some clues for artificially synthesizing a NOS2 inhibitor. Additionally, the Hum-mPLoc online tool revealed that the cytoplasm was the main subcellular localization of NOS2 (Figure 9E). Functional analysis indicated that genes highly correlated with NOS2 (R≥0.5) may participate in multiple biological processes, such as ‘response to hypoxia’ and ‘GTPase binding’ (Figure 9F).

NOS2 acts as a biological mediator in human diseases by generating nitric oxide (NO). Numerous studies have demonstrated its roles in various types of cancers (Table 2). Through the above analyses, we found that NOS2 was not only closely associated with clinical features of HB, but also with its prognosis, anti-tumor immunity, and therapeutic effects. Therefore, we further ascertained the biological functions of NOS2 in HB cells.

Through RT-qPCR (n=25) and Western blot (n=4) analysis on clinical samples, we successfully verified the upregulation of NOS2 in HB (Figures 10A, B). Further experiments involving transfection with specific plasmids (i.e., sh-NOS2 and OE-NOS2), effectively altered NOS2 expression in HepG2 and HuH-6 cells (Figure 10C). Colony formation assays showed that NOS2 knockdown significantly inhibited cell clone formation, whereas NOS2 overexpression markedly promoted this process (Figure 10D). Similar results were also found in the quantitative analyses (Figures 10E, F). Hence, NOS2 could enhance the proliferative abilities of HepG2 and HuH-6 cells.

Next, we verified the effects of NOS2 on the migration and invasion of HB cells through Transwell assays. In HepG2 cells, the migrative and invasive cells were significantly increased in the overexpression group (OE-NOS2) compared to the control group (Figures 11A, B). However, migration and invasion were markedly decreased in the NOS-silenced group (sh-NOS2) (Figures 11A, B). The quantitative analyses, based on cell counts, were also consistent with the above findings (Figures 11C, D). Similarly, in HuH-6 cells, NOS2 overexpression increased migration and invasion of cells, whereas NOS2 deletion exerted an inhibitory effect (Figures 11E–H). Collectively, NOS2 was shown to possess cancer-promoting abilities within HB cells, significantly influencing their migration and invasion capacities.

Using the TRRUST, hTFtarget, and GeneCards databases, we predicted upstream TFs of NOS2.By merging the results using a Venn diagram, we identified serval candidate TFs, including FOS, JUND, STAT1, and YY1 (Figure 12A). The JASPAR database was used to quantify the probabilities of these TFs binding to the NOS2 promoter via predictive scores (Figure 12B). Among these, STAT1 had the highest predictive score of 12.366, indicating that STAT1 may play a role in transcriptional regulation of NOS2. Curiously, only JUND had a positive expression correlation with NOS2, while STAT1 demonstrated a negative correlation (Figure 12C). There were no significant differences in the expressions of FOS, STAT1, and YY1 between patients with HB that had deceased or surviving outcomes, except for JUND (Figure 12D). Therefore, despite STAT1 having the highest predictive score, JUND seemed to be more likely to regulate NOS2. The motif sequence of JUND’s potential binding to NOS2 is presented in Figure 12E. The binding site was identified as 5’-TTCTGACTCTTTT-3’, located between the 458^th base to the 446^th base positions upstream of the NOS2 transcription start site (TSS) (Figure 12F).

Of note, the JASPAR database can only infer the binding possibility between TFs and the promoter regions. It fails to predict interactions between TFs and enhancer regions. Therefore, the above analyses could not exhibit a comprehensive view of the upstream regulatory landscape of NOS2. Thus, to address the above shortfall, we preliminarily investigated the enhancers that regulated NOS2 through the EnhancerDB database (Figure 12G). The results revealed that a total of 12 enhancers may be involved in NOS2 regulation. The distance between enhancers and the NOS2 TSS ranged from 3,040 kb to 88,494 kb. This approach offers a broader perspective on the potential regulatory elements influencing NOS2, capturing both promoter-based and enhancer-based interactions.

To date, some studies have developed new biomarkers for HB clinical assessments. In comparison with these investigations, NOS2 offers several distinct advantages (Table 3). First, NOS2 exhibits the excellent diagnostic performance. Recently, Liu et al. and Jiao et al. have reported that exosomal miR-21 and miR-34 were upregulated in HB samples, thereby acting as potential diagnostic biomarkers (34, 35). NOS2 diagnostic accuracy is comparable to that of above biomarkers (0.823 versus 0.861 and 0.837). Second, NOS2 comprehensively demonstrates the clinical significance of biomarkers. Previous studies primarily focused on establishing the relationships between biomarkers and HB prognosis. For instance, Zhao et al. exclusively confirm the expression pattern of apolipoprotein A-I in HB without performing any clinical correlation analysis (36). By contrast, our study not only investigated the prognostic value of NOS2, but also elaborated its associations with immune response, metabolism, and therapeutic effects. Third, NOS2 provides insights into the biological functions of biomarkers. In the present study, we employed a series of in vitro experiments to ascertain the effects of NOS2 on HB proliferation, migration, and invasion. These findings provide new clues for the identification of novel therapeutic targets. In summation, despite the existence of identified HB biomarkers, NOS2 significantly enriches our understanding of HB clinical assessments, offering valuable information with its distinct advantages.