Clinicopathologic characteristics and survival outcomes
According to the criteria, 289 AYA patients and 257 elderly patients were finally involved in this study. The characteristics of all the patients were listed in supplementary table 2. Among the preoperative indicators, except for the hemoglobin count, platelet count, alanine aminotransferase, albumin, and AFP, there was no significant difference. In pathological indexes, only MVI had significant difference. The median OS time of the young group was 41 months, while the median survival time of the elderly group did not reach the midpoint. The 1-, 3-, and 5-year OS rates of the young group were 72.9%, 52.8%, and 44.9%, respectively, significantly lower than that in the elderly group (83.8%, 66.6%, and 62.3%, respectively, P < 0.01) (Fig. 1A). What’s more, the RFS of the two groups showed a more significant difference with 56.1%, 45.2%, and 38.0% of 1-, 3-, 5- year RFS in the young group and 73.8%, 65.5%, and 63.1% of 1-, 3-, 5- year RFS in the elderly group. And the median RFS time of AYA HCC patients was only 20 months (Fig. 1B).
All significant factors in univariate analysis and other clinical meaningful data like age, AFP and MVI were entered into the multivariate analysis, which was expressed as hazard ratio (HR) and 95% confidence interval (CI). Age < 40 years, abnormal AST, max tumor size > 5cm, positive MVI and low differentiated tumor were the main impacts of OS (Fig. 1C), while Age < 40 years, AFP ≥ 400 ng/mL, max tumor size > 5cm, positive MVI and low differentiated tumor were the major contributing factors of RFS (Fig. 1D). Similar to previous studies, we found the OS and RFS of young group were worse than elderly group, and youth was independent risk factor of both.
Metabolic Pathways Were Critical For Aya Hcc
To explore the HBV-AYA HCC specific molecular features, whole transcriptome sequencing was carried out. Firstly, differential mRNAs in tumor and adjacent tissues were sorted out in the young group and the elderly group, respectively (Fig. 2A and B). There were 4280 differential mRNAs in the young group and 3255 differential mRNAs in the elderly group, among which there were 1771 unique differential genes in the young group (Fig. 2C). Considering that the role of genes with extremely low expression is negligible to some extent, we further set the screening criteria as corrected p value < 0.05 and FPKM > 15, and finally 241 unique differentially expressed mRNAs were found in HBV-related AYA HCC patients (Fig. 2D).
KEGG functional analysis of young unique differential genes revealed many metabolism pathways, such as lipid metabolism, carbon metabolism, amino acid metabolism, PPAR signaling pathway, etc. In addition, it is also enriched in ribosomes, protein processing in endoplasmic reticulum (ER), etc. Additionally, GSEA analysis enriched similar pathways about metabolism and protein translation (Fig. 2E and F). The age-related GSEA analysis revealed that metabolism-related pathways and peroxisomes were enriched in KEGG database, which showed a clear downward trend (Fig. 2G). For Reactome database, the down-regulated genes were enriched in fatty acid metabolism and biological oxidation, while the up-regulated genes were enriched in protein translation-related pathways (Fig. 2G). These results were basically consistent with the previous enrichment analysis results, both of which inferred that HBV-AYA HCC patients had a stronger metabolic level.
Mfuzz Clustering Analysis Divided Differential Genes Into Five Clusters
In order to further analyze the functions of differential genes, and deeply understand how the expression of these genes changes with age, the trend of expression of genes according to different age groups was divided into five clusters using Mfuzz analysis. The number of genes in clusters 1 to 5 was 57, 33, 37, 40, 74, respectively (Fig. 3A). Clusters 1 and 5 were genes highly expressed in HBV-AYA HCC patient tissues but decreased in elderly group. Importantly, the expression of genes in cluster 1 showed a downward trend with age increasing in the young group. Clusters 2, 3, and 4 were genes with low expression in young group and relatively high expression in the elderly group.
The genes in clusters 1 and 5 were mainly enriched on ribosomes, and the genes in cluster 5 was also associated with phagosomes and endoplasmic reticulum. The genes in cluster 2 was mainly enriched in fatty acid metabolism and other metabolic pathways, while those in cluster 3 was more concentrated in carbon metabolism and protein-related metabolism in addition to lipid metabolism (Fig. 3B). The detail genes in cluster 2 and 3 were displayed with the heatmap (Fig. 3C).
The Co-expression Of Part Of Up-regulated Genes Was Correlated With Age
To further explore the age-related genes, WGCNA algorithm was used to perform the gene expression correlation analysis on 20 liver tumor samples. It was found that the No. 9 sample in the elderly group was obviously outlier, so it was eliminated. The remaining 19 samples were included for further analysis, and the sample correlation dendrogram and heatmap of clinical phenotypes were drawn (Fig. 4A). These differential genes were divided into 3 modules, namely blue module, brown module and turquoise module, and the co-expressed gene network heatmap of the differential genes was shown in Fig. 4B.
The WGCNA analysis defined the data of different modules as eigengene. By calculating the correlation coefficient between eigengene and clinical phenotype data, the correlation diagram between modules and phenotype was obtained. As shown in Fig. 4C, the blue module had a significant negative correlation with age (P = 0.002), which showed that with the increase of age, the expression of these genes in tumors was down-regulated. In addition, the turquoise module and age had a certain positive correlation, but no significant difference had been reached. Intriguingly, these genes had a strong positive correlation with the grade of pathological differentiation (P < 0.001), that was, the lower the gene expression, the worse the pathological grade (Fig. 4C). Functional analysis of the blue module gene showed that it was mainly concentrated in the pathway of protein processing in the ER. Moreover, the functional analysis of the turquoise group showed that these down-regulated genes were mainly related to peroxisomes, carbon metabolism, amino acid metabolism, lipid metabolism, and PPAR signaling pathway (Fig. 4D).
Low Expression Of Five Hub Metabolism Related Genes Led To Worse Prognosis Of Hcc
Although Mfuzz and WGCNA analysis had proved the AYA HCC unique genes were metabolic pathway related, it was unable to judge which genes were more critical. Sixty-four metabolism-related genes among the down-regulated differential genes were screened according to the metabolism pathways in the KEGG and Reactome databases. Combining the metabolic related products in databases such as KEGG, Reactome, Human-GEM and BRENDA, a metabolite-protein interaction network (MPIs) was constructed (Fig. 5A), in which there were 1118 protein-metabolite pairs in the network, including 531 metabolites and 63 metabolism-related genes. In addition, based on the down-regulated genes, the PPI network was constructed. The top 12 key genes were obtained by using the MCC algorithm in CytoHubba[14], including ACOX1, SCP2, EHHADH, HMGCL, ACOX2, CAT, BAAT, HAO1, LONP2, ECHS1, ACSL1 and ACAT1. The metabolism-related genes in the MPIs were sorted by degree in CytoNCA, and the top 10 genes were BAAT, EHHADH, ACSL1, CYP3A7, UGT1A1, ECHS1, ACOX1, GSTA1, MGST1 and MAOB. Taking the intersection of the two networks, five key metabolism-related genes were finally sorted out, namely EHHADH, ECHS1, ACOX1, ACSL1 and BAAT (Fig. 5A). The functions of the five genes were mainly concentrated in lipid metabolism, especially fatty acid degradation, unsaturated fatty acid biosynthesis and peroxisome.
To explore the relationship between this metabolic abnormality and the prognosis of patients with HCC, principal component (PCA) analysis of TCGA liver cancer samples based on the expression of 5 key genes was performed. All TCGA tumor samples were divided into two groups (Fig. 5B). Combined with the gene expression and clinical data of the two groups, the OS curve based on the C1 and C2 groups, as well as the difference analysis between the groups on the age distribution and the expression levels of 5 hub metabolism-related genes were drawn (Fig. 5C, D). The OS rate of C1 was significantly lower than that of C2 (P = 0.02). There were significant differences in the age distribution and the expression of five hub gene between the two groups (P < 0.01). Intriguingly, the age distribution of the C1 group was relatively younger and the expression of five key genes was lower. Therefore, it could be speculated that abnormal lipid metabolism had critical impact on the pathogenesis and prognosis of HCC, and when the age distribution was younger, the expression of lipid metabolism-related genes was relatively lower, which was line with the results obtained from our sequencing data.
Functional Analysis And Immune Infiltration Landscape Of Two Subtypes Of Hcc
Differential gene analysis was performed on the two groups of tumor samples, and a total of 1285 differentially expressed mRNAs were screened with |log2 fold change (FC)|>2 and the corrected p value < 0.01 as the threshold (Fig. 6A). GSEA analysis based on C1 and C2 was performed on all mRNAs according to the fold changes between the two groups. Compared with C2, C1 was enriched and down-regulated in metabolism-related pathways, such as fatty acid elongation, fat and carbohydrate digestion and absorption, conversely NF-κB pathway was enriched with many up-regulated genes in C1 (Fig. 6B). The PPI network was constructed for the differential genes enriched in the above pathways (Fig. 6C). Among the down-regulated pathways in C1, the fatty acid elongation pathway was at the core of these pathways and highly associated with other pathways except the fat digestion and absorption pathway, including genes such as PPT1, MECR, HADHA, HADHB, ELOVL1, ACOT7, THEM4, and THEM5. The differential genes of the NF-κB pathway were related to TNF receptors and apoptosis. And the ICAM1, PLCG1, PLCG2 and other genes in the pathway were also related to tumor migration[15].
The tumor microenvironment (TME) was consisted of a variety of nontumor cells including immune cells and stromal cells, such as endothelial cells and fibroblasts, which profoundly impacts cancer growth and invasion.[16] By ESTIMATE analysis[17], we noticed that stromal score showed no difference between the C1 and C2 groups. However, the immune score of C1 was significantly higher than that of C2 (Fig. 6D). The above analyses indicated that more pronouncedly abnormal fatty acid metabolism in HCC might have more profound impact on components of the TME. The cellular and functional TME properties can be reflected by the expression of 29 functional gene signatures (Fges) described by Bagaev et al.[18] Based on their description, it showed that the patients in C1 group were characterized by the elevated expression of Fges associated with cancer-associated fibroblast (CAF) activation, pro-tumor immune infiltration, myeloid-derived suppressor cells (MDSC), Treg, M1 signature and a higher proliferation rate. Some components of anti-tumor immune infiltrates like MHCII, coactivation molecules, B cells, T cells and checkpoint inhibition showed a significantly higher enrichment in C1, whereas no difference was observed in the remaining antitumor immune infiltrate-related terms like antitumor cytokines (Fig. 6E). We conducted the following analysis by using CIBERSORT[19], it showed that B cells memory, T cells CD4 memory activated, T cells follicular helper, Tregs, Macrophages M0, Dendritic cells resting were significantly up-regulated in C1, while T cells gamma delta, NK cells resting, Monocytes, Macrophage M2 and Mast cell resting were significantly decreased (Fig. 6F). It could be seen that there was no significant difference in CD8 + T cells and CD4 + T cells naive and B cells naive. Although there were rises and falls in different types of macrophages, M2 macrophage with anti-tumor function increased more significantly in C2, and the overall amount of macrophage and dendritic cell and M1 signature was increased in C1 with no difference in tumor associated macrophage (Fig. 6E, 6F). In generally, upregulation of anti-tumor immune infiltration was not as significant as the pro-tumor components.
The Cerna Network Construction For Hbv-aya Hcc
To further explore the molecular features of AYA HCC, we also screened young unique differential lncRNAs and miRNAs. Through preliminary screening, 1504 differential lncRNAs in the young group and 1103 differential lncRNAs in the elderly group were obtained, of which 982 differential lncRNAs were unique to the HBV-AYA patients (Fig. S1A-C). Further defined adjust p value < 0.05, a total of 445 young unique differential lncRNAs were obtained (Fig. 7A). According to the condition that |log2 fold change (FC)|> 1 and p value < 0.01, 134 differential miRNAs were screened out in the young group, and 88 differential miRNAs were screened out in the old group (Fig. S1D-F). Finally, 34 young unique differential miRNAs were screened, of which 18 were up-regulated and 16 were down-regulated in tumor tissue (Fig. 7B).
Furtherly, we calculated the correlation coefficient between differential lncRNAs and differential mRNAs to roughly predict the function of target genes that were trans-regulated by these lncRNAs. We set the threshold as the correlation coefficient > 0.85 and p value < 0.05, and 144 differential target mRNAs were obtained. It could be found that there are two major types of functions were enriched, one is the ribosome related genes, and the other is related to various metabolic pathways including fatty acid metabolism (Fig. 7C). Target mRNA of miRNA was screened by Target Scan and miRanda, and a total of 4305 differential miRNA-mRNA pairs were obtained. Functional analysis of target genes was performed respectively in up-regulated and down-regulated differential miRNA groups. The up-regulated group mainly enriched in fatty acid degradation, adipocytokine and AMPK signaling pathway (Fig. 7D). The AMPK signaling pathway is closely related to cell nutritional status, glucose metabolism and lipid metabolism.[20] Meanwhile, taking the advantage of the mirPath website and the experimental TarBase database, 34 differential miRNAs were summarized by KEGG. It was found that fatty acid metabolism, biosynthesis of unsaturated fatty acid, fatty acid biosynthesis, ECM-receptor interaction, steroid biosynthesis, and Hippo signaling pathway were enriched (Fig. 7E).
Firstly, miRNAs with a p value < 0.05 in the elderly group were excluded. Then combined with the KEGG analysis results based on the Tarbase database, the miRNAs related to lipid metabolism were retained. Finally, 22 differential miRNAs were obtained for the ceRNA network. Focusing on these miRNAs, the differential mRNAs and differential lncRNAs targeted by them were predicted under the conditions that the TargetScan score ≥ 90, the Tot Score of miRanda > 140 and Tot Energy <-20. Since mRNAs and lncRNAs often have the same expression trend in ceRNA networks, the differential lncRNA-mRNA relationship pairs with a correlation coefficient > 0.65 and a p value < 0.05 were selected, and a lncRNA-miRNA-mRNA ceRNA network was finally constructed. And the ceRNA network had a total of 28 relationship pairs (Fig. 7F).
There was a total of 10 differential miRNAs in the ceRNA network. The functional analysis of the differential mRNAs unique to youth targeted by these miRNAs revealed that these genes were enriched in fatty acid degradation as well as glycine, serine and threonine metabolism. And ACOX1 was included in the five-hub metabolism-related genes previously screened (Fig. 7G). Finally, we verified the differentially expression of CDC42SE1 and miR-378a-5p in another HBV-HCC clinical samples (Fig. 7H).