Combined DNA Methylation and Transcriptomic Assessments to Determine a Prognostic Model for PD-1-Negative Hepatocellular Carcinoma

Hepatocellular carcinoma (HCC) has the highest incidence and mortality of any malignancy in the world. Immunotherapy has been a major breakthrough for HCC treatment, but immune checkpoint inhibitors (ICIs) are effective in only a small percentage of HCC patients. In the present study, we screened programmed cell death protein 1 (PD-1) -negative HCC samples, which are frequently resistant to ICIs, and identified their methylation and transcription characteristics through the assessment of differential gene methylation and gene expression. We also screened for potential targeted therapeutic drugs using the DrugBank database. Finally, we used a LASSO (least absolute shrinkage and selection operator) regression analysis to construct a prognostic model based on three differentially methylated and expressed genes (DMEGs). The results showed that ESTIMATE (Estimation of Stromal and Immune Cells in Malignant Tumors using Expression Data) scores for the tumor samples were significantly lower compared to normal sample ESTIMATE scores. In addition, we identified 31 DMEGs that were able to distinguish PD-1-negative samples from normal samples. A functional enrichment analysis showed that these genes were involved in a variety of tumor-related pathways and immune-related pathways, and the DrugBank screening identified potential therapeutic drugs. Finally, the prognostic model based on three DMEGs (UBD, CD5L, and CD213A2) demonstrated good predictive power for HCC prognosis and was verified using an independent cohort. The present study demonstrated the methylation characteristics of PD-1-negative HCC samples, identified several potential therapeutic drugs, and proposed a prognostic model based on UBD, CD5L, and CD213A2 methylation expression. In conclusion, this work provides an in-depth understanding of methylation in HCC samples that are not sensitive to ICIs.

Introduction

Hepatocellular carcinoma (HCC) represents the malignancy with the highest incidence and fatality rate in the world (especially in East Asia and Southern Africa), resulting in approximately 800,000 deaths per year (Forner et al., 2018). Although a variety of treatment methods (e.g., radiotherapy, chemotherapy, surgical treatment, and liver transplantation) have made great progress, the prognoses of patients with HCC have not improved significantly, mainly due to continued difficulty with early diagnosis, a high rate of recurrence, and limited indicators (El-Serag et al., 2008; Kawamura et al., 2020; Kudo, 2020b; Schoenberg et al., 2020). In recent years, immunotherapy has been a major breakthrough for HCC treatment (Li et al., 2019; Zhang Q. et al., 2020). With the development and application of immune checkpoint inhibitors (ICIs) targeting programmed cell death protein 1 (PD-1) and programmed death-ligand 1, the prognoses for some HCC patients have significantly improved (Kalathil et al., 2016; Mahn et al., 2020). However, immune escape and other related mechanisms within the tumor microenvironment have not been adequately studied, and only a small proportion of patients respond to immunotherapy (Sakuishi et al., 2010; Ng et al., 2020). Therefore, it is both urgent and necessary to further explore immune mechanisms related to HCC occurrence and development to identify new targets for immunotherapy.

A well-known epigenetic modification, DNA methylation occurs mainly in mammalian CpG islands and can regulate gene transcription to ensure cell-specific programmed gene expression (Bird, 2002; Oe et al., 2021; Singh and Edwards, 2021). Many studies have shown that abnormal DNA methylation is related to the occurrence of various diseases, including cancers (Meng et al., 2018; Zhang M. et al., 2020). The latest research has also confirmed that DNA methylation is closely related to tissue immune status (Lee et al., 2001; Thomas et al., 2012; Ghoneim et al., 2017). Delacher et al. found that an important feature of differentiated regulatory T-cell populations and lymphoid T cells in different tissues was the gain (or loss) of DNA methylation (Delacher et al., 2017). In addition, the blocking of DNA methylation has also been reported to maintain the effector functions of CD8⁺ T cells during chronic infections (Ghoneim et al., 2017). Considering this close connection between DNA methylation and immune function, we reasoned that this epigenetic modification may be involved in the ICI process against HCC. However, few studies have extensively analyzed the relationship between DNA methylation and the effects of ICIs on HCC.

Previous studies have shown that the expression of PD-1 is closely related to ICI treatment (El-Khoueiry et al., 2017; Zhu et al., 2018). There have been many studies based on the expression of PD-1 in the immune microenvironment of liver cancer to explore its role, but these studies mainly focus on samples with high PD-1 expression (Tan et al., 2019; Voutsadakis, 2019). Here, we have demonstrated the DNA-methylation characteristics of PD-1-negative HCC samples and identified 31 differentially methylated and expressed genes (DMEGs) using combined analyses of methylation, transcriptomes, and prognostic information in concert with a functional-enrichment analysis to determine the potential functions of these genes. Furthermore, we identified potential HCC therapeutic drugs based on the DrugBank database. Based on these findings, we used a LASSO (least absolute shrinkage and selection operator) regression analysis to determine a prognostic model, based on three of these DMEGs, with good predictive ability. This research provides new insights for in-depth studies of methylation in PD-1-negative HCC.

Materials and Methods

Data Acquisition and Processing

Data from TCGA was obtained through the TCGA Genomic Data Commons application programming interface. We obtained the most current (October 2, 2020) TCGA-LIHC expression profile data, DNA methylation data, and clinical follow-up information. Both normal samples (n = 50) and tumor samples (n = 371) were represented in TCGA data set, and 50 normal samples and 177 tumor samples were also represented in the DNA-methylation data set. HCC samples with PD-1 expression levels lower than the average found in normal samples were regarded as PD-1-negative samples. Using this criteria, 177 PD-1-negative samples were identified.

Processing of Gene-Expression Data and DNA-Methylation Data

Gene expression data from the PD-1-negative samples was log2-converted and then analyzed for differential expression using the “limma” package in R software (Phipson et al., 2016). The p-values were converted to FDR values based on the Benjamini and Hochberg method. FDR > 0.01 and log2FC > 1 were considered to be up-regulated gene expression; FDR > 0.01 and log2FC < 1 were considered to be down-regulated gene expression.

The same R software package was applied to the DNA methylation data set from TCGA-LIHC to identify differentially methylated CpG genes (DMGs). Methylation intensities were represented by β values, and the threshold for DMG recognition was FDR < 0.05 and an absolute delta β-value > 0.3. We subsequently calculated average β-values for different regions, including the 5′- untranslated region (5′-UTR), first exon, gene body, 3′-UTR, TSS1500, and TSS200.

Immune Infiltration Analysis

The “ESTIMATE” R software package was used to determine ESTIMATE scores, stromal scores, and immune scores for both the HCC and normal samples (Verhaak et al., 2010). These scores were used to describe the overall immune-cell infiltration of the microenvironment.

Identification of DEGs and DMGs

All gene identifications based on differences in FDR < 0.01 were considered to be credible. We then performed a joint analysis of both DEGs and DMGs and divided the resulting DMEGs genes into four groups: HypoUp (β-value < -0.3 and log2FC > 1); HypoDown (β-value < -0.3 and log2FC < 1); HyperUp (β-value > -0.3 and log2FC > 1); and HyperDown (β-value > -0.3 and log2FC < 1).

Functional Enrichment Analysis

We used the “clusterProfiler” package in R software to perform a gene ontology functional enrichment analysis and a KEGG pathway-annotation analysis of the DMGs, DEGs, and DMEGs to identify the important biological processes and pathways related to these differentially expressed genes (Yu et al., 2012).

Screening for Potential Target Drugs

We screened the DrugBank database¹ to identify potential drugs capable of up-regulating DMEGs. NetworkAnalyst 3.0², a web-based visual analysis platform for analyzing and interpreting systems-level gene expression data, was used to analyze protein-drug interactions from the DMEGs based on the DrugBank database (Zhou et al., 2019). The HCC-drug proximity was calculated using the following formula:

$d(S,T)=\mathit{1}{\displaystyle \sum t}\in \mathit{Tmins}\in S(d(s,t)+\omega )$

where S represents the DMEGs; D represents the degree of BPH-related gene-set nodes in the PPIs; T represents the drug-target gene set; the distance d (s,t) represents the shortest path between the s node and the t node; and ω is the weight of the target gene. If the target gene was a gene in the BPH-related gene set, the calculation method was ω = -ln(D + 1), otherwise ω = 0.

Analysis of DMEG-Related Prognostic Signature Genes

For DMEGs, we used the Principal Component Analysis method to distinguish between HCC and adjacent samples. Linear Discriminant Analysis was used to classify the samples using DMEG expression-profile data and methylation data, and the leave-one-out cross-validation method was used for verification.

To determine relationships between DMEG expressions and prognoses, we first randomly divided the PD-1-negative samples into two groups: a training set (n = 88) and a validation set (n = 89). For the DMEG-expression and clinical-survival data, we performed 1000 LASSO regression analyses, using 10-fold cross-validation, summarized the dimensionality reduction results each time, and then counted the number of times each probe appeared per 100 times.

Results

PD-1 Expression in HCC and Normal Tissue Samples, and Their Microenvironment Characteristics

A differential-expression analysis indicated that the expression of PD-1 in HCC samples was significantly higher compared to its expression in normal samples (Figure 1A). As expected, stromal and immune scores based on the ESTIMATE (Estimation of Stromal and Immune Cells in Malignant Tumors using Expression Data) analysis were both significantly reduced in HCC samples compared to normal samples (Figure 1B). These results indicate that stromal and immune infiltrations in the tumor microenvironment were significantly inhibited, which is consistent with previous studies (Binnewies et al., 2018; Ruf et al., 2021).

FIGURE 1

Figure 1. Differences in immune infiltration between normal and HCC samples. (A) PD-1 expression differences between normal and HCC samples. (B) Differences in stromal scores, immune scores, and ESTIMATE scores between normal and HCC samples, ^∗∗∗p < 0.001.

Gene Analysis for Differential DNA Methylation

To identify differentially methylated genes, we analyzed TSS200 [transcription start site (TSS) to 200 nucleotides upstream of the TSS], TSS1500 (200 to 1500 nucleotides upstream of the TSS), and gene body methylation levels in both the PD-1-negative HCC and normal samples. The results indicated a total of 1,700 differentially methylated genes in the HCC samples. Specifically, 52 hypermethylated and 775 hypomethylated genes were identified in gene body regions, 150 hypermethylated and 407 hypomethylated genes were identified in TSS200 regions, and 90 hypermethylated and 606 hypomethylated genes were identified in TSS1500 regions (Figure 2A). The number of hypomethylated genes in these regions was far greater than the number of hypermethylated genes, especially for gene body regions (Figure 2B). A set-distribution analysis of the results showed that six hypermethylated genes were represented in all three of the regions, with 206 genes represented in only one of the regions, and that 40 hypomethylated genes were represented in all three of the regions, with 1162 genes represented in only one of the regions (Figures 2C,D), suggesting that DNA methylation levels are region-specific. The functional enrichment analysis showed that these differentially methylated genes (DMGs) were mainly enriched in 11 biological processes, four cellular components, four molecular functions, and two Kyoto Encyclopedia of Genes and Genomes (KEGG) pathways (including mRNA binding, structural constituents of the epidermis, and mRNA binding involved in posttranscriptional gene silencing) (Figure 2E).

FIGURE 2

Figure 2. Identification and characteristics of genes with differential DNA methylation in PD-1 negative HCC samples. (A) A volcano map of methylation differences of the gene body, TSS200, and TSS1500 regions in PD-1 negative HCC samples. (B) The differentially methylated genes in these three regions were mostly hypomethylated. (C) A Venn diagram of hypermethylated genes in the three regions. (D) A Venn diagram of hypomethylated genes in the three regions. (E) Gene ontology functional-enrichment analysis and KEGG annotation analysis of differentially methylated genes. The most-intense red color represents the smallest FDR, and dot size is proportional to enrichment number.

Functional Enrichment Analysis of Differentially Expressed Genes (DEGs)

Gene-expression profile data from 177 PD-1-negative tumor samples and 50 normal samples were used to determine differential gene expression. The screening criteria were false discovery rate (FDR) < 0.01 and log2FC (fold change) > 1. A total of 2249 differentially expressed genes (DEGs) were identified, of which 1659 were up-regulated in tumors and 590 were down-regulated (Figure 3A). A unsupervised hierarchical cluster analysis showed that these DEGs could distinguish between tumor and normal samples (Figure 3B). We performed a functional-enrichment analysis of these genes, and the results showed that a variety of tumor and immune-related pathways were enriched (e.g., cell cycle, DNA replication, complement and coagulation cascades, complement activation, and humoral immune responses) (Figures 3C–F).

FIGURE 3

Figure 3. Identification and characteristics of the DEGs based on PD-1-negative HCC samples. (A) A volcano map of the DEGs. (B) A differential-expression heat map of the normal and HCC samples. (C) KEGG annotation analysis of the DEGs. (D) Gene ontology biological-function enrichment analysis of the DEGs. (E) DEG analysis for gene ontology cell-component enrichment. (F) Molecular function enrichment analysis of the DEGs.

Identification of DMEGs

Given the importance of both methylation and transcription in HCC occurrence and development, we jointly analyzed DMGs and DEGs to more fully explore the relationship between these two processes with the idea that genes demonstrating both differential DNA methylation and expression may play crucial HCC roles (Hu et al., 2020). A set-distribution analysis showed that 11 DMEGs were identified in gene body regions, 18 DMEGs were identified in TSS200 regions, and 15 DMEGs were identified in TSS1500 regions (Figures 4A–C). In addition, the extent to which these genes showed both differential methylation and expression is shown in Figure 4D. The observation that most of these DMEGs were hypomethylated in tumors is consistent with previous reports (Fabianowska-Majewska et al., 2021). Interestingly, several DMEGs (e.g., TBX15, REG1A, and HBB) were found in all three regions, indicating that these genes may have transcriptional differences caused by differential methylation. Previous reports have shown that TBX15 expression can be used as a prognostic marker for HCC (Morine et al., 2020), and HBB has been reported to play a key role in prostate cancer differentiation and in a variety of important biological pathways (e.g., iron metabolism) (Chen and Sun, 2021; Lin et al., 2021). Based on the present expression data, we identified four different regulatory gene sets in the three DNA regions, for a total of 31 DMEGs (Figure 4E).

FIGURE 4

Figure 4. Identification and expression patterns of DMEGs. (A) A Venn diagram of the DEGs and the DMGs in gene body regions. (B) A Venn diagram of the DEGs and the DMGs in TSS200 regions. (C) A Venn diagram of the DEGs and the DMGs in TSS1500 regions. (D) Distribution of the four DMEG expression modes in the three regions. (E) A histogram of the four DMEG regulatory modes in the three regions.

To better understand the potential functions of these genes in HCC, an in-depth chromosome (chr) distribution analysis was carried out (Figure 5A). We found that chr7, chr18, and chrX had the most DMEGs. Interestingly, the methylation patterns of the DMEGs in adjacent chromosomal gene regions were roughly consistent, suggesting that these genes may have cooperative expression patterns and functions. In order to verify gene effectiveness, we used the expression profiles of these 31 DMEGs and the methylation data from the DMEGs in each of the three regions to construct a linear judgment classification model. The results showed that these DMEGs could effectively distinguish PD-1-negative HCC samples from normal samples (Figure 5B). The corresponding receiver operating characteristic (ROC) curve analysis showed an area-under-the-curve (AUC) value ≥0.99 (Figure 5C). The functional enrichment analysis showed that DMEGs were mainly enriched in amyloid and apoptotic-cell clearance, as well as in the activation of signaling pathways, such as mitogen-activated protein kinase pathways (Figure 5D).

FIGURE 5

Figure 5. Characteristics of DMEGs and drug distance. (A) DMEG distributions on chromosomes. (B) Principal component analysis of DMEG expression and methylation in the three DNA regions. (C) Construction of the linear discriminant model based on the expression profiles and methylation levels of DMEGs for predicting the ROC curves for HCC and normal samples. (D) The results of KEGG pathway-annotation and gene ontology enrichment analyses for the DMEGs, with different colors representing different pathways, and the lines indicating relationships between genes and pathways. (E) A drug distance-density plot for the differentially methylated and expressed gene set.

Identification of Potential Target Drugs Based on PD-1-Negative DMEGs

A DMEG protein-drug interaction analysis was performed using the NetworkAnalyst 3.0 tool based on the DrugBank database, and 6 genes were found to interact with drugs (Table 1). Among them, cystic fibrosis transmembrane conductance regulator (CFTR) and hemoglobin subunit beta (HBB) proteins had the most drug interactions: 12 with HBB (e.g., pentaerythritol tetranitrate, and 4-carboxycinnamic acid) and 7 with CFTR (e.g., colforsin, crofelemer, and lonidamine). Some of these identified drugs may be effective against HCC.

TABLE 1

Table 1. The interaction between DMEGs and drugs.

Based on these drug-target pairs from DrugBank, as well as the string key protein-protein interaction (PPI) network (threshold score was set at 600), we calculated drug-HCC proximities. For both our DMEGs as samples, and for randomly selected genes as samples, we found that the number of drugs was significantly reduced when drug distance was less than 0.8 (Figure 5E). This suggests that when HCC-drug proximity is less than 0.8, the drug may have a targeted impact on the disease.

Molecular Docking Analysis Verifies the Affinity of Candidate Drugs

Considering the accuracy of molecular docking, we chose SRD5A2 with a moderate molecular weight as a representative to perform molecular docking analysis in order to clarify the binding model between drug candidates and gene targets. We first downloaded the 3D model of SRD5A2 protein (PDB ID: 7BW1) from the PDB database for molecular docking experiments (Figure 6A). Autodock Vina molecular docking results show that the compound can bind tightly to the active site of the SRD5A2 protein, with a molecular docking score of -6.5kcal/mol (Figure 6B). In addition, we found that the compound could generate favorable hydrogen bonds with the important amino acid residues GLU57, GLN56 and TYR91 in the SRD5A2 protein, as shown in Figure 6C. The above results suggest that the drug Azelaic acid can interact closely with the SRD5A2 protein, thereby affecting the activity of the SRD5A2 protein. Meanwhile, we used molecular dynamics simulation to further evaluate the stability of the protein model combined with the drug, and used the RMSD method to estimate the stability of the protein model (Figure 6D). During the 100ns molecular dynamics simulation, we can find that the protein-backbone is maintained in a relatively stable state as a whole, indicating that the protein is relatively stable during the molecular dynamics simulation.

FIGURE 6

Figure 6. SRD5A2- Azelaic acid molecular docking analysis. (A) 3D model of SRD5A2 protein structure. (B) 2D interaction diagram of Azelaic acid and SRD5A2 protein. (C) The binding domain of Azelaic acid and SRD5A2 protein, in which hydrogen bonds are shown in green and hydrophobic interactions are shown in pink. Amino acid residues are displayed as steel blue. (D) The curve of the RMSD value of the protein backbone during molecular dynamics in 100ns.

Prognostic Genetic Signature of DMEGs in PD-1-Negative Samples

We used a LASSO regression analysis to reduce the dimensionality of the expression and prognostic data for these DMEGs, and obtained a combined maximum frequency for three probe genes (Figure 7A; ENSG00000073754, ENSG00000123496, and ENSG00000213886). The trajectories for these three genes with different lambdas are shown in Figure 7B, and the standard-deviation distribution of the different lambdas is shown in Figure 7C. The survival-model results demonstrated that, with a median cutoff, the high-expression group was significantly different from the low-expression group using these three genes, indicating a highly effective model (Figures 7D–F). According to the LASSO analysis, the determination formula was:

$\mathit{RiskScore}=-\mathit{0.253}\times \mathit{ENSG00000073754}+\mathit{0.111}\times \mathit{ENSG00000213886}-\mathit{0.843}\times \mathit{ENSG00000123496}$

FIGURE 7

Figure 7. Prognostic gene signatures related to DMEGs. (A) LASSO-regression frequencies for each combination of genes. (B) Change-coefficient trajectories for each of the three probe genes using different lambdas. (C) The distribution of standard deviations in this model under different lambda conditions. (D) Kaplan-Meier curves for the CD5L high- and low-expression groups. (E) Kaplan-Meier curves for the CD213A2 high- and low-expression groups. (F) Kaplan-Meier curves for the UBD high- and low-expression groups.

Detailed information about these three genes is presented in Table 2. Both CD5L and CD213A2 were determined to be HCC protective factors, and UBD was determined to be a risk factor.

TABLE 2

Table 2. LASSO identifies 3 prognostic-related DMEGs.

According to the above formula, we calculated the risk score for each sample and determined corresponding survival status and expression changes for the three different signature genes as risk values increased (Figure 8A). We found that most of the training-set samples had higher risk scores, and that samples with high-risk scores had worse prognoses. In addition, UBD expression was found to increase with increasing risk-score values, while CD5L and CD213A2 expressions decreased with increasing risk scores. Furthermore, we conducted ROC analyses for the prognostic classification of risk scores, and the prognostic prediction classification efficiencies for 1 year, 3 years, and 5 years were 0.64, 0.74, and 0.94, respectively (Figure 8A). These data indicate that this model was highly predictive for long-term survival. In addition, based on risk-score value z-scores, we divided the samples into high- and low-risk groups. The survival curves showed very significant differences between these groups (log rank P = 0.0011), in which 57 of the samples were classified as high-risk and 19 samples were classified as low-risk. We also conducted an analysis of the validation set using the same model and coefficients as for the training-set analysis. The validation results showed that the expression trends of these three signature genes were consistent with the training set (Figure 8B). Compared to the low-risk group, the overall survival rate in the high-risk group was worse, but this difference was not statistically significant.

FIGURE 8

Figure 8. Characteristics of prognostic gene signatures based on 3 DMEGs. (A) Training-set relationships between risk scores, survival status, and gene expressions for the three DMEGs and prognostic-signature ROC analysis and Kaplan-Meier curves. (B) Validation-set relationships between risk scores, survival status, and gene expressions for the three DMEGs and prognostic-signature ROC analysis and Kaplan-Meier curves. (C) The total TCGA-LIHC-set PD-1 negative samples relationships between risk scores, survival status, and gene expressions of the three DMEGs and prognostic-signature ROC analysis and Kaplan-Meier curves. (D) The GSE10141-set PD-1 negative samples relationships between risk scores, survival status, and gene expressions of the three DMEGs and prognostic-signature ROC analysis and Kaplan-Meier curves.

In order to further verify the effectiveness of this gene-signature model, we applied it to all PD-1 low-expression samples in the Cancer Genome Atlas (TCGA) and GSE10141 HCC data set. Similar to the previous results, most of these samples had higher risk scores that related to lower survival status (Figure 8C). The ROC analysis showed that the AUC values at 1, 3, and 5 years were 0.63, 0.65, and 0.82, respectively (Figure 8C). At the same time, the overall survival rates for the high-risk group and the low-risk group were also significantly different (P = 0.0016, Figure 8C). The GSE10141 data set results showed that more samples were identified as having lower-risk scores, and that the 1-year, 3-year, and 5-year AUC values were 0.78, 0.78, and 0.67, respectively (Figure 8D). In addition, the prognoses for the high-risk group and the low-risk group were significantly different (P = 0.00099). These differences in prognostic scores and expression levels may have been due to the batch effect on the different platforms.

Discussion

With the development of ICIs, immunotherapy has become the new focus of attention in the field of tumor treatments (Bersanelli et al., 2021). However, only a small proportion of HCC patients actually respond to ICIs, and one of the important reasons for this may be HCC expression of PD-1 and the infiltration of CD4⁺ and CD8⁺ T cells (Sia et al., 2017; Kurebayashi et al., 2018; Kudo, 2020a). Previous studies have shown that DNA methylation can affect the immune status of the tumor microenvironment and tumor responses to ICIs, and that a lack of DNA methylation is related to immune-evasion characteristics (Duruisseaux et al., 2018; Jung et al., 2019). Therefore, it is of instructive significance to investigate the changes in DNA methylation in PD1-negative HCC samples on the loss of the anti-tumor effect of ICIs and to further search for other therapeutic targets.

Here, using the Liver Hepatocellular Carcinoma (LIHC) cohort from TGCA, we found that stromal-score, immune-score, and ESTIMATE-score values for tumor samples were significantly lower compared to normal samples, and that PD-1 expression was higher compared to normal samples. These differences highlight the immunosuppressive state in these tumors. Using PD-1-negative samples for screening, we identified DEGs and DMGs in three DNA regions, and a functional enrichment analysis showed that they were related to tumor immunity and cancer-related cell pathways (Hu et al., 2020). As there is evidence that DMEGs play a key role in tumorigenesis, we conducted a joint DEG-DMG analysis and identified 31 gene candidates. We then divided these DMEGs into four groups based on modes of expression. Interestingly, most DMEGs demonstrated hypomethylation. Consistent with these results, studies have shown that overall demethylation and hypomethylation of both oncogenes and metastasis-promoting genes are characteristics common to almost all cancers, including HCC (Zhao et al., 2020; Fabianowska-Majewska et al., 2021). The chromosomal distribution of these genes showed that DMEGs in adjacent gene regions had similar expression patterns, suggesting that adjacent regions may be regulated by the same methylases/demethylases. In addition, these 31 DMEGs could distinguish between tumor and normal samples based on methylation and expression levels, indicating their potential importance.

The protein-drug interaction analysis provided another perspective for evaluating the potential therapeutic effects of DMEGs on HCC. Conductance regulator and HBB were identified as the DMEGs having the most drug interactions. Conductance regulator encodes chloride ion and bicarbonate ion channels and has been implicated in a variety of cancers. It has also been identified as a molecular biomarker for early HCC diagnosis (Hogan, 1999; Moribe et al., 2009). Hemoglobin subunit beta has also been reported to be a diagnostic biomarker in cancers (Shi et al., 2018). Here, we found that these genes interacted with a variety of drugs, including the commonly used anti-tumor drug Lonidamine, a mitochondrial hexokinase inhibitor, which can inhibit the glycolysis of tumor cells. The drugs identified using the analysis above may provide alternative ways to treat HCC. In addition, we found that when the drug-DMEG distance was less than 0.8, drug interactions were significantly reduced, suggesting that 0.8 may represent an important threshold. With distances <0.8, the corresponding drugs may have more precise targeting effects for HCC treatment.

In order to evaluate the predictive power of these DMEGs for HCC prognosis, we used a LASSO regression to determine a prognostic model for HCC based on three genes (UBD, CD5L, and CD213A2) that can modulate immune responses. Ubiquitin D (UBD), a ubiquitin-like protein modifier, binds to target proteins by covalent bonding and then causes them to be degraded by the 26s proteasome (Guarascio et al., 2020). Studies have shown that UBD can regulate the activation of the tumor necrosis factor alpha (TNF-α) -induced, and lipopolysaccharide-mediated, innate immune response central mediator nuclear factor kappa B (NF-κB), by promoting TNF-α-mediated ubiquitinated-I-κB-α proteasome degradation (Kawamoto et al., 2019). We speculate that increased UBD expression may promote an immune response in the tumor microenvironment, which suppresses tumor growth. The cysteine-rich inflammatory regulator CD5L has been shown to promote proliferation and activate autophagy in HCC by binding heat shock-A5 proteins (Armengol et al., 2013; Sanjurjo et al., 2015; Aran et al., 2018). Consistent with previous studies, the present results show that CD5L expression is a risk factor for HCC and may help researchers to reinterpret its role from the new perspective of methylation and immunity. CD213A2 has been shown to bind to interleukin-13 (IL-13) and activate its immunomodulatory function (Tabata and Khurana Hershey, 2007); however, no direct tumor-related role has been found, so further study of CD213A2 and its effects on HCC, including through methylation and immunity, is warranted.

The tumor’s response to ICIs largely depends on the state of the tumor microenvironment. As mentioned above, the lack of DNA methylation is related to the immune evasion characteristics of the tumor microenvironment. Here, we explored the DNA methylation characteristics and potential functional pathways of PD1-negative HCC patients, and identified the genes that play a key role in this process. These genes participate in the tumor immune microenvironment through possible DNA methylation regulation and further affect the anti-tumor effect of ICIs. We have proposed a model for determining PD-1-negative HCC prognoses based on these three genes. This model was clearly able to divide PD-1-negative HCC samples into high- and low-risk groups, with clear trends for DMEG expressions, and significantly different prognoses between these two groups. The use of an independent verification queue also confirmed its effectiveness. Therefore, further exploring the role of these hub genes in this process will help guide researchers to have a deeper understanding of the PD1-negative tumor microenvironment. These hub genes are also expected to become potential targets for enhancing the efficacy of ICIs.

Conclusion

The present research has revealed the methylation/transcription characteristics of PD-1-negative HCC samples and identified potential therapeutic targets and drugs. Most importantly, we have demonstrated the effectiveness of a prognostic model for HCC based on three DMEGs. These results provide insights into potential treatment strategies for HCC that are not sensitive to PD-1 inhibitors and into mechanisms by which methylation may affect HCC.

Data Availability Statement

Publicly available datasets were analyzed in this study. The names of the repository/repositories and accession number(s) can be found in the article/supplementary material.

Author Contributions

WG designed the study, wrote this manuscript, and collected samples. LZ searched the articles and made figures. Both authors read and approved the final manuscript.

Funding

This study was supported by funds from the National Natural Science Foundation of China (81671958). Leading Talents of Zhongyuan Science and Technology Innovation (214200510027), Henan Provincial Medical Science and Technology Research Plan (SBGJ2018002), Science and Technology Innovation Talents in Henan Universities (19HASTIT003), and Outstanding Foreign Scientist Studio in Henan Province (GZS2020004).

Conflict of Interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Publisher’s Note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

Acknowledgments

We thank the patients and investigators who participated in TCGA and GSE10141 for providing data.

Abbreviations

HCC, hepatocellular carcinoma; ICIs, immune checkpoint inhibitors; PD-1, programmed cell death protein 1; LASSO, least absolute shrinkage and selection operator; DMEGs, differentially methylated and expressed genes; ESTIMATE, Estimation of Stromal and Immune Cells in Malignant Tumors using Expression Data; DMGs, differentially methylated genes; KEGG, Kyoto Encyclopedia of Genes and Genomes; DEGs, differentially expressed genes; TCGA, Cancer Genome Atlas; LIHC, liver hepatocellular carcinoma; UBD, Ubiquitin D; TNF-α, tumor necrosis factor alpha; NF-κB, nuclear factor kappa B; IL-13, interleukin-13.

