We obtained the original data and corresponding clinical information of 371 liver hepatocellular carcinoma (LIHC) cases from The Cancer Genome Atlas (TCGA) database (https://portal.gdc.com). A total of 42 NRGs were identified from a previous comprehensive review (7). Four NRGs, ADRB3, CHRM1, CHRM2, and CHRNA9, were eliminated since they were not expressed in LIHC. The R software package ConsensusClusterPlus (v1.54.0) (8) was applied for consensus analysis. The PAC structure of the tool identifies the default cluster number and is repeated 100 times to extract 80% of the sample, clusterAlg = “hc”, and innerLinkage= ‘ward. D2’. Heat map clustering was analyzed using the R software package pheatmap (v1.0.12). The genes with a variance above 0.1 were retained in the gene expression heat map. If the number of input target genes was more than 1,000, the top 25% genes were extracted and displayed, after sequencing, from a large to small variance.

Two subtypes were gained following a consensus analysis, and clinical information was downloaded. For CHRNE, GFRA2, GRIN2D, cutoff-high (top 25%) and cutoff-low (bottom 25%) and for GFRA3, cutoff-high (25%) and cutoff-low (25%) were used as thresholds. R software package (v4.0.3) ggplot2 (9) and pheatmap were applied. Significance p-values were analyzed by the chi-square test, where the size of the value was taken -log10 (p-value), and, if marked with *, it represents a significant difference in the distribution of this clinical characteristic in the corresponding two groups (p < 0.05).

The R software package Limma (v3.40.2) was used to study the differential expression of mRNA (10). Adjusted p-values were analyzed in TCGA or GTEx to correct for false- positive results. Adjusted P<0.05 with log2(FC) (multiple change)>1 or log2(FC) (multiple change) <-1 was defined as the threshold for the differential expression of mRNA. The differential expression analysis between subtypes was commonly used: C1 vs. C2. CHRNE, GFRA2, and GRIN2D used cutoff-high (top 25%) and cutoff-low (bottom 25%) as expression thresholds. GFRA3 used cutoff-high (25%) and cutoff-low (25%) as thresholds, since they had significantly low expression in LIHC.

The R software package ClusterProfiler (11) was used to carry out Gene Ontology (GO) (12), and the Kyoto Encyclopedia of Genes and Genomes (KEGG) (13) was used for the analysis of potential mRNA.

Based on the clinical information from 371 LIHC cases from the TCGA database, we generated Kaplan–Meier (KM) plots with log-rank P-value and time-dependent receiver operating characteristic (ROC) analysis to compare the predictive accuracy and risk score of 38 neural-related genes. The least absolute shrinkage and selection operator (LASSO) regression algorithm was used for feature selection, and 10-fold cross validation was used, lambda. min=0.071 (14). The R software package glmnet was used for the above analysis. For KM curves, p-values and hazard ratios with a 95% confidence interval (CI) were determined using a log rank test and univariate Cox proportional hazard regression. All the above analyses were performed using the R software package (v4.0.3). P<0.05 was considered as a statistical difference.

The one-class logistic regression (OCLR) algorithm, invented by Malta et al., was used to evaluate the stemness of mRNA (15).

CIBERSORT and EPIC from the R software package immunedeconv (https://grst.github.io/immunedeconv) were applied to analyze the immune infiltration of different subtypes (16). The R software packages (v4.0.3) ggplot2 and pheatmap were used for visualization.

We tested the correlation of specific neural-related gene expression and 8 commonly used immune checkpoint genes (ICGs) (e.g., CD274, CTLA4, HAVCR2, LAG3, PDCD1, PDCD1LC2, SIGLEC15 and TIGIT). The R software packages (v4.0.3) ggplot2 and pheatmap were used for visualization. The TIDE algorithm was applied to predict the therapeutic response to immune checkpoint blockades (ICBs) (16).

The data of 38 NRGs in LIHC were collected from the TCGA database. As mentioned previously, ADRB3, CHRM1, CHRM2, and CHRNA9 were eliminated since they were not expressed in LIHC. For identification of subtypes, the R software package ConsensusCluster Plus (V1.5.4.0) was applied to obtain 4 subtypes: group 1 (C1), group 2 (C2), group 3 (C3), and group 4 (C4) (Figures 1A–C). The groups C3 and C4 were excluded from this study due to limited sample size, which might lead to a confined biological or clinical value. Ten years post-follow-up, the overall survival (OS) and progression-free survival (PFS) of the four groups had no statistic difference (P>0.05) (Figure 1D). However, the 5-year OS and 3-year PFS of C1 were better than C2 (P=0.03; P=0.03) (Figure 1E). They could be attributable to the withdrawal of C1 during the sixth-year follow-up. In addition, C1 had a better T category, clinical staging, and pathological grading with a significant statistical difference (P<0.05) (Figure 2 and Table 1).

We used the R software package Limma (v3.40.2) to analyze differentially expressed genes in C1 and C2. Those with FC>2 and P<0.05 were interpreted as differentially expressed genes (DEGs). Compared to C2, C1 had 622 downregulated and 262 upregulated DEGs. Among 38 NRGs, C1 possessed 27 downregulated (e.g., CHIN3A), 5 upregulated (e.g., CHRNA4), and 6 unregulated (e.g., CHRM3) DEGs compared to C2 (Figures 3A, B, S1). Meanwhile, KEGG and GO analyses were carried out to further assess signaling pathways in C1 and C2 (Figures 3C, D). KEGG analysis revealed that compared to C2, C1 showed the activation of signaling pathways for liver metabolism, including bile secretion, cholesterol metabolism, glycolysis/gluconeogenesis, and chemical carcinogenesis (DNA adducts and receptor activation), while a suppression of cancer-associated signaling pathways such as bladder cancer, the cell cycle, central carbon metabolism in cancer, and PI3K-Akt signaling pathway (Figure 3C). GO analysis suggested that compared to C2, C1 had activated the liver metabolic process (e.g., alcohol metabolic process, fatty acid metabolic process) and a suppression of growth factor beta stimulus, negative regulation of cell adhesion, extracellular matrix organization, which were regarded as critical processes in promoting cancer growth and metastasis (Figure 3D).

For the evaluation of the immune infiltration of C1 and C2, we used immunedeconv, an R software package based on the integration of CIBERSORT, EPIC, MCP-counter, quanTIseq, TIMER, and xCell (16). During the analysis, we used CIBERSORT and EPIC, which were the most used algorithms. CIBERSORT analysis revealed significant differences in the naïve B cell (P<0.01), memory B cell (P<0.001), regulatory T cell (Tregs)(P<0.001), monocyte (P<0.05), and macrophage M0 (P<0.001), suggesting that C1 exhibited stronger immunosuppression when compared with C2 (Figure 4A). EPIC further confirmed that C1 and C2 had significant difference in macrophage infiltration (P<0.001); however, there were insufficient data on macrophage subtypes (Figure 4B). We also applied R software packages ggplot2 and pheatmap to analyze ICGs in the two subtypes and showed that CTLA4, HAVCR2, LAG3, PDCD1, SIGLEC15 (P<0.001), and TIGIT (P<0.01) were downregulated in C1 when compared to C2 (Figure 4C). The TIDE algorithm, used to predict cancer immune response, was also used herein (16). The TIDE score was higher in C2 than C1, with significant difference (P<0.001), indicating that C1 might achieve more clinical benefit upon ICBs (Figure 4D). The OCLR algorithm showed that the stemness of C1 and C2 had no statistical difference (P>0.05) (16) (Figure 4E).

We attempted to elucidate the relationship between NRGs and LIHC prognosis using the LASSO regression algorithm, an R software package conducive to dimension reduction analysis and prognostic gene model construction (14). The results showed that the increased expression of neural-related genes was positively associated with the poor prognosis of LIHC (P<0.05) (Figures 5A–D). The expressions of 38 NRGs were potentially prognostic biomarkers for LIHC patients: the area under the curve (AUC) was 0.653 for 1-year, 0.646 for 3-year, and 0.665 for 5-year ROC curves (Figure 5E). Individual prognosis analysis showed that 4 of the 38 NRGs had a correlation with LIHC prognosis (P<0.05) (Figure 5F). The genes CHRNE and GFRA2, and GFRA3 and GRIN2D correlated positively and negatively with LIHC prognosis, respectively (Figure 5F). When the four aforementioned NRGs were included in the prognostic model, the AUC value was 0.672 for 1-year, 0.618 for 3-year, and 0.679 for 5-year ROC curves (Figure S2E). These results demonstrated that the expressions of specific neural-related genes including CHRNE, GFRA2, GFRA3, and GRIN2D might be potential biomarkers for LIHC prognosis.

We classified the expression of CHRNE, GFRA2, GFRA3, and GRIN2D into high and low expression groups based on the RNA-sequencing (RNA-seq) and clinical data collected from the TCGA database. A cutoff-high (top 25%) for high expression, while cutoff-low (bottom 25%) for low expression for CHRNE, GFRA2 and GRIN2D was defined. Moreover, a cutoff high (50%) and cutoff low (50%) was used as the expression threshold for GFRA3 due to its relatively low expression. The results showed that high expression of CHRNE was positively correlated with the T category and TNM staging (Figure S3A and Table 2); high expression of GFRA2 was positively correlated with the T category (Figure S3B and Table 3); high expression of GFRA3 was negatively correlated with the T category and pathological grading (Figure S3C and Table 4); high expression of GRIN2D was negatively correlated with the T category, TNM staging, and pathological grading (Figure S3D and Table 5).

The LIHC data was further classified according to the expression levels of the NRGs: CHRNE, GFRA2, GFRA3 and GRIN2D. Furthermore, GO and KEGG analyses were performed for each group. The criterion for low and high expression remains the same as mentioned above.

In LIHC with high CHRNE expression, 150 genes were unregulated and 16 genes were downregulated (FC>2, P<0.05) (Figures 6A, B), while a low expression of CHRNE had an upregulation of genes associated with physiological functions of the liver (e.g., bile secretion, cholesterol metabolism, drug metabolism) and activation processes, such as the regulation of coagulation; downregulation of tumor-promoting pathways, including the PI3K-Akt signaling pathway; and suppression of processes like transition metal ion homeostasis (Figures 6C, D).

In LIHC with a high GFRA2 expression, 195 genes were upregulated and 19 genes were downregulated (FC>2, P<0.05) (Figures S4A, B), while a low expression of GFRA2 in LIHC-activated cytokine and cytokine receptors, Th1 and Th2 cell differentiation, and Th17 cell differentiation signaling pathways, which are closely related to immune regulation, were observed. Pro-tumor pathways, such as the HIF-1 signaling pathway and p53 signaling pathway, were suppressed in the group with high GFRA2 expression. GO analysis exhibited similar results: activated immune regulation-associated processes like T-cell activation, leukocyte proliferation, regulation of T-cell activation, and inhibited cell maturation (Figures S4C, D).

In LIHC with a high GFRA3 expression, 144 genes were upregulated and 56 genes were downregulated (FC>2, P<0.05) (Figures S5A, B) For the LIHC group with a high GFRA3 expression, widely recognized critical oncogenes in liver cancer, such as AFP, IGF2, and liver cancer-associated pathways (e.g., cell cycle, forkhead box O, signaling pathway, hepatitis B, MAPK signaling pathway, VEGF signaling pathway, and p53 signaling pathway) were unregulated, while the cell proliferation-related processes (e.g., chromosome segregation, mitotic nuclear division, spindle organization) were activated. However, pathways like bile secretion and processes such as alcohol metabolism were inhibited in the high GFRA3 expression group, when compared with the low GFRA3 expression group (Figures S5C, D).

In LIHC with high GRIN2D expression, 1,971 genes were upregulated and 302 genes were downregulated (FC>2, P<0.05) (Figures S6A, B). On the other hand, the high GRIN2D expression group had an upregulation of pathways like proteoglycans in cancer, the PI3K-Akt signaling pathway, cell adhesion molecules, and activation of processes like the positive regulation of cell activation, regulation of leukocyte proliferation, downregulation of bile secretion, cholesterol metabolism, and suppression of processes like the alcohol metabolic process and lipid homeostasis (Figures S6C, D).

These findings suggest that GHRNE and GFRA2 expression in LIHC might be beneficial in maintaining the liver physiological function and suppressing tumor growth and metastasis; the effect of GFRA3 and GRIN2D was antagonistic to GHRNE and GFRA2.

The R software package, immunedeconv, was used to obtain the immune infiltration data of high/low-expression CHRNE, GFRA2, GFRA3, and GRIN2D groups of LIHC. CIBERSORT and EPIC algorithms were used herein. The CIBERSORT algorithm showed that in C1, unlike C2, high CHRNE expression was positively correlated with memory B cell (P<0.05) and mast cell (activated/resting) infiltration (P<0.05), while being negatively correlated with macrophage M0 (P<0.05) (Figure 7A). The EPIC algorithm showed that in C1, compared with C2, high CHRNE expression was positively correlated with B cell (P<0.05) (Figure 7B). In terms of high GFRA2 expression, the CIBERSORT algorithm showed that in C1, compared with C2, it was positively correlated with the CD4+ memory resting T cell (P<0.001), while it was negatively correlated with monocytes (P<0.05), macrophage M0 (P<0.001), eosinophils (P<0.05), and neutrophils (P<0.01) (Figure S7A). The EPIC algorithm showed that high GFRA2 expression was positively correlated with the CD4+ T cell (P<0.001) (Figure S7B). The CIBERSORT algorithm showed that in C1, compared with C2, high GFRA3 expression was positively correlated with the memory resting B cell (P<0.01) and T cell follicular helper (P<0.01), while it was negatively correlated with monocytes (P<0.001) (Figure S8A). The EPIC algorithm showed that low GFRA3 expression was positively correlated with macrophages (P<0.001) (Figure S8B). The CIBERSORT algorithm showed that in C1, compared with C2, high GRIN2D expression was positively correlated with Tregs (P<0.05) and macrophage M0 (P<0.001), while it was negatively correlated with naïve B cells (P<0.01), resting natural killer (NK) cells (P<0.001), monocytes (P<0.001), and activated mast cells (P<0.01) (Figure S9A). The EPIC algorithm also showed that high GRIN2D expression was positively correlated with the CD4+ T cell (P<0.001), but negatively correlated with macrophage (P<0.001) (Figure S9B).

In addition, we analyzed the correlations between ICGs and the expression of CHRNE, GFRA2, GFRA3, and GRIN2D. When compared with C2, CHRNE expression in C1 was positively correlated with SIGLEC15 (P<0.05) (Figure 7C); GFRA2 and GRIN2D expression was positively correlated with 7 of 8 IGCs including CD274, CTLA4, HAVCR2, LAG3, PDCD1, PDCD1LC2, and TIGIT with significant difference (P<0.001) (Figures S7C, S9C); GFRA3 was positively correlated with CTLA4 (P<0.001), HAVCR2 (P<0.01), LAG3 (P<0.001), PDCD1 (P<0.001), and TIGIT (P<0.01) (Figure S8C). The TIDE algorithm showed that high expression of CHRNE, GFRA3, and GRIN2D correlated with a poor immune response (Figures 7D, S7D–S9D). According to the Spearman correlation analysis of the OCLR score, the CHRNE (Figure 7E), GFRA2 (Figure S7E), and GRIN2D (Figure S9E) high-expression groups show a lower stemness score than the low-expression groups, whereas GFRA3 (Figure S8E) has an opposite result.

We obtained mutational, transcriptomic, and clinical data of LIHC patients from the TCGA database and performed visualization analysis with R software package maftools (17) and found no significant mutations for CHRNE, GFRA2, GFRA3, and GRIN2D in LIHC (Figure 8).