Comprehensive Analysis of CDCAs Methylation and Immune Infiltrates in Hepatocellular Carcinoma

Background: As essential components of cycle growth, the cell division cycle-associated family genes (CDCAs) have crucial roles in tumor development and progression, especially in hepatocellular carcinoma (HCC). However, due to the tumor heterogeneity of HCC, little is known about the methylation variability of CDCAs in mediating phenotypic changes (e.g., immune infiltrates) in HCC. Presently, we aim to comprehensively explore the expression and prognosis of CDCAs methylation with regard to immune infiltrates of HCC.

Methods: We first identified the correlating differentially expressed genes (co-DEGs) among 19 different types of cancer cohorts (a total of 7,783 patients) and then constructed the weighted gene co-expressed and co-methylated networks. Applying the clustering analysis, significant modules of DEGs including CDCAs were selected and their functional bioinformatics analyses were performed. Besides, using DiseaseMeth and TIMER, the correlation between the methylation levels of CDCAs and tumor immune infiltrates was also analyzed. In final, to assess the influence of CDCAs methylation on clinical prognosis, Kaplan-Meier and Cox regression analysis were carried out.

Result: A total of 473 co-DEGs are successfully identified, while seven genes of CDCAs (CDCA1–3 and CDCA5–8) have significant over-expression in HCC. Co-expressed and co-methylated networks reveal the strong positive correlations in mRNA expression and methylation levels of CDCAs. Besides, the biological enrichment analysis of CDCAs demonstrates that they are significantly related to the immune function regulation of infiltrating immune cells in HCC. Also, the methylation analysis of CDCAs depicts the strong association with the tumor immunogenicity, i.e., low-methylation of CDCA1, CDCA2, and CDCA8 dramatically reduced the immune infiltrate levels of T cells and cytotoxic lymphocytes. Additionally, CDCA1–6 and CDCA8 with low-methylation levels significantly deteriorate the overall survival of patients in HCC.

Conclusions: The co-expressed and co-methylated gene networks of CDCAs show a powerful association with immune function regulation. And the methylation levels of CDCAs suggesting the prognostic value and infiltrating immune differences could be a novel and predictive biomarker for the response of immunotherapy.

Introduction

In accordance with the global cancer statistics released in 2018, hepatocellular carcinoma (HCC) has become the sixth most common cancer worldwide and the second leading cause of cancer-related death in men (1). Although several established treatments, including surgical resection, chemotherapy, and radiofrequency ablation, have been applied to the patient with HCC, the long-term survival of patients remains poor (2). Therefore, research into novel and effective treatment and prognostic signatures will need to be further performed.

Currently, the evolution and progression of tumors are considered to be related to the abnormal expression of various cancer-related genes, which involves different aspects of the life process in cells (e.g., cell cycle controlling, cell growth, cell apoptosis). Among them, the cell division cycle-associated family genes (CDCAs), consisting of eight members (CDCA1–8), are essential components of cell cycle growth, and they have been found malfunctions in various cancers, including HCC (3–6). For instance, CDCA1, also known as NUF2 (NUF2 component of NDC80 kinetochore complex), could inhibit tumor growth and induce apoptosis in HCC through silencing itself (7). And CDCA2 promoted the proliferation of tumor cells via activating the AKT/CCND1 pathway (8). Through enhancing cell proliferation with the prevention of G1 phase arrest, the upregulation of CDCA3 promoted cancer progression (9). With respect to CDCA6, also known as CBX2 (chromobox 2), it could become the regulator of the proliferation and apoptosis by the phosphorylation of YAP (yes1 associated transcriptional regulator) in HCC (10). Besides, CDCA8, as the key mediator of estrogen-stimulated cell proliferation in cancer cells, could inhibit cancer cell survival and growth by cell cycle G1 phase arrest (11).

Additionally, cancer immunotherapy has made tremendous success in cancer treatment, and a more detailed understanding of the immune infiltrate may be beneficial in developing the rationale for immunotherapy (12, 13). Increasing researches have indicated the relation between cancer immunotherapy and DNA methylation (14–16). For example, Gallagher et al. have reported that DNA methylation could contribute to cancer immunotherapy via modulating immune cell differentiation and function (15). Notably, several genes of CDCAs have also been involved in the immune infiltrates. Taking CDCA4 as an example, it is the E2F transcription factor, which plays an important role in immunity regulation (17, 18). Also, the expression of CDCA6 in immune cells (e.g. macrophages and T cells) could suggest the intensity of innate immune response (19). Moreover, CDCA7, associated with DNA methylation and cellular functions, has a regulation in the apoptosis of T cells (20, 21). However, despite the closed relation of the CDCAs in the occurrence and development of HCC, the association between the methylation of these genes and immune infiltrates in HCC has not been systematically investigated.

In this study, we obtained plenty of cancer datasets from TCGA and found significantly up-regulated CDCAs in HCC. Through constructing the co-expressed and co-methylated networks, we explored the biological enrichment of CDCAs and found that they are significantly related to immune function. Further, the immune infiltration and prognostic value between the high and low methylation groups of CDCAs were also compared. All the results demonstrate that the methylation levels of CDCAs could be served as predictive biomarkers for determining prognosis and immune infiltration in HCC. The flow diagram of this study was illustrated in Figure 1.

FIGURE 1

Figure 1 Workflow of this study. The high and low methylation groups are classified by the median methylation levels of CDCAs.

Materials and Methods

Differentially Expressed Genes (DEGs) Analysis