Footnotes

1. ^ https://www.drugbank.com/
2. ^ http://www.networkanalyst.ca/

References

Aran, G., Sanjurjo, L., Barcena, C., Simon-Coma, M., Tellez, E., Vazquez-Vitali, M., et al. (2018). CD5L is upregulated in hepatocellular carcinoma and promotes liver cancer cell proliferation and antiapoptotic responses by binding to HSPA5 (GRP78). FASEB J. 32, 3878–3891. doi: 10.1096/fj.201700941RR

PubMed Abstract | CrossRef Full Text | Google Scholar

Armengol, C., Bartoli, R., Sanjurjo, L., Serra, I., Amezaga, N., Sala, M., et al. (2013). Role of scavenger receptors in the pathophysiology of chronic liver diseases. Crit. Rev. Immunol. 33, 57–96.

Google Scholar

Bersanelli, M., Buti, S., Giannatempo, P., Raggi, D., Necchi, A., Leonetti, A., et al. (2021). Outcome of patients with advanced upper tract urothelial carcinoma treated with immune checkpoint inhibitors: a systematic review and meta-analysis. Crit. Rev. Oncol. Hematol. 159:103241. doi: 10.1016/j.critrevonc.2021.103241

PubMed Abstract | CrossRef Full Text | Google Scholar

Binnewies, M., Roberts, E. W., Kersten, K., Chan, V., Fearon, D. F., Merad, M., et al. (2018). Understanding the tumor immune microenvironment (TIME) for effective therapy. Nat. Med. 24, 541–550. doi: 10.1038/s41591-018-0014-x

PubMed Abstract | CrossRef Full Text | Google Scholar

Bird, A. (2002). DNA methylation patterns and epigenetic memory. Genes Dev. 16, 6–21. doi: 10.1101/gad.947102

PubMed Abstract | CrossRef Full Text | Google Scholar

Chen, D., and Sun, X. (2021). [Research progress of iron metabolism in phenotype modification of beta-thalassemia]. Zhonghua Yi Xue Yi Chuan Xue Za Zhi 38, 27–31. doi: 10.3760/cma.j.cn511374-20190819-00418

PubMed Abstract | CrossRef Full Text | Google Scholar

Delacher, M., Imbusch, C. D., Weichenhan, D., Breiling, A., Hotz-Wagenblatt, A., Trager, U., et al. (2017). Genome-wide DNA-methylation landscape defines specialization of regulatory T cells in tissues. Nat. Immunol. 18, 1160–1172. doi: 10.1038/ni.3799

PubMed Abstract | CrossRef Full Text | Google Scholar

Duruisseaux, M., Martinez-Cardus, A., Calleja-Cervantes, M. E., Moran, S., Castro de Moura, M., Davalos, V., et al. (2018). Epigenetic prediction of response to anti-PD-1 treatment in non-small-cell lung cancer: a multicentre, retrospective analysis. Lancet Respir. Med. 6, 771–781. doi: 10.1016/S2213-2600(18)30284-4

CrossRef Full Text | Google Scholar

El-Khoueiry, A. B., Sangro, B., Yau, T., Crocenzi, T. S., Kudo, M., Hsu, C., et al. (2017). Nivolumab in patients with advanced hepatocellular carcinoma (CheckMate 040): an open-label, non-comparative, phase 1/2 dose escalation and expansion trial. Lancet 389, 2492–2502. doi: 10.1016/S0140-6736(17)31046-2

CrossRef Full Text | Google Scholar

El-Serag, H. B., Marrero, J. A., Rudolph, L., and Reddy, K. R. (2008). Diagnosis and treatment of hepatocellular carcinoma. Gastroenterology 134, 1752–1763. doi: 10.1053/j.gastro.2008.02.090

PubMed Abstract | CrossRef Full Text | Google Scholar

Fabianowska-Majewska, K., Kaufman-Szymczyk, A., Szymanska-Kolba, A., Jakubik, J., Majewski, G., and Lubecka, K. (2021). Curcumin from turmeric rhizome: a potential modulator of DNA methylation machinery in breast cancer inhibition. Nutrients 13:332. doi: 10.3390/nu13020332

PubMed Abstract | CrossRef Full Text | Google Scholar

Forner, A., Reig, M., and Bruix, J. (2018). Hepatocellular carcinoma. Lancet 391, 1301–1314. doi: 10.1016/S0140-6736(18)30010-2

CrossRef Full Text | Google Scholar

Ghoneim, H. E., Fan, Y., Moustaki, A., Abdelsamed, H. A., Dash, P., Dogra, P., et al. (2017). De Novo Epigenetic programs inhibit PD-1 blockade-mediated T cell rejuvenation. Cell 170, 142–157e19. doi: 10.1016/j.cell.2017.06.007 142-157 e19

PubMed Abstract | CrossRef Full Text | Google Scholar

Guarascio, R., Salih, D., Yasvoina, M., Edwards, F. A., Cheetham, M. E., and van der Spuy, J. (2020). Negative regulator of ubiquitin-like protein 1 modulates the autophagy-lysosomal pathway via p62 to facilitate the extracellular release of tau following proteasome impairment. Hum. Mol. Genet. 29, 80–96. doi: 10.1093/hmg/ddz255

PubMed Abstract | CrossRef Full Text | Google Scholar

Hogan, D. B. (1999). Did Osler suffer from “paranoia antitherapeuticum baltimorensis”? A comparative content analysis of the principles and practice of medicine and Harrison’s principles of internal medicine, 11th edition. CMAJ 161, 842–845.

Google Scholar

Hu, W., Wang, G., Yarmus, L. B., and Wan, Y. (2020). Combined methylome and transcriptome analyses reveals potential therapeutic targets for EGFR wild type lung cancers with low PD-L1 expression. Cancers (Basel) 12:2496. doi: 10.3390/cancers12092496

PubMed Abstract | CrossRef Full Text | Google Scholar

Jung, H., Kim, H. S., Kim, J. Y., Sun, J. M., Ahn, J. S., Ahn, M. J., et al. (2019). DNA methylation loss promotes immune evasion of tumours with high mutation and copy number load. Nat. Commun. 10:4278. doi: 10.1038/s41467-019-12159-9

PubMed Abstract | CrossRef Full Text | Google Scholar