Publicly available expression data of HCC were collected from The Cancer Genome Atlas (TCGA, https://cancergenome.nih.gov/). After removing the datasets without enough normal examples, a total of 19 different types of cancer datasets (the detailed descriptions of cancers were listed in Figure 2) and 7,783 patients were finally obtained. Meanwhile, to identify the DEGs, we compared the tumor and normal samples of each dataset employing the “limma” R package (22), and the thresholds were set to FDR<0.01 and |log2FC|>2, respectively. Using the R package “RobustRankAggreg” (RRA) (23), the DEGs of these datasets were integrated to explore the correlating differentially expressed genes (co-DEGs). Notably, the genes with adjusted P < 0.05 were sorted by their log2FC into the up-regulated and down-regulated gene lists in the RRA analysis. Additionally, since members of the gene family have similar biological functions, we used the Hugo Gene Nomenclature Committee (HGNC) online database (https://www.genenames.org/) to build the gene families of DEGs and validate the results through literature checking.

FIGURE 2

Figure 2 Result of significant DEGs via the RRA analysis. Red shows up-regulation and green reveals down-regulation. The value is the logFC from the DEG analysis of each dataset. Additionally, seven genes of CDCAs (NUF2, CDCA2, CDCA3, CDCA5, CBX2, CDCA7, and CDCA8) are up-regulated in each dataset. And the full names of cancer types in x-arias are Bladder Urothelial Carcinoma (BLCA), Breast invasive carcinoma (BRCA), Cervical squamous cell carcinoma and endocervical adenocarcinoma (CESC), Cholangiocarcinoma (CHOL), Colon adenocarcinoma (COAD), Esophageal carcinoma (ESCA), Glioblastoma multiforme (GBM), Head and Neck squamous cell carcinoma (HNSC), Kidney renal clear cell carcinoma (KIRC), Kidney renal papillary cell carcinoma (KIRP), Lung adenocarcinoma (LUAD), Lung squamous cell carcinoma (LUSC), Pancreatic adenocarcinoma (PRAD), Rectum adenocarcinoma (READ), Sarcoma (SARC), Stomach adenocarcinoma (STAD), Thyroid carcinoma (THCA), Uterine Corpus Endometrial Carcinoma (UCEC), and Liver hepatocellular carcinoma (LIHC). logFC, logarithmic fold change; RRA, robust rank aggregation.

Weighted Gene Co-Expressed and Co-Methylated Networks Construction

The gene expression profile matrix file of patients with HCC, including 374 tumor samples, was publicly obtainable from TCGA. The co-expressed and co-methylated genes were analyzed by the “WGCNA” R package (24). Through setting the parameters of the network analysis (e.g., soft-thresholding = 12, scale-free R^2 = 0.90, and cut-height = 0.25), all the genes were divided into different gene modules according to the mRNA expression matrix. Meanwhile, using the cutoff of |log2FC|>2 and FDR<0.01 for the DEGs analysis in the liver hepatocellular carcinoma (LIHC) dataset, we collected the significant expression genes from the obtained gene modules.

Additionally, the methylation microarray data containing 380 tumor samples were collected from the TCGA-LIHC dataset. After removing the invalid and low-quality expression value, we corrected the deviation of methylation expression via the R package ChAMP v2.14.0. With the same parameters of the co-expressed network analysis, the methylation genes were divided into different gene modules. Employing the threshold values of |logFC|>0.2 and p < 0.05, we explored the differential methylation probes (DMP) to choose the significant methylation genes via the R package minfi v1.30.0. All the significant genes overlapping in the key expression module and methylation module (25) were used to perform further analysis.

Protein−Protein Interactions (PPI) Network Construction and Functional Enrichment Analysis

We established the functional protein interaction network for 71 co-expressed and co-methylated genes by using known and predicted interactions through the STRING database (26). Besides, we calculated the betweenness and degree (27) of all significant genes via the R package “igraph” v1.2.6. Furthermore, the Database for Annotation Visualization and Integrated Discovery (DAVID) (28) was performed to explore the biological enrichment of these genes. With setting the statistical significance as 0.05, the Gene Ontology (GO) terms and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway were visualized by R package “GOplot” (29).

Methylation Analysis

To investigate the methylation levels of CDCAs in normal and tumor samples of HCC, we used the online tools DiseaseMeth (30, 31). It integrates methylation data from several public sources and literature and provides information on the associations between diseases and genes. Herein, we set the technology experimental platform as 450k (Illumina Infinium HumanMethylation 450 BeadChip) and the significant p-value as 0.05. Moreover, to indicate whether the abnormal methylation has an effect on the expression level of CDCAs, we explored the correlation between their expression and methylation levels and visualized this result via the R package corrplot v0.84.

Immune Cell Infiltration Analysis

To explore the association between the mRNA expressions of CDCAs and immune cell infiltration, we used the website service TIMER (32, 33). The TIMER database validated by pathological estimations can systematically analyze the tumor immune infiltration across different cancer types.

Analyzing the Relationship Between Immune Infiltrates and CDCAs Methylation

Based on the median methylation levels of each CDCAs, the tumor samples were divided into the high and low methylation groups. To explore the tumor immune infiltration of HCC, we further compared the respective associations between the neoantigen load (33), seven classical immune signatures (34), immune-related gene expressions (35), and the methylation groups of CDCAs. We applied the R package MCP-counter v1.1.0 to assess abundances of tissue-infiltrating immune. And the R package GSVA v1.34.0 was applied to investigate single sample GSEA (ssGSEA) enrichment scores of these immune signatures.

Survival Analysis

To evaluate the associations between DNA methylation and survival of HCC patients, 370 samples with both the clinical and methylation information from the obtained HCC dataset were subjected to further perform cox regression analysis. Additionally, the Kaplan-Meier method was also applied to explore the differences in the survival time of the methylation groups, with the R packages survival v3.1-8 and survminer v0.4.6.

Statistical Analysis

The multivariate cox regression analysis was employed to investigate the independent prognostic value of other survival features. Log-rank and cox p-value analysis was performed to explore the statistical significance of observed differences in methylation groups. Additionally, the statistical significance was set as 0.05.

Result

mRNA Expression Levels of the CDCAs in HCC

According to the result of RRA analysis, we obtained 159 up-regulation and 314 down-regulation differential genes. Among them, seven genes of CDCAs (34) (CDCA1/NUF2: P = 2.73E-22, CDCA2: P = 1.52E-16, CDCA3: p = 1.50E-20, CDCA5: P = 1.77E-20, CDCA6/CBX2: P = 6.87E-19, CDCA7: P = 5.92E-16, and CDCA8: P = 1.67E-18) were all up-regulated in 19 cancer datasets (Figure 2). Especially, CDCAs are significantly up-regulated in HCC: CDCA1 (log2FC = 3.72, adjusted P = 6.70E-46), CDCA2 (log2FC = 2.76, adjusted P = 1.4E-29), CDCA3 (log2FC = 2.92, adjusted P = 1.96E-39), CDCA5 (log2FC = 3.15, adjusted P = 2.603E-45), CDCA6 (log2FC = 2.12, adjusted P = 3.79E-21), CDCA7 (logF2C = 2.29, adjusted P = 3.40E-10), and CDCA8 (log2FC = 2.86, adjusted P = 3.21E-39).

Subsequently, the mRNA expression levels of CDCAs were explored in the tumor and normal samples from the LIHC dataset. The results indicated that CDCA1–8 were over-expressed in cancer tissues than in normal tissues, and these genes all had significant differences (Figure 3).

FIGURE 3

Figure 3 The CDCAs have significant differences between HCC and normal tissue. Meanwhile, eight genes of CDCAs have over-expression in the tumor sample. T-test, ****p < 0.0001.

Analysis of Co-Expressed and Co-Methylated Networks

To find the co-expression genes of CDCAs, we applied the WGCNA on the LIHC dataset. From the analysis of the scale independence and mean connectivity, we chose the number 12 as the best soft-threshold power to construct the co-expression module (Figure 4A). Through the cut-height and minimal module size, we clustered the genes and cut the tree into nine gene modules (Figures 4B, C). Furthermore, the turquoise model contained a total of 2,961 genes and eight genes of CDCAs were all in this model. According to the DEGs analysis, a total of 178 genes were obtained by extracting the turquoise module. After taking 178 genes as the center nodes, we calculated the number of edges connected with CDCAs from the network of the turquoise module. As a result, a total of 126 genes with higher connections with CDCAs were obtained and then considered as the co-expression genes. Notably, the Gleason score (correlation coefficient = 0.41 and P = 2E-120) indicated that there was a reliably positive correlation between the module membership and gene significance in the turquoise module for HCC (Figure 4D).

FIGURE 4

Figure 4 Identify the key module of the LIHC dataset by WGCNA (A) Analysis of soft-threshold powers with the scale-free index and the mean connectivity. (B) Dendrogram of all filtered genes enrich according to a dissimilarity measure and the cluster module colors in the mRNA expression (minimum number of genes in the module is 30). (C) Adjust cut-height as 0.25 (the red line) to cluster the module eigengene. (D) Correlation of the co-expressed genes in the turquoise model. Each dot is one gene of the module.

In addition, the co-methylation network analysis reveals that the methylation genes of HCC were divided into 11 modules (Figure 5A). Particularly, five genes of CDCAs were found in the module colored as “brown” (a total of 1,859 genes). Through the DMP and connection analysis, we got 117 co-methylated genes from the brown module. Furthermore, 71 genes overlapping in the turquoise and brown modules, are regarded as the co-expressed and co-methylated genes. The heatmap showed that there were positive correlations of these genes in the mRNA expression and DNA methylation values (Figures 5B, C).

FIGURE 5

Figure 5 (A) Cluster dendrogram of DNA methylation values into 11 modules. (B) Heatmap shows the correlation of the top 35 genes from co-expressed and co-methylated genes according to the standard deviation of the mRNA expression. (C) Heatmap shows the correlation of the top 35 genes from co-expressed and co-methylated genes according to the standard deviation of the DNA methylation values. (D) The protein-protein interaction network for co-expressed and co-methylated genes. (E) The degree and betweenness of 29 significant genes. The bar graph shows the degree, and the line graph indicates the betweenness.

Function Enrichment Analysis

After removing the disconnected nodes, the protein interaction network of the co-expressed and co-methylated genes was constructed, consisting of 29 genes and 243 edges (Figure 5D). The mean value of the targets per gene is 8.3, which indicated the diversity of these significant genes in biological targeting. Additionally, it is clear that genes with higher degrees would also have higher betweenness (Figure 5E), displaying the strong positive correlation between the degree and betweenness (Supplementary Figure 1). Notably, NUF2, CDCA5, and CDCA8 had the highest degree and betweenness, which are considered as the center nodes of this network. The GO enrichment analysis reveals that the functions of these target genes are associated with the biological process, cellular component, and molecular functions. Notably, we obtained several significant enrichment terms from the analysis of biological process, such as cell cycle checkpoint and mitotic nuclear division (Figure 6A). In addition, the cellular component analysis showed the enrichment of these genes in the chromosome region and condensed chromosome (Figure 6B). Also, it is clear that protein serine/threonine kinase activity was the most significant enrichment term in the analysis of molecular functions (Figure 6C). As to the KEGG pathway analysis, cell cycle, p53 signaling pathway, HEpatitis B, and viral carcinogenesis had strong associations with these genes (Figure 6D). The results of the PPI network and enrichment analysis show the potential association between CDCAs and the immune function. For example, CDCAs have the higher connectivity with the immune-related genes [e.g., CDK1 (36, 37) and PLK1 (38, 39)] in the network, and also significantly enrich in the immune-related terms [e.g., p53 signaling pathway (40–42) and viral carcinogenesis (43–45)].

FIGURE 6

Figure 6 Go and KEGG analysis of 71 co-expressed and co-methylated genes. (A) Terms for biological process. (B) Terms for cellular component. (C) Terms for molecular functions. (D) Pathways for KEGG analysis.

Methylation Difference of CDCAs

In accordance with the result of DiseaseMeth (Figure 7), the methylation levels of CDCA1 (P = 6.59E-13), CDCA3 (P = 2.78E-08), CDCA4 (P = 3.22e-02), CDCA5 (P = 1.10E-06), CDCA6 (P = 1.13E-03), and CDCA8 (P = 8.77e-13) are all significantly higher in normal samples than the levels in disease samples, whereas CDCA7 (P = 1.66E-15) was significantly higher in disease samples. Besides, CDCA2 (P = 5.04E-02) had no significant difference in the sample groups. Additionally, further analysis demonstrates that CDCAs have a consistently negative correlation between the expression and methylation levels (Supplementary Figure 2).

FIGURE 7

Figure 7 Methylation analysis of CDCAs via DiseaseMeth 2.0. Methylation values for (A) CDCA1, (B) CDCA2, (C) CDCA3, (D) CDCA4, (E) CDCA5, (F) CDCA6, (G) CDCA7, and (H) CDCA8.

Transcriptional Levels of CDCAs Have Strongly Positive Correlations With the Immune Infiltration

The online tool TIMER was applied to investigate potential relations between the mRNA expression levels of CDCAs and the states of the tumor purity and immune cell infiltration. The strongly positive associations (Supplementary Figure 3) were observed between CDCA1–8 and six types of immune infiltrates (e.g., B cells and dendritic cells). Meanwhile, CDCA1–5 and CDCA8 demonstrated the weak correlations in the tumor purity, whereas the associations of CDCA6 and CDCA7 were weak and negative (Figure 8).

FIGURE 8

Figure 8 Correlation heatmap between the expression levels of CDCAs and the immune infiltration levels. The number shows the correlation value. ns, no statistical significance; *p < 0.05; **p < 0.01; ***p < 0.001.

Association of CDCAs Methylation With Enhancing Tumor Immunogenicity

To explore cancer immune infiltrates in HCC, we investigated the tumor immunogenicity between methylation groups of CDCAs. The neoantigen load had significant differences in CDCA1, CDCA2, and CDCA8, which indicated that these genes with low-methylation had powerful associations with weakening tumor immunogenicity (Figure 9A). Additionally, several tumor-infiltrating parameters, especially the T cells and cytotoxic lymphocytes, were generally more abundant in the high-methylation samples than those in low-methylation samples across CDCAs (Figure 9B). And the results of the immune signature analysis suggest that the type I and type II IFN responses are almost higher in all the CDCAs with high-methylation (Figure 9C). Moreover, it also showed that CDCA1, CDCA2, and CDCA8 with high-methylation were more abundant in seven immune signatures. To further estimate the immune profile, we profoundly analyzed the different expression patterns of immune-related genes in methylation states of CDCAs (Figure 9D). In consonance with the result of immune infiltrates and gene signatures, different stimulatory immunomodulators had the general up-regulation in CDCA1, CDCA2, and CDCA8 with high-methylation, such as chemokines (CCL5, CX3CL1, CXCL10, CXCL9) and human leucocyte antigen (HLA-A, HLA-DPA1, HLA-DQA1). All these results reveal that the methylation states of CDCA1, CDCA2, and CDCA8 have strong associations in the tumor immunogenicity.

FIGURE 9

Figure 9 The high-methylation CDCAs have associations with enhancing tumor immunogenicity. (A) The difference of neoantigen load between methylation groups of CDCAs in the LIHC cohorts (only show the result with significant differences). (B) Heatmap visualizes the log2FC between methylation groups in the MCP-counter scores across CDCAs. (C) Heatmap visualizes the log2FC between methylation groups in the ssGSEA scores across CDCAs. (D) Heatmap visualizes the log2FC between high and low methylation groups in the immune-related mRNA expressions across CDCAs. The y-axis indicates different immune-related genes and their function and category types. Red suggests up-regulation, while blue suggests down-regulation.

CDCAs With Low-Methylation Deteriorate the Overall Survival in HCC

The cox regression analyses were conducted to explore the contribution of CDCAs with methylation groups as the independent prognostic signatures of patient survival. Based on the clinical factors, including age, gender, and pathological stage, the result confirms that the predictive ability of CDCAs methylation, including CDCA1–6 and CDCA8, were independent of other clinical characteristics for overall survival in HCC (Table 1).

TABLE 1

Table 1 Cox regression analysis of CDCAs methylation groups and clinical characteristics of patients in the LIHC dataset.

The Kaplan-Meier plotter was applied to visualize the survival of CDCAs methylation groups in the LIHC cohorts. Furthermore, the log-rank and cox tests were applied to evaluate the result. As shown in Figure 10, the patients with the high-methylation levels of CDCAs, including CDCA1–6 and CDCA8, extensively had a longer OS than the low-methylation counterparts. Besides, the p-value of the log-rank and cox analysis indicated that the result had significant differences and effectiveness.

FIGURE 10

Figure 10 Kaplan-Meier curves for CDCAs methylation groups. Survival analyses for (A) CDCA1, (B) CDCA2, (C) CDCA3, (D) CDCA4, (E) CDCA5, (F) CDCA6, (G) CDCA7, and (H) CDCA8. Patient samples are dichotomized based on the median methylation levels of each CDCAs. P-value refers to the cox and log-rank tests. HR shows the hazard ratio and 95% confidence interval.

Discussion

The success of immune therapy in cancer treatment has been widely reported (12, 13). As an important strategy for immunotherapy, DNA methylation could be served as a promising target in different cancers, especially in HCC (14, 46, 47). In addition, CDCAs play crucial roles in the development and progression of cancers, such as HCC (7, 10). However, little is known about the relationship between the immune infiltrates and DNA methylation of CDCAs.

Presently, we first explored the mRNA expression data in 19 different types of cancer datasets. As shown in Figures 2 and 3, the results indicate that most genes of CDCAs, including CDCA1–3 and CDCA5–8, are significantly up-regulated in diverse malignant cancers, especially in HCC. Indeed, the abnormal expression of CDCAs has been reported in the tumor tissue of HCC (48–50). Additionally, the co-expressed and co-methylated gene networks analysis shows that the biological enrichments of CDCAs are significantly associated with immune functions, such as p53 signaling pathway and viral carcinogenesis (Figures 4–6). Actually, it is reported that p53 could stimulate the innate immune system to maintain tissue homeostasis and suppress tumorigenesis (40), and CDCA2 also plays a critical role in the activation of p53 (41). Besides, CDCA4 has involvement in regulating the transcriptional activities of p53, and its over-expression could lead to p53-independent growth inhibition (42). All the results demonstrate that CDCAs have a potential association with immune functions.

Moreover, to explore the potential relation of CDCAs and infiltrating immune cells, we analyzed their expression modes and found the significant and positive correlations between the expression levels of CDCAs and proportions of infiltrating immune cells (Figure 8). Consistently, it is reported that the knockdown of CDCA5 shows strong inhibition in tumor proliferation, which could be an attractive target for immunotherapy (51). Besides, CBX2 could promote the production of type I interferon in macrophages, which shows the regulation in antiviral immune response (19). Also, the expression of CDCA6 in immune cells (e.g. macrophages and T cells) suggests the intensity of innate immune response (19). In addition, tumor heterogeneity could be a significant obstacle for the tumor diagnosis and immunotherapy treatment in HCC, especially the methylations variability (46, 52–54).

Additionally, our study indicated that CDCAs (i.e., CDCA1, CDCA3–6, and CDCA8) have significant hypomethylation levels in tumor tissues of HCC (Figure 7). Notably, the consistently negative correlation between the expression and methylation levels of CDCAs (Supplementary Figure 2) shows that the abnormal methylation can cause the expression complexity of CDCAs, which play important roles in regulating the biological functions of CDCAs. In fact, the methylation of CDCAs has been partly reported in HCC, but the association with infiltrating immune would still be unclear (20, 55, 56). For instance, Cai et al (55). have demonstrated that CDCA5, as a differentially methylated gene, can be served as a new novel therapeutic target for HCC. In another study, it is also revealed that the abnormal methylation and over-expression of CDCA8 plays an essential role in the regulation of cell cycle (56, 57), while the dysregulated cell cycle is a significant hallmark of immune checkpoint pathways for immunotherapy (58). With respect to CDCA7, it has been reported that CDCA7 can be involved in the regulation of T cells, through dynamic DNA methylation (20, 21). Thus, we explored the link between CDCAs methylation and immune states in HCC and found that the methylation levels of CDCAs have strong associations in tumor immunogenicity (Figure 9). For example, the neoantigen load, as the mutated protein from the somatic mutation of tumors, participates in the function of immune responses, which would be considered as a novel marker for immunity therapy (59–61). Herein, it is indicated that CDCAs, especially CDCA1, CDCA2, and CDCA8, have higher levels of neoantigen load in the high methylation group (Figure 9A). In addition, the immune cell infiltration and significant markers of immune cells (Figures 9B–D) suggest the improvement of immune activity under the high methylation status of CDCAs, which is consistent with the obvious result of neoantigen load. Besides, the survival curve analysis also shows that CDCA1–6 and CDCA8 with low-methylation levels significantly deteriorate the overall survival of patients in HCC (Figure 10), further suggesting that methylation levels of CDCAs can be regarded as the predictive biomarkers for determining prognosis and immune infiltration in HCC. All the results provide valuable evidence to the combination of the DNA methylation of CDCAs and immune infiltration in HCC.

However, our study also has some limitations. Due to the data type requirements, including mRNA expression, methylation expression, and neoantigen load calculation, we only obtained the data from TCGA, which may cause the data bias of this investigation. Therefore, more tumor samples and further experimental validation are necessary to perform for evaluating the biological roles of CDCAs in HCC.

Conclusion

Collectively, we identified 71 co-expressed and co-methylated genes of CDCAs that were enriched in biological terms related to immune function. And we integratively analyzed the association between the immune infiltrates and methylation of CDCAs. Our results demonstrate that the methylation groups of CDCAs are valuable in predicting prognosis. Thus, the methylation of CDCAs could be novel predictive biomarkers for immune infiltrates and be used in the clinical diagnosis and treatment of HCC.

Data Availability Statement

The datasets presented in this study can be found in online repositories. The names of the repository/repositories and accession number(s) can be found in the article/Supplementary Material.

Author Contributions

YW conducted the bioinformatics analysis, prepared the figures, and wrote the paper. YW, assisted by HG, TO, LZ, JH, and SH, contributed to the experiments of the manuscript. YY and HK conceived the idea, designed the study, and coordinated the project. All authors contributed to the article and approved the submitted version.

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.

Acknowledgments

Thanks for the High-level talent introduction project of Anhui University of Chinese Medicine (Grant no. DT19200011) and the Natural Science Key Project of Anhui University of Chinese Medicine (Grant no. 2020zrzd15), and the Natural Science Key Project of Anhui universities (Grant no. KJ2020A0419).

Supplementary Material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fonc.2020.566183/full#supplementary-material

References

1. Bray F, Ferlay J, Soerjomataram I, Siegel RL, Torre LA, Jemal A. Global cancer statistics 2018: GLOBOCAN estimates of incidence and mortality worldwide for 36 cancers in 185 countries. CA: A Cancer J Clin (2018) 68(6):394–424. doi: 10.3322/caac.21492

CrossRef Full Text | Google Scholar

2. Heimbach JK, Kulik LM, Finn RS, Sirlin CB, Abecassis MM, Roberts LR, et al. AASLD guidelines for the treatment of hepatocellular carcinoma: Heimbach et al. Hepatology (2018) 67(1):358–80. doi: 10.1002/hep.29086

PubMed Abstract | CrossRef Full Text | Google Scholar

3. Xie HJ, Noh JH, Kim JK, Jung KH, Eun JW, Bae HJ, et al. HDAC1 Inactivation Induces Mitotic Defect and Caspase-Independent Autophagic Cell Death in Liver Cancer. PloS One (2012) 7(4):e34265. doi: 10.1371/journal.pone.0034265

PubMed Abstract | CrossRef Full Text | Google Scholar

4. Yoshino Y, Ishioka C. Inhibition of glycogen synthase kinase-3 beta induces apoptosis and mitotic catastrophe by disrupting centrosome regulation in cancer cells. srep (2015) 5:13249. doi: 10.1038/srep13249

CrossRef Full Text | Google Scholar

5. Khongkow P, Gomes AR, Gong C, Man EPS, Tsang JW-H, Zhao F, et al. Paclitaxel targets FOXM1 to regulate KIF20A in mitotic catastrophe and breast cancer paclitaxel resistance. Oncogene (2016) 35(8):990–1002. doi: 10.1038/onc.2015.152

PubMed Abstract | CrossRef Full Text | Google Scholar

6. Osman AA, Monroe MM, Alves MVO, Patel AA, Katsonis P, Fitzgerald AL, et al. Wee-1 Kinase Inhibition Overcomes Cisplatin Resistance Associated with High-Risk TP53 Mutations in Head and Neck Cancer through Mitotic Arrest Followed by Senescence. Mol Cancer Ther (2015) 14(2):608–19. doi: 10.1158/1535-7163.MCT-14-0735-T

PubMed Abstract | CrossRef Full Text | Google Scholar

7. Liu Q, Dai S-J, Li H, Dong L, Peng Y-P. Silencing of NUF2 Inhibits Tumor Growth and Induces Apoptosis in Human Hepatocellular Carcinomas. Asian Pacific J Cancer Prev (2014) 15(20):8623–9. doi: 10.7314/APJCP.2014.15.20.8623

CrossRef Full Text | Google Scholar

8. Feng Y, Qian W, Zhang Y, Peng W, Li J, Gu Q, et al. CDCA2 promotes the proliferation of colorectal cancer cells by activating the AKT/CCND1 pathway in vitro and in vivo. BMC Cancer (2019) 19(1):576. doi: 10.1186/s12885-019-5793-z

PubMed Abstract | CrossRef Full Text | Google Scholar

9. Uchida F, Uzawa K, Kasamatsu A, Takatori H, Sakamoto Y, Ogawara K, et al. Overexpression of cell cycle regulator CDCA3 promotes oral cancer progression by enhancing cell proliferation with prevention of G1 phase arrest. BMC Cancer (2012) 12(1):321. doi: 10.1186/1471-2407-12-321

PubMed Abstract | CrossRef Full Text | Google Scholar

10. Mao J, Tian Y, Wang C, Jiang K, Li R, Yao Y, et al. CBX2 Regulates Proliferation and Apoptosis via the Phosphorylation of YAP in Hepatocellular Carcinoma. J Cancer (2019) 10(12):2706–19. doi: 10.7150/jca.31845

PubMed Abstract | CrossRef Full Text | Google Scholar

11. Bu Y, Shi L, Yu D, Liang Z, Li W. CDCA8 is a key mediator of estrogen-stimulated cell proliferation in breast cancer cells. Gene (2019) 703:1–6. doi: 10.1016/j.gene.2019.04.006

PubMed Abstract | CrossRef Full Text | Google Scholar

12. Pardoll DM. The blockade of immune checkpoints in cancer immunotherapy. Nat Rev Cancer (2012) 12(4):252–64. doi: 10.1038/nrc3239

PubMed Abstract | CrossRef Full Text | Google Scholar

13. Mandal R, Şenbabaoğlu Y, Desrichard A, Havel JJ, Dalin MG, Riaz N, et al. The head and neck cancer immune landscape and its immunotherapeutic implications. JCI Insight (2016) 1(17):e89829. doi: 10.1172/jci.insight.89829

PubMed Abstract | CrossRef Full Text | Google Scholar

14. Jones PA, Ohtani H, Chakravarthy A, De Carvalho DD. Epigenetic therapy in immune-oncology. Nat Rev Cancer (2019) 19(3):151–61. doi: 10.1038/s41568-019-0109-9

PubMed Abstract | CrossRef Full Text | Google Scholar

15. Gallagher SJ, Shklovskaya E, Hersey P. Epigenetic modulation in cancer immunotherapy. Curr Opin Pharmacol (2017) 35:48–56. doi: 10.1016/j.coph.2017.05.006

PubMed Abstract | CrossRef Full Text | Google Scholar

16. Accomando WP, Wiencke JK, Houseman E, Nelson HH, Kelsey KT. Quantitative reconstruction of leukocyte subsets using DNA methylation. Genome Biol (2014) 15(3):R50. doi: 10.1186/gb-2014-15-3-r50

PubMed Abstract | CrossRef Full Text | Google Scholar

17. Hayashi R, Goto Y, Ikeda R, Yokoyama KK, Yoshida K. CDCA4 Is an E2F Transcription Factor Family-induced Nuclear Factor That Regulates E2F-dependent Transcriptional Activation and Cell Proliferation. J Biol Chem (2006) 281(47):35633–48. doi: 10.1074/jbc.M603800200

PubMed Abstract | CrossRef Full Text | Google Scholar

18. Ren Z, Kang W, Wang L, Sun B, Ma J, Zheng C, et al. E2F1 renders prostate cancer cell resistant to ICAM-1 mediated antitumor immunity by NF-κB modulation. Mol Cancer (2014) 13(1):84. doi: 10.1186/1476-4598-13-84

PubMed Abstract | CrossRef Full Text | Google Scholar

19. Sun D, Cao X, Wang C. Polycomb chromobox Cbx2 enhances antiviral innate immunity by promoting Jmjd3-mediated demethylation of H3K27 at the Ifnb promoter. Protein Cell (2019) 10(4):285–94. doi: 10.1007/s13238-018-0581-0

PubMed Abstract | CrossRef Full Text | Google Scholar

20. Qin X-Y, Feng J, Chen G, Dou X-W, Dai X-Q, Dong H-L, et al. ZBTB24 regulates the apoptosis of human T cells via CDCA7/TRAIL-receptor axis. Biochem Biophys Res Commun (2019) 514(1):259–65. doi: 10.1016/j.bbrc.2019.04.147

PubMed Abstract | CrossRef Full Text | Google Scholar

21. Thompson JJ, Kaur R, Sosa CP, Lee J-H, Kashiwagi K, Zhou D, et al. ZBTB24 is a transcriptional regulator that coordinates with DNMT3B to control DNA methylation. Nucleic Acids Res (2018) 46(19):10034–51. doi: 10.1093/nar/gky682

PubMed Abstract | CrossRef Full Text | Google Scholar

22. Ritchie ME, Phipson B, Wu D, Hu Y, Law CW, Shi W, et al. limma powers differential expression analyses for RNA-sequencing and microarray studies. Nucleic Acids Res (2015) 43(7):e47–7. doi: 10.1093/nar/gkv007

PubMed Abstract | CrossRef Full Text | Google Scholar

23. Kolde R, Laur S, Adler P, Vilo J. Robust rank aggregation for gene list integration and meta-analysis. Bioinformatics (2012) 28(4):573–80. doi: 10.1093/bioinformatics/btr709

PubMed Abstract | CrossRef Full Text | Google Scholar

24. Langfelder P, Horvath S. WGCNA: an R package for weighted correlation network analysis. BMC Bioinf (2008) 9:559. doi: 10.1186/1471-2105-9-559

CrossRef Full Text | Google Scholar

25. Wang J, Li Y, Yang Y, Du J, Zhao M, Lin F, et al. Systems Pharmacology Dissection of Multiscale Mechanisms of Action for Herbal Medicines in Treating Rheumatoid Arthritis. Mol Pharmaceutics (2017) 14(9):3201–17. doi: 10.1021/acs.molpharmaceut.7b00505

CrossRef Full Text | Google Scholar

26. Szklarczyk D, Gable AL, Lyon D, Junge A, Wyder S, Huerta-Cepas J, et al. STRING v11: protein-protein association networks with increased coverage, supporting functional discovery in genome-wide experimental datasets. Nucleic Acids Res (2019) 47(D1):D607–13. doi: 10.1093/nar/gky1131

PubMed Abstract | CrossRef Full Text | Google Scholar

27. Yang Y, Li Y, Wang J, Sun K, Tao W, Wang Z, et al. Systematic Investigation of Ginkgo Biloba Leaves for Treating Cardio-cerebrovascular Diseases in an Animal Model. ACS Chem Biol (2017) 12(5):1363–72. doi: 10.1021/acschembio.6b00762

PubMed Abstract | CrossRef Full Text | Google Scholar

28. Huang DW, Sherman BT, Lempicki RA. Systematic and integrative analysis of large gene lists using DAVID bioinformatics resources. Nat Protoc (2009) 4(1):44–57. doi: 10.1038/nprot.2008.211

PubMed Abstract | CrossRef Full Text | Google Scholar

29. Walter W, Sánchez-Cabo F, Ricote M. GOplot: an R package for visually combining expression data with functional analysis. Bioinformatics (2015) 31(17):2912–4. doi: 10.1093/bioinformatics/btv300

PubMed Abstract | CrossRef Full Text | Google Scholar

30. Lv J, Liu H, Su J, Wu X, Liu H, Li B, et al. DiseaseMeth: a human disease methylation database. Nucleic Acids Res (2012) 40(Database issue):D1030–1035. doi: 10.1093/nar/gkr1169

PubMed Abstract | CrossRef Full Text | Google Scholar

31. Xiong Y, Wei Y, Gu Y, Zhang S, Lyu J, Zhang B, et al. DiseaseMeth version 2.0: a major expansion and update of the human disease methylation database. Nucleic Acids Res (2017) 45(D1):D888–95. doi: 10.1093/nar/gkw1123

PubMed Abstract | CrossRef Full Text | Google Scholar

32. Li B, Severson E, Pignon J-C, Zhao H, Li T, Novak J, et al. Comprehensive analyses of tumor immunity: implications for cancer immunotherapy. Genome Biol (2016) 17(1):174. doi: 10.1186/s13059-016-1028-7

PubMed Abstract | CrossRef Full Text | Google Scholar

33. Li T, Fan J, Wang B, Traugh N, Chen Q, Liu JS, et al. TIMER: A Web Server for Comprehensive Analysis of Tumor-Infiltrating Immune Cells. Cancer Res (2017) 77(21):e108–10. doi: 10.1158/0008-5472.CAN-17-0307

PubMed Abstract | CrossRef Full Text | Google Scholar

34. Rooney MS, Shukla SA, Wu CJ, Getz G, Hacohen N. Molecular and Genetic Properties of Tumors Associated with Local Immune Cytolytic Activity. Cell (2015) 160(1-2):48–61. doi: 10.1016/j.cell.2014.12.033

PubMed Abstract | CrossRef Full Text | Google Scholar

35. Thorsson V, Gibbs DL, Brown SD, Wolf D, Bortone DS, Ou Yang T-H, et al. The Immune Landscape of Cancer. Immunity (2018) 48(4):812–30.e14. doi: 10.1016/j.immuni.2018.03.023

PubMed Abstract | CrossRef Full Text | Google Scholar

36. Laphanuwat P, Jirawatnotai S. Immunomodulatory Roles of Cell Cycle Regulators. Front Cell Dev Biol (2019) 7:23. doi: 10.3389/fcell.2019.00023

PubMed Abstract | CrossRef Full Text | Google Scholar

37. Enserink JM, Kolodner RD. An overview of Cdk1-controlled targets and processes. Cell Div (2010) 5(1):11. doi: 10.1186/1747-1028-5-11

PubMed Abstract | CrossRef Full Text | Google Scholar

38. Li M, Liu Z, Wang X. Exploration of the Combination of PLK1 Inhibition with Immunotherapy in Cancer Treatment. J Oncol (2018) 2018:1–13. doi: 10.1155/2018/3979527

CrossRef Full Text | Google Scholar

39. Raab M, Strebhardt K, Rudd CE. Immune adaptor SKAP1 acts a scaffold for Polo-like kinase 1 (PLK1) for the optimal cell cycling of T-cells. Sci Rep (2019) 9(1):10462. doi: 10.1038/s41598-019-45627-9

PubMed Abstract | CrossRef Full Text | Google Scholar

40. Raj N, Attardi LD. Tumor Suppression: p53 Alters Immune Surveillance to Restrain Liver Cancer. Curr Biol (2013) 23(12):R527–30. doi: 10.1016/j.cub.2013.04.076

PubMed Abstract | CrossRef Full Text | Google Scholar

41. Shang D, Han T, Xu X, Liu Y. Decitabine induces G2/M cell cycle arrest by suppressing p38/NF-κB signaling in human renal clear cell carcinoma. Int J Clin Exp Pathol (2015) 8(9):11140–8.

PubMed Abstract | Google Scholar

42. Tategu M, Nakagawa H, Hayashi R, Yoshida K. Transcriptional co-factor CDCA4 participates in the regulation of JUN oncogene expression. Biochimie (2008) 90(10):1515–22. doi: 10.1016/j.biochi.2008.05.014

PubMed Abstract | CrossRef Full Text | Google Scholar

43. Chen Y, Williams V, Filippova M, Filippov V, Duerksen-Hughes P. Viral Carcinogenesis: Factors Inducing DNA Damage and Virus Integration. Cancers (2014) 6(4):2155–86. doi: 10.3390/cancers6042155

PubMed Abstract | CrossRef Full Text | Google Scholar

44. Tarocchi M. Molecular mechanism of hepatitis B virus-induced hepatocarcinogenesis. WJG (2014) 20(33):11630. doi: 10.3748/wjg.v20.i33.11630

PubMed Abstract | CrossRef Full Text | Google Scholar

45. Dandri M, Petersen J. Mechanism of Hepatitis B Virus Persistence in Hepatocytes and Its Carcinogenic Potential. Clin Infect Dis (2016) 62(suppl 4):S281–8. doi: 10.1093/cid/ciw023

PubMed Abstract | CrossRef Full Text | Google Scholar

46. Pardee AD, Butterfield LH. Immunotherapy of hepatocellular carcinoma. OncoImmunology (2012) 1(1):48–55. doi: 10.4161/onci.1.1.18344

PubMed Abstract | CrossRef Full Text | Google Scholar

47. Zhang Y, Petropoulos S, Liu J, Cheishvili D, Zhou R, Dymov S, et al. The signature of liver cancer in immune cells DNA methylation. Clin Epigenet (2018) 10(1):8. doi: 10.1186/s13148-017-0436-1

CrossRef Full Text | Google Scholar

48. Jin B, Wang W, Du G, Huang G-Z, Han L-T, Tang Z-Y, et al. Identifying hub genes and dysregulated pathways in hepatocellular carcinoma. Eur Rev Med Pharmacol Sci (2015) 19(4):592–601.

PubMed Abstract | Google Scholar

49. Wei YM, Li YY, Zhang YC, Nie YQ. Aberrant gene expression profiles in hepatocellular carcinoma detected by microdissection. Genet Mol Res (2013) 12(4):5527–36. doi: 10.4238/2013.November.18.3

PubMed Abstract | CrossRef Full Text | Google Scholar

50. Wang J, Xia C, Pu M, Dai B, Yang X, Shang R, et al. Silencing of CDCA5 inhibits cancer progression and serves as a prognostic biomarker for hepatocellular carcinoma. Oncol Rep (2018) 40(4):1875–84. doi: 10.3892/or.2018.6579

PubMed Abstract | CrossRef Full Text | Google Scholar

51. Xu J, Zhu C, Yu Y, Wu W, Cao J, Li Z, et al. Systematic cancer-testis gene expression analysis identified CDCA5 as a potential therapeutic target in esophageal squamous cell carcinoma. EBioMedicine (2019) 46:54–65. doi: 10.1016/j.ebiom.2019.07.030

PubMed Abstract | CrossRef Full Text | Google Scholar

52. Prieto J, Melero I, Sangro B. Immunological landscape and immunotherapy of hepatocellular carcinoma. Nat Rev Gastroenterol Hepatol (2015) 12(12):681–700. doi: 10.1038/nrgastro.2015.173

PubMed Abstract | CrossRef Full Text | Google Scholar

53. Lin D-C, Mayakonda A, Dinh HQ, Huang P, Lin L, Liu X, et al. Genomic and Epigenomic Heterogeneity of Hepatocellular Carcinoma. Cancer Res (2017) 77(9):2255–65. doi: 10.1158/0008-5472.CAN-16-2822

PubMed Abstract | CrossRef Full Text | Google Scholar

54. Sceusi EL, Loose DS, Wray CJ. Clinical implications of DNA methylation in hepatocellular carcinoma. HPB (2011) 13(6):369–76. doi: 10.1111/j.1477-2574.2011.00303.x

PubMed Abstract | CrossRef Full Text | Google Scholar

55. Cai C, Wang W, Tu Z. Aberrantly DNA Methylated-Differentially Expressed Genes and Pathways in Hepatocellular Carcinoma. J Cancer (2019) 10(2):355. doi: 10.7150/jca.27832

PubMed Abstract | CrossRef Full Text | Google Scholar

56. Dai C, Miao C-X, Xu X-M, Liu L-J, Gu Y-F, Zhou D, et al. Transcriptional Activation of Human CDCA8 Gene Regulated by Transcription Factor NF-Y in Embryonic Stem Cells and Cancer Cells. J Biol Chem (2015) 290(37):22423–34. doi: 10.1074/jbc.M115.642710

PubMed Abstract | CrossRef Full Text | Google Scholar

57. Lu T, Klein KO, Colmegna I, Lora M, Greenwood CMT, Hudson M. Whole-genome bisulfite sequencing in systemic sclerosis provides novel targets to understand disease pathogenesis. BMC Med Genomics (2019) 12(1):144. doi: 10.1186/s12920-019-0602-8

PubMed Abstract | CrossRef Full Text | Google Scholar

58. Hanahan D, Weinberg RA. Hallmarks of Cancer: The Next Generation. Cell (2011) 144(5):646–74. doi: 10.1016/j.cell.2011.02.013

PubMed Abstract | CrossRef Full Text | Google Scholar

59. Giannakis M, Mu XJ, Shukla SA, Qian ZR, Cohen O, Nishihara R, et al. Genomic Correlates of Immune-Cell Infiltrates in Colorectal Carcinoma. Cell Rep (2016) 15(4):857–65. doi: 10.1016/j.celrep.2016.03.075

PubMed Abstract | CrossRef Full Text | Google Scholar

60. McGranahan N, Furness AJS, Rosenthal R, Ramskov S, Lyngaa R, Saini SK, et al. Clonal neoantigens elicit T cell immunoreactivity and sensitivity to immune checkpoint blockade. Science (2016) 351(6280):1463–9. doi: 10.1126/science.aaf1490

PubMed Abstract | CrossRef Full Text | Google Scholar

61. Matsushita H, Sato Y, Karasaki T, Nakagawa T, Kume H, Ogawa S, et al. Neoantigen Load, Antigen Presentation Machinery, and Immune Signatures Determine Prognosis in Clear Cell Renal Cell Carcinoma. Cancer Immunol Res (2016) 4(5):463–71. doi: 10.1158/2326-6066.CIR-15-0225

PubMed Abstract | CrossRef Full Text | Google Scholar