Kalathil, S. G., Lugade, A. A., Miller, A., Iyer, R., and Thanavala, Y. (2016). PD-1(+) and Foxp3(+) T cell reduction correlates with survival of HCC patients after sorafenib therapy. JCI Insight 1:e86182. doi: 10.1172/jci.insight.86182

PubMed Abstract | CrossRef Full Text | Google Scholar

Kawamoto, A., Nagata, S., Anzai, S., Takahashi, J., Kawai, M., Hama, M., et al. (2019). Ubiquitin D is upregulated by synergy of notch signalling and TNF-alpha in the inflamed intestinal epithelia of IBD Patients. J. Crohns Colitis 13, 495–509. doi: 10.1093/ecco-jcc/jjy180

PubMed Abstract | CrossRef Full Text | Google Scholar

Kawamura, Y., Kobayashi, M., Shindoh, J., Kobayashi, Y., Okubo, S., Tominaga, L., et al. (2020). Lenvatinib-transarterial chemoembolization sequential therapy as an effective treatment at Progression during lenvatinib therapy for advanced hepatocellular carcinoma. Liver Cancer 9, 756–770. doi: 10.1159/000510299

PubMed Abstract | CrossRef Full Text | Google Scholar

Kudo, M. (2020a). Gd-EOB-DTPA-MRI could predict WNT/beta-catenin mutation and resistance to immune checkpoint inhibitor therapy in hepatocellular carcinoma. Liver Cancer 9, 479–490. doi: 10.1159/000509554

PubMed Abstract | CrossRef Full Text | Google Scholar

Kudo, M. (2020b). Recent advances in systemic therapy for hepatocellular carcinoma in an aging society: 2020 update. Liver Cancer 9, 640–662. doi: 10.1159/000511001

PubMed Abstract | CrossRef Full Text | Google Scholar

Kurebayashi, Y., Ojima, H., Tsujikawa, H., Kubota, N., Maehara, J., Abe, Y., et al. (2018). Landscape of immune microenvironment in hepatocellular carcinoma and its additional impact on histological and molecular classification. Hepatology 68, 1025–1041. doi: 10.1002/hep.29904

PubMed Abstract | CrossRef Full Text | Google Scholar

Lee, P. P., Fitzpatrick, D. R., Beard, C., Jessup, H. K., Lehar, S., Makar, K. W., et al. (2001). A critical role for Dnmt1 and DNA methylation in T cell development, function, and survival. Immunity 15, 763–774. doi: 10.1016/s1074-7613(01)00227-8

CrossRef Full Text | Google Scholar

Li, W., Wang, H., Ma, Z., Zhang, J., Ou-Yang, W., Qi, Y., et al. (2019). Multi-omics analysis of microenvironment characteristics and immune escape mechanisms of hepatocellular carcinoma. Front. Oncol. 9:1019. doi: 10.3389/fonc.2019.01019

PubMed Abstract | CrossRef Full Text | Google Scholar

Lin, X. D., Lin, N., Lin, T. T., Wu, Y. P., Huang, P., Ke, Z. B., et al. (2021). Identification of marker genes and cell subtypes in castration-resistant prostate cancer cells. J. Cancer 12, 1249–1257. doi: 10.7150/jca.49409

PubMed Abstract | CrossRef Full Text | Google Scholar

Mahn, R., Vogt, A., Kupczyk, P., Sadeghlar, F., van Beekum, K., Huneburg, R., et al. (2020). Programmed cell death protein 1 (PD-1)-inhibition in hepatocellular carcinoma (HCC): a single center experience. Scand J. Gastroenterol. 55, 1057–1062. doi: 10.1080/00365521.2020.1794539

PubMed Abstract | CrossRef Full Text | Google Scholar

Meng, C., Shen, X., and Jiang, W. (2018). Potential biomarkers of HCC based on gene expression and DNA methylation profiles. Oncol. Lett. 16, 3183–3192. doi: 10.3892/ol.2018.9020

PubMed Abstract | CrossRef Full Text | Google Scholar

Moribe, T., Iizuka, N., Miura, T., Kimura, N., Tamatsukuri, S., Ishitsuka, H., et al. (2009). Methylation of multiple genes as molecular markers for diagnosis of a small, well-differentiated hepatocellular carcinoma. Int. J. Cancer 125, 388–397. doi: 10.1002/ijc.24394

PubMed Abstract | CrossRef Full Text | Google Scholar

Morine, Y., Utsunomiya, T., Saito, Y., Yamada, S., Imura, S., Ikemoto, T., et al. (2020). Reduction of T-Box 15 gene expression in tumor tissue is a prognostic biomarker for patients with hepatocellular carcinoma. Oncotarget 11, 4803–4812. doi: 10.18632/oncotarget.27852

PubMed Abstract | CrossRef Full Text | Google Scholar

Ng, H. H. M., Lee, R. Y., Goh, S., Tay, I. S. Y., Lim, X., Lee, B., et al. (2020). Immunohistochemical scoring of CD38 in the tumor microenvironment predicts responsiveness to anti-PD-1/PD-L1 immunotherapy in hepatocellular carcinoma. J. Immunother. Cancer 8:e000987. doi: 10.1136/jitc-2020-000987

PubMed Abstract | CrossRef Full Text | Google Scholar

Oe, M., Ojima, K., and Muroya, S. (2021). Difference in potential DNA methylation impact on gene expression between fast- and slow-type myofibers. Physiol. Genomics 53, 69–83. doi: 10.1152/physiolgenomics.00099.2020

PubMed Abstract | CrossRef Full Text | Google Scholar

Phipson, B., Lee, S., Majewski, I. J., Alexander, W. S., and Smyth, G. K. (2016). Robust hyperparameter estimation protects against hypervariable genes and improves power to detect differential expression. Ann. Appl. Stat. 10, 946–963. doi: 10.1214/16-AOAS920

CrossRef Full Text | Google Scholar

Ruf, B., Heinrich, B., and Greten, T. F. (2021). Immunobiology and immunotherapy of HCC: spotlight on innate and innate-like immune cells. Cell Mol. Immunol. 18, 112–127. doi: 10.1038/s41423-020-00572-w

PubMed Abstract | CrossRef Full Text | Google Scholar

Sakuishi, K., Apetoh, L., Sullivan, J. M., Blazar, B. R., Kuchroo, V. K., and Anderson, A. C. (2010). Targeting Tim-3 and PD-1 pathways to reverse T cell exhaustion and restore anti-tumor immunity. J. Exp. Med. 207, 2187–2194. doi: 10.1084/jem.20100643

PubMed Abstract | CrossRef Full Text | Google Scholar

Sanjurjo, L., Aran, G., Roher, N., Valledor, A. F., and Sarrias, M. R. (2015). AIM/CD5L: a key protein in the control of immune homeostasis and inflammatory disease. J. Leukoc. Biol. 98, 173–184. doi: 10.1189/jlb.3RU0215-074R

PubMed Abstract | CrossRef Full Text | Google Scholar

Schoenberg, M. B., Ehmer, U., Umgelter, A., Bucher, J. N., Koch, D. T., Borner, N., et al. (2020). Liver transplantation versus watchful-waiting in hepatocellular carcinoma patients with complete response to bridging-therapy. Transpl. Int. 34, 465–473. doi: 10.1111/tri.13808

PubMed Abstract | CrossRef Full Text | Google Scholar

Shi, F., Wu, H., Qu, K., Sun, Q., Li, F., Shi, C., et al. (2018). Identification of serum proteins AHSG, FGA and APOA-I as diagnostic biomarkers for gastric cancer. Clin. Proteomics 15:18. doi: 10.1186/s12014-018-9194-0

PubMed Abstract | CrossRef Full Text | Google Scholar

Sia, D., Jiao, Y., Martinez-Quetglas, I., Kuchuk, O., Villacorta-Martin, C., Castro de Moura, M., et al. (2017). Identification of an immune-specific class of hepatocellular carcinoma. based on molecular features. Gastroenterology 153, 812–826. doi: 10.1053/j.gastro.2017.06.007

PubMed Abstract | CrossRef Full Text | Google Scholar

Singh, M. K., and Edwards, J. R. (2021). Integrated analysis of DNA methylation, hydroxymethylation, and gene expression data USING ME-Class2. Methods Mol. Biol. 2198, 467–489. doi: 10.1007/978-1-0716-0876-0_34

CrossRef Full Text | Google Scholar

Tabata, Y., and Khurana Hershey, G. K. (2007). IL-13 receptor isoforms: breaking through the complexity. Curr. Allergy Asthma Rep. 7, 338–345. doi: 10.1007/s11882-007-0051-x

PubMed Abstract | CrossRef Full Text | Google Scholar

Tan, K. W., Chacko, A. M., and Chew, V. (2019). PD-1 expression and its significance in tumour microenvironment of hepatocellular carcinoma. Transl. Gastroenterol. Hepatol. 4:51. doi: 10.21037/tgh.2019.06.08

PubMed Abstract | CrossRef Full Text | Google Scholar

Thomas, R. M., Gamper, C. J., Ladle, B. H., Powell, J. D., and Wells, A. D. (2012). De novo DNA methylation is required to restrict T helper lineage plasticity. J. Biol. Chem. 287, 22900–22909. doi: 10.1074/jbc.M111.312785

PubMed Abstract | CrossRef Full Text | Google Scholar

Verhaak, R. G., Hoadley, K. A., Purdom, E., Wang, V., Qi, Y., Wilkerson, M. D., et al. (2010). Integrated genomic analysis identifies clinically relevant subtypes of glioblastoma characterized by abnormalities in PDGFRA. IDH1, EGFR, and NF1. Cancer Cell 17, 98–110. doi: 10.1016/j.ccr.2009.12.020

PubMed Abstract | CrossRef Full Text | Google Scholar

Voutsadakis, I. A. (2019). PD-1 inhibitors monotherapy in hepatocellular carcinoma: meta-analysis and systematic review. Hepatobiliary Pancreat. Dis. Int. 18, 505–510. doi: 10.1016/j.hbpd.2019.09.007

PubMed Abstract | CrossRef Full Text | Google Scholar

Yu, G., Wang, L. G., Han, Y., and He, Q. Y. (2012). clusterProfiler: an R package for comparing biological themes among gene clusters. OMICS 16, 284–287. doi: 10.1089/omi.2011.0118

PubMed Abstract | CrossRef Full Text | Google Scholar

Zhang, M., Wang, Y., Wang, Y., Jiang, L., Li, X., Gao, H., et al. (2020). Integrative analysis of DNA methylation and gene expression to determine specific diagnostic biomarkers and prognostic biomarkers of breast cancer. Front. Cell Dev. Biol. 8:529386. doi: 10.3389/fcell.2020.529386

PubMed Abstract | CrossRef Full Text | Google Scholar

Zhang, Q., Yu, X., Zheng, Q., He, Y., and Guo, W. (2020). A molecular subtype model for liver HBV-related hepatocellular carcinoma patients based on immune-related genes. Front. Oncol. 10:560229. doi: 10.3389/fonc.2020.560229

PubMed Abstract | CrossRef Full Text | Google Scholar

Zhao, Y., Li, K., Sun, J., He, N., Zhao, P., Zang, C., et al. (2020). Genomic DNA methylation profiling indicates immune response following thermal ablation treatment for HBV-associated hepatocellular carcinoma. Oncol. Lett. 20, 677–684. doi: 10.3892/ol.2020.11636

PubMed Abstract | CrossRef Full Text | Google Scholar

Zhou, G., Soufan, O., Ewald, J., Hancock, R. E. W., Basu, N., and Xia, J. (2019). NetworkAnalyst 3.0: a visual analytics platform for comprehensive gene expression profiling and meta-analysis. Nucleic Acids Res. 47, W234–W241. doi: 10.1093/nar/gkz240

PubMed Abstract | CrossRef Full Text | Google Scholar

Zhu, A. X., Finn, R. S., Edeline, J., Cattan, S., Ogasawara, S., Palmer, D., et al. (2018). Pembrolizumab in patients with advanced hepatocellular carcinoma previously treated with sorafenib (KEYNOTE-224): a non-randomised, open-label phase 2 trial. Lancet Oncol. 19, 940–952. doi: 10.1016/S1470-2045(18)30351-6

CrossRef Full Text | Google Scholar