Assessment of alterations in histone modification function and guidance for death risk prediction in cervical cancer patients

Background: Cervical cancer is the second most lethal malignancy among women, and histone modification plays a fundamental role in most biological processes, but the prognostic value of histone modification in cervical cancer has not been evaluated.

Methods: A total of 594 cervical cancer patients from TCGA-CESC, GSE44001, and GSE52903 cohorts were enrolled in the current study, along with the corresponding clinicopathological features. Patients with a follow-up time less than one month were removed. A total of 122 histone modification-associated signaling pathways were obtained from the MSigDB. The activation scores of these pathways were evaluated using the “GSVA” package, differentially expressed genes were identified by the “limma” package, and pathway enrichment was conducted using the “clusterProfiler 4.0” package. The subsequent least absolute shrinkage and selection operator (LASSO) regression analysis was performed using the “glmnet” package, and a prognostic nomogram was established using the “regplot” package. For the prediction of potential therapeutic drugs, we used the data from GDSC2016 and visualized them via “MOVICS”.

Results: Nine of 23 histone modification-associated prognostic genes were identified to construct the prognostic signature by LASSO analysis, named the histone modification-associated gene (HMAG) signature. Cervical patients with HMAG-H in TCGA-CESC cohort showed a 2.68-fold change of death risk, with the 95% CI from 1.533 to 4.671 (p < 0.001), as well as the increased death risk of HMAG-H in the GSE44001 cohort (HR: 2.83, 95% CI: 1.370–5.849, p = 0.005) and GSE44001 cohort (HR: 4.59, 95% CI: 1.658–12.697, p = 0.003). We observed the preferable AUC values of the HMAG signature in TCGA-CESC cohort (1-year: 0.719, 3-year: 0.741, and 5-year: 0.731) and GSE44001 cohort (1-year: 0.850, 3-year: 0.781, and 5-year: 0.755). The C-index of the nomogram showed a prognostic value as high as 0.890, while the C-index for age was only 0.562, and that for grade was only 0.542. Patients with high HMAG scores were more suitable for the treatment of CHIR-99021, embelin, FTI-277, JNK-9L, JQ12, midostaurin, PF-562271, pyrimethamine, and thapsigargin, and patients with low HMAG scores were more suitable for the treatment of BMS-536924, CP466722, crizotinib, PHA-665752, rapamycin, and TAE684.

Conclusion: We comprehensively evaluated the histone modification status in cervical cancer patients and revealed histone modification-associated prognostic genes to construct the HMAG signature, aiming to provide a new insight into prognosis prediction and precise clinical treatment.

Introduction

Currently, cervical carcinoma has become the second most lethal malignancy among women worldwide (Bray et al., 2018), with 527,624 new cases and 265,672 tumor-specific deaths annually (Shrestha et al., 2018). For the accurate prediction of prognosis, the International Federation of Gynecology and Obstetrics put forward the cervical cancer staging standard in 2018, according to the depth, greatest dimension of stromal invasion, and the extension of tumor on adjacent regions (such as vagina and pelvis), separating cervical cancer into stages I, II, III, and IV and further substages (Balcacer et al., 2019; Salvo et al., 2020). For different stages of the tumor, the treatments are diverse. Surgical interventions, including trachelectomy and radical hysterectomy, are performed in most early cervical cancers. In addition, radiation and chemoradiation are applied in most advanced or metastatic cervical cancers (Johnson et al., 2019). Unfortunately, these treatments still result in a low response rate and poor prognosis. Therefore, it is essential to develop new predictive prognostic models for optimizing treatment strategies.

Histone proteins, a type of abundant cellular protein, are surrounded by DNA to make nucleosomes (Taylor et al., 2020). The N-terminal tail of each histone protein is the site of posttranslational modifications (PTMs) and can make contact with adjacent nucleosomes (Bannister and Kouzarides, 2011; Bartke et al., 2013; Taylor et al., 2020). There are numerous types of histone PTMS, including acetylation, methylation, ubiquitinoylation, and phosphorylation (Bannister and Kouzarides, 2011; Taylor et al., 2020; Zhang et al., 2021). They play fundamental roles in most biological processes (Bannister and Kouzarides, 2011; Taylor et al., 2020). Histone methylation functions in many levels of transcriptional regulation from the chromatin architecture to specific locus regulation and RNA processing (Greer and Shi, 2012). Histone acetylation influences myriad cellular and physiological processes, including transcription, phase separation, autophagy, mitosis, differentiation, and neural function (Shvedunova and Akhtar, 2022). Histone ubiquitinoylation works in organization of repair to preserve genomic integrity after the breaking of the DNA double strand (Uckelmann and Sixma, 2017). Histone phosphorylation provides a rapid and reversible physiological response to DNA damage, nutritional stress, or an altered metabolic state (Uckelmann and Sixma, 2017). More importantly, histone modification affects the accessibility of DNA and recruitment of DNA-binding proteins and thereby regulates gene transcription, which controls transcriptional regulation and corresponding disorders. For example, modification of histones by acetylation/deacetylation influences gene expression and is therefore related to the carcinogenic process (Audia and Campbell, 2016) (Fang et al., 2014).

Many molecular mechanisms of histone modification have also been found. A high level of histone acetylation is associated with the expression of proto-oncogenes; nonetheless, a low level of histone acetylation is linked with the silencing of tumor suppressor genes (Armenta-Castro et al., 2020). Moreover, it has been reported that the repressive expression of osteoprotegerin (OPG) genes is associated with histone modification and further causes chromatin architecture alterations and results in transcriptional repression of OPG genes. OPG plays a tumor suppressive role in tumorigenesis; therefore, histone modification-associated OPG suppression can accelerate the deterioration of cervical cancer (Lu et al., 2009). Some studies also suggest that histone acetylation is linked with the expression of the LGALS9 gene. Histone acetylation is an activator modification for the LGALS9 gene, which can be observed in active promoters. This modification can promote LGALS9 gene transcription. Galectin-9, which is encoded by the LGALS9 gene, can have a positive effect on the apoptosis of tumor cells. Histone acetylation of LGALS9 reduces deterioration in cervical cancer (Armenta-Castro et al., 2020). Relevant molecular mechanisms have been studied to some extent. However, the impact of changes in signaling pathways caused by histone modification on the progression of cervical cancer is not clear. Accordingly, we carried out the current study to illuminate the prognostic value and function of histone modification in cervical cancer patients.

In this study, we gathered information on the signaling pathways related to histone modification. According to the activation status of the pathways, we identified cervical cancer patients with activated and suppressed histone modification. Then, we sought differentially expressed genes (DEGs) to reflect the inner alteration of the tumor caused by histone modification. Among these DEGs, we screened for relevant genes that best reflected prognosis and constructed a histone modification-associated gene (HMAG) signature to predict the prognosis of patients with cervical cancer.

Methods

Data collection

We obtained the gene expression profile and corresponding clinical information of CESC patients from three independent clinical cohorts: TCGA-CESC, GSE44001, and GSE52903. We downloaded the TCGA-CESC dataset from the GDC TCGA project via the R package “TCGAbiolinks.” The transcripts per million (TPM) format of the gene expression file was chosen and then transformed to the log2 (TPM+1) type to make it comparable with the sequencing results from the microarray. Patient samples in the GSE52903 cohort were collected from Mexico City, which were HPV16-positive and fresh samples with more than 70% tumor tissues. The clinical outcome recorded in TCGA-CESC and GSE52903 was the overall survival (OS). Samples in the GSE44001 dataset were collected from South Korea, and gene expression values were detected using the GPL14951 Illumina HumanHT-12 WG-DASL V4.0 R2 expression beadchip platform, along with the prognosis information of disease-free survival (DFS). We further filtered the enrolled patients to exclude those without paired gene expression data and clinical information and those with a follow-up time of less than one month to decrease the potential bias (Table 1).

TABLE 1

TABLE 1. Clinicopathological information of the enrolled cohorts.

Removal of batch effects between cohorts

For the three enrolled cohorts, the potential nonbiological bias was eliminated to make the gene expression profiles of different cohorts more comparable. The “sva” R package was applied with its ComBat algorithms to remove the batch effects, and the gene expression profiles were all adjusted. Then, GSE44001 and GSE52903 were combined as the GEO-combined cohort. Three cohorts were detached before the removal of the batch effect (SupplementartyFigure S1A) and fused together after removal (SupplementartyFigure S1B).

Collection of the histone modification pathways

To comprehensively reveal the diverse distribution of histone modification in cervical cancer patients, we collected a total of 122 histone modification-associated signaling pathways from the Molecular Signatures Database (MSigDB)-C5 (Liberzon et al., 2011): ontology gene sets, including the process of histone-mediated phosphorylation, methylation, ubiquitination, and acetylation.

Gene set variation analysis

The activation of 122 histone modification-associated signaling pathways was assessed by the “GSVA, v.3.5” R package, the enrichment score for a specific gene set in each sample was calculated, and GSVA quantified the total gene set activation results (Hanzelmann et al., 2013). Therefore, the gene expression profiles were transferred to gene set activation profiles, including the activated scores of 122 signaling pathways for each cohort.

Revealing the differentially expressed genes and pathways

To identify the diverse downstream altered biological processes, we carried out the differentially expressed genes (DEGs) among histone modification-activated and -suppressed subgroups via the “limma” package, with the preset threshold value of p < 0.05 and |log2fc| > 0.4. Pathway enrichment was performed by the “clusterProfiler 4.0” R package (Wu et al., 2021), with the employment of 50 HALLMARK pathways and KEGG pathways. In addition, Metascape (http://metascape.org) (Zhou et al., 2019) was also used to annotate the DEGs to reveal the activated pathways.

Construction and calculation of the HMAG signature

The prognostic value of the aforementioned identified DEGs was further assessed by univariate Cox regression analysis. Gene expression was first separated into high and low groups and then entered into the univariate Cox regression analysis. The prognostic genes in both the TCGA-CESC cohort and GEO-combined cohort were selected for a subsequent least absolute shrinkage and selection operator (LASSO) regression analysis, which was conducted using the “glmnet” package. The selected genes were used to calculate the risk score by adding the gene expression multiplied by the corresponding coefficient and named the HMAG score. The HMAG signature was trained in the TCGA-CESC cohort and validated in the GSE44001 and GSE52903 cohorts. The HMAG score was calculated for each patient and then applied for the subsequent analysis.

Multivariate analysis and establishment of the nomogram

The prognostic value of the HMAG score in each cohort was assessed by the K-M plot and receiver operating characteristic (ROC) curve. Multivariate Cox regression analysis was also performed to adjust the potential impact of other clinical features and is presented with a forest plot. Furthermore, we employed the “regplot” package to establish the clinical prognostic nomogram, which could provide a quantitative method for the individualized prediction of cervical cancer. Factors that emerged from the multivariate Cox regression analysis were enrolled for the establishment of the nomogram. C-index curve, calibration curve, and decision curve analyses were all performed to validate the clinical usefulness and accuracy of the nomogram via the “rms” and “rmda” packages.

Revealing activated pathways and potential therapeutic drugs

We assessed the activated pathways by fast gene set enrichment analysis (fgsea). First, GSEA is performed with the ranking of input molecular readouts, and then the pathway enrichment score is calculated by walking down the list of features, which means that if a feature is felled into the target pathway, the running-sum statistic increases; otherwise, it decreases. The final score is the maximum deviation from zero encountered in the random walk and normalized by computing the z-score of the estimate compared to a null distribution obtained from a random permutation. Moreover, the highly expressed genes in high-HMAG and low-HMAG groups were also annotated to reveal the activated biological process by clusterProfiler 4.0 and visualized via the treeplot model.

For the prediction of potential therapeutic drugs, we used the data from GDSC2016 (https://www.cancerrxgene.org/) and visualized them via “MOVICS” (Lu et al., 2020). Ridge regression analysis was performed to predict the potential response result to the chemotherapy drug of each patient and represented by the estimated inhibitory concentration (IC₅₀); the lower the IC₅₀ was, the higher the effectiveness of the drug treatment was (Geeleher et al., 2014). Moreover, we also searched for potential therapeutic new drugs through the Gene Set Cancer Analysis (GSCA) online website (http://bioinfo.life.hust.edu.cn/GSCA).

Comparison of the HMAG signature with published signatures

To assess the prognostic value of the HMAG signature carried out in the current study, we searched published articles to collect signatures from other studies. Finally, we used four signatures, namely, immune-associated genes (Yu et al., 2021), DNA damage repair-associated genes (Zhou et al., 2022), ferroptosis-related genes (Qi et al., 2021), and autophagy-related genes (Chen et al., 2020). The risk score of each signature was calculated in the TCGA-CESC cohort, and the ROC value and C-index were used to evaluate the prognostic value.

Statistics

The log-rank test was used to compare the survival outcome between two groups, Student’s t-test was used to compare the distribution between two groups, and Fisher’s exact test was performed to distinguish the difference in categorical data. All statistical analyses were performed using R (Version: 4.1.2). A two-tailed p-value < 0.05 was recognized as statistically significant.

Results

Characterizing the diversified histone modification-activated status and altered signaling pathways

A total of 122 histone modification-associated pathways were obtained from MSigDB-C5, and the activated score was generated based on the GSVA quantification. Then, the 248 cervical cancer patients from TCGA-CESC cohort were divided into three subtypes based on the distanceMatrix function of the “ClassDiscovery” package, with the preset parameters of “euclidean” and “ward.D.” We observed diversified histone modification activation in C1, C2, and C3. The patients in C2 contained the most activated status of histone modification, while the C1 patients contained the suppressed histone modification (Figure 1A). Then, we compared the DEGs between C2 and C1 (Figure 1B), with a preset threshold value of p < 0.05 and |log2fc| > 0.4. A total of 2040 genes were upregulated in C2, and the other 200 genes were upregulated in C1. We combined the DEGs between C1 and C2, annotated the enriched signaling pathways, and observed the altered biological processes of the cell cycle, cellular response to DNA damage stimulus, chromatin organization, and DNA metabolic process (Figure 1C). In addition, we also validated the DEGs in the HALLMARK and KEGG pathways. The DEGs targeted the G2/M checkpoint, mitotic spindle, and E2F targets (Figure 1D) and were also linked with proteoglycans in cancer, axon guidance, nucleocytoplasmic transport, and cell cycle pathways in KEGG pathways (Figure 1E).

FIGURE 1

FIGURE 1. Identification of the histone modification-altered genes and pathways. (A) Heatmap showing the activation status and clinical features in cervical cancer patients. (B) Differentially expressed genes (DEGs) between the histone modification-activated and -suppressed subgroups. (C) Enriched signaling pathways of the DEGs using Metascape. (D) Enriched signaling pathways of the DEGs by HALLMARK pathways. (E) Enriched signaling pathways of the DEGs by KEGG pathways.

Distinguishing the prognostic genes and constructing the signature

Based on the 2240 DEGs, we assessed and selected the prognostic DEGs in TCGA-CESC cohort (Figure 2A) and GEO-combined cohort (Figure 2B). Genes with an HR higher than 1 and a p-value less than 0.05 were regarded as risky genes, while genes with an HR less than 1 and a p-value less than 0.05 were regarded as protective genes. Ninety-six risk genes from TCGA-CESC cohort and 127 risk genes from the GEO-combined cohort merged out six common risk genes; 90 protective genes from TCGA-CESC cohort and 162 protective genes from the GEO-combined cohort merged out 17 common protective genes (Figure 2C). Therefore, 23 genes were finally enrolled for LASSO analysis, and nine genes were identified under the best optimal lambda value of 0.055, leading to the prognostic genes SOC21, HLF, FGFR2, MYLIP, ZDHHC11, FASN, PDK1, MYO10, and TNFRSF12A (Figures 2D,E). The HMAG signature was calculated using the following formula: HMAG score = 0.115713151 * TNFRSF12A + 0.032286676 * MYO10 + 0.300717617 * PDK1 + 0.131853347 * FASN + -0.015900246 * MYLIP - 0.033923872 * FGFR2 - 0.099329649 * HLF - 0.034241118 * SOX21 - 0.042765486 * ZDHHC11.

FIGURE 2

FIGURE 2. Construction of the histone modification-associated gene signature. (A) Selection of the prognostic genes in TCGA-CESC cohort; (B) selection of the prognostic genes in the GEO-combined cohort; (C) Venn plot showing the risk genes and protective genes in both TCGA-CESC and GEO-combined cohorts; (D) optimal tuning parameter (lambda) in the LASSO analysis selected with 10-fold cross-validation and one standard error rule; and (E) LASSO coefficient profiles of the 23 candidate genes.

Prognostic effectiveness of the HMAG signature

We calculated the HMAG score of each patient based on the aforementioned formula. The overall distributions of the risk score, survival status, and gene expression profiles of the nine-gene signature in TCGA-CESC cohort (left), the GSE44001 cohort (middle), and the GSE52903 cohort (right) are shown in Figure 3A. Patients in each cohort were separated into HMAG-L and HMAG-H based on the median HMAG score, and the prognostic value was assessed by Cox regression analysis. Cervical patients with HMAG-H in TCGA-CESC cohort showed a 2.68-fold change of death risk, with the 95% CI from 1.533 to 4.671 (p < 0.001), as well as the increased death risk of HMAG-H in the GSE44001 cohort (HR: 2.83, 95% CI: 1.370–5.849, p = 0.005) and GSE44001 cohort (HR: 4.59, 95% CI: 1.658–12.697, p = 0.003) (Figure 3B).

FIGURE 3

FIGURE 3. Risk map and K-M plot for TCGA-CESC, GSE44001, and GSE52903 cohorts. (A) Risk score distribution (upper), clinical outcome distribution (middle), and gene expression distribution (bottom) for TCGA-CESC, GSE44001, and GSE52903 cohorts. (B) K-M plot of the HMAG-L score and HMAG-H score subgroups for TCGA-CESC, GSE44001, and GSE52903 cohorts.

The prognostic accuracy was evaluated by ROC curves. We observed the preferable AUC values of the HMAG signature in TCGA-CESC cohort (1-year: 0.719, 3-year: 0.741, and 5-year: 0.731) and GSE44001 cohort (1-year: 0.850, 3-year: 0.781, and 5-year: 0.755) (Figure 4A). In addition, we also conducted multivariate Cox regression analysis to adjust for the impact of clinicopathological features. We observed that the HMAG signature (HR: 3.509, 95% CI: 1.869–6.590, p < 0.001) and tumor status (tumor vs. tumor free: HR: 37.094, 95% CI: 15.073–91.290, p < 0.001) were independent prognostic factors for cervical cancer patients in TCGA-CESC cohort but age, menopausal status, tumor stage, and tumor grade were not (Figure 4B). In the GSE44001 cohort, the HMAG signature also showed an independent risk factor after adjusting for tumor size and tumor stage (HR: 2.43, 95% CI: 1.166–5.070, p = 0.0178, Figure 4C). Similar results were also shown in the GSE52903 cohort (HR: 3.02, 95% CI: 1.021–8.920, p = 0.0459, Figure 4D).

FIGURE 4

FIGURE 4. Evaluation of the prognostic value of the HMAG signature. (A) 1-year, 3-year, and 5-year AUC values in TCGA-CESC, GSE44001, and GSE52903 cohorts; (B) forest plot showing the prognostic value of the HMAG signature after adjusting for clinical features in TCGA-CESC cohort; (C) forest plot showing the prognostic value of the HMAG signature after adjusting for clinical features in the GSE44001 cohort; and (D) forest plot showing the prognostic value of the HMAG signature after adjusting for clinical features in the GSE52903 cohort.

Prognostic nomogram established by the HMAG score and tumor status

With the aforementioned results, we established a prognostic nomogram model containing the factors HMAG score and tumor status (Figure 5A). For a specific patient, the HMAG score and tumor status correspond to the point, and the summary of two points is the total point, with a straight line from the total point site to the bottom line of 1-year, 3-year, and 5-year death risks, indicating the risk of death. With a cumulative incidence plot, we visualized the estimated probability of the death event prior to a specified time. The patients were separated into a low-point subgroup and a high-point subgroup, and we observed significantly diverse cumulative events in the two groups (HR: 2.68, 95% CI: 1.533–4.671, p < 0.001, Figure 5B). The AUC values changed over the follow-up time and presented a better prognostic value of the nomogram than regardless of age, grade, or stage (Figure 5C). The overall C-index of the nomogram showed a prognostic value as high as 0.890, while the C-index for age was only 0.562, for grade was only 0.542, and for stage was only 0.614 (Figure 5D). A p-value of 0.854 calculated by the Hosmer‒Lemeshow analysis in the calibration plot indicated that the prediction performance of this nomogram might be equivalent to an ideal predictive model (Figure 5E). DCA was performed to demonstrate a high clinical net benefit that was almost over the entire threshold probability of the nomogram model compared with other features (Figure 5F).

FIGURE 5

FIGURE 5. Prognostic nomogram constructed by the HMAG signature and tumor status. (A) Establishment of a nomogram combining the tumor status and HMAG signature; (B) K-M plot showing the diverse cumulative death event in low points and high points subgroups; (C) time ROC lines showing the prognostic value of nomogram, age, grade, and tumor stage; (D) bar plots showing the C-index of nomogram, age, grade, and tumor stage; (E) calibration plot for the nomogram. The dashed line represents the ideal nomogram, the solid line represents our nomogram, and a p-value of 0.854 indicates that the nomogram-predicted probability is very close to the actual death events; (F) DCA showed that our nomogram had the greatest net benefit compared with the single factor of age, stage, or grade.

Revealing activated pathways and potential therapeutic drugs

The fgsea analysis was performed on TCGA-CESC cohort based on the tumor HALLMARK pathways. We revealed that patients with high HMAG scores contained the activated tumor pathways of angiogenesis, G2/M checkpoint, MYC targets, and hypoxia, while patients with low HMAG scores contained the activated tumor pathways of interferon gamma response, interferon alpha response, and IL2/STAT5 signaling (Figure 6A). The DEGs between the high and low HMAG score subgroups also indicated that the pathways of fatty acid peptidyl-tyrosine biosynthetic, positive regulation angiogenesis vasculature, decreased oxygen levels, extracellular matrix encapsulating disassembly, and endoderm formation differentiation terms were activated in high HMAG score patients (Figure 6B), and the pathways of apoptotic antigen presentation processing, lymphocyte mononuclear differentiation, arachidonic eicosanoid unsaturated fatty, positive cell‒cell activation adhesion, and cellular hormone alcohol compound terms were activated in low HMAG score patients (Figure 6C).

FIGURE 6

FIGURE 6. Revealing HMAG signature-regulated biological pathways. (A) fgsea showing the diverse activation of tumor biological pathways in the high-HMAG-score group and the low-HMAG-score group; (B) hierarchical clustering of enriched terms in HMAG upregulated genes; and (C) hierarchical clustering of enriched terms in HMAG downregulated genes.

We identified suitable therapeutic drugs based on the GDSC2016 database and revealed that patients with high HMAG scores were more suitable for the treatment of CHIR-99021, embelin, FTI-277, JNK-9L, JQ12, midostaurin, PF-562271, pyrimethamine, and thapsigargin, and patients with low HMAG scores were more suitable for the treatment of BMS-536924, CP466722, crizotinib, PHA-665752, rapamycin, and TAE684 (all p < 0.05, Figure 7A). Similar results were also validated in the GEO-combined cohort (all p < 0.05, Figure 7B). We also searched for new drugs from the GSCA dataset to treat the specific alteration of a single gene, and a total of 41 components were identified that might be functional in the clinical treatment of cervical cancer patients (Table 2).

FIGURE 7

FIGURE 7. Prediction of therapeutic drugs for the HMAG low- and high-risk groups. (A) Potential therapeutic chemo drugs identified in TCGA-CESC cohort and (B) potential therapeutic chemo drugs validated in the GEO cohort.

TABLE 2

TABLE 2. Predicting potential chemo drugs via the Cancer Therapeutics Response Portal dataset.

The HMAG signature showed better prognostic value than other published signatures

We collected the published signatures of the genes and corresponding indexes from four studies, referring to the genes in immunity, DNA damage repair, ferroptosis, and autophagy. The risk score of each signature was calculated along with the HMAG score using the formula. We observed that the HMAG signature showed more excellent prognostic value than any other signature based on the AUC value (HMAG: 0.711, Yu et al.: 0.632, Zhou et al.: 0.572, Qi et al.: 0.652, and Chen et al.: 0.661, Figure 8A) and C-index (HMAG: 0.757, Yu et al.: 0.723, Zhou et al.: 0.626, Qi et al.: 0.637, and Chen et al.: 0.718, Figure 8B).

FIGURE 8

FIGURE 8. Comparison of the prognostic value between HMAG and proposed signatures. (A) ROC curves showing the prognostic value of five signatures and (B) C-index showing the prognostic value of five signatures.

Discussion

Epigenetic modifications are reversible and do not change the DNA sequence of genetic material but can alter how our body reads the DNA sequence, including the basic forms of DNA methylation, posttranslational modifications on histones and noncoding RNAs. However, regulation of the epigenome can distinctly lead to gene malfunction, regulating cell differentiation, proliferation, and even apoptosis and causing disordered cell growth and tumorigenesis. The abnormal control of posttranslational modification-related enzymes, including histone methylase, demethylase, acetylase, and acetyltransferase, acts as a pivotal risk factor for tumors, and these epigenetic alterations may lead to the reprogramming of genomes, activation of oncogenes, or silencing of tumor suppressors (Berry and Janknecht, 2013). Epigenetic regulators can impact the level of histone modification to the enhancer activity of genes via histone methyltransferases or acetylases (Xiao et al., 2022). Histone deacetylases (HDACs) are essential for maintaining the balance of cell processes by altering histone deacetylation. The abnormal expression of HDASs is tightly linked with several cancers, and their inhibitors have been applied in the clinic to treat several cancers (Patel et al., 2022). Moreover, histone modification impacts the process of epithelial–mesenchymal transition (EMT), and histone acetylation can modulate the acetylation levels of distinct histones at the promoters of EMT-related markers, EMT-inducing transcription factors, and EMT-related long noncoding RNAs to control EMT (Kong et al., 2022).

Histone modification function also emerges in the carcinogenesis of cervical cancer. Higher staining of H3K9ac indicates low grading, negative N-status, and low T-status in cervical cancer. Moreover, the increased expression of H3K4me3 in the cytoplasm was observed to be associated with advanced T stage and unfavorable prognosis in cervical cancer patients (Beyer et al., 2017). Galectin-9, which is encoded by LGALS9, is evidently detected in normal epithelium and endocervical glands but not in cervical intraepithelial neoplasia and cervical squamous cell carcinoma, which indicates that decreased Galectin-9 is a biomarker for the malignant potential of cervical cancer (Liang et al., 2008). Armenta-Castro et al. (2020) reported that histone H3K9 and H3K14 acetylation at the promoter of LGALS9 genes is tightly associated with the protein level of Galectin-9, and the LGALS9 gene presented higher levels of histone acetylation in normal cervical cells than in cancer cells. Sun et al. (2022) also reported that the inhibition of HDACs can activate mitophagy by mediating Parkin acetylation and leading to the inhibition of cervical cancer cell proliferation. Valproic acid/sodium valproate (VPA), a well-known antiepileptic agent, inhibits histone deacetylases, induces histone hyperacetylation, promotes DNA demethylation, and affects the histone methylation status in some cell models. Rocha et al. demonstrated that VPA promotes the abundance of H3K4me2/me3 and increases methyltransferase KMT2D gene expression in HeLa cells. Meanwhile, VPA can also induce the hypomethylation of H3K9me2 and concomitant with the increased gene expression of KDM3A (Rocha et al., 2022). Based on the aforementioned evidence, we sensed the important role of histone modification in cervical cancer tumorigenesis; therefore, we evaluated the histone modification pathway activation status in cervical cancer and generated a prognostic model to predict clinical prognosis.

In the current study, we enrolled a total of 594 cervical cancer patients. We first evaluated the histone modification activation status of 248 patients from TCGA-CESC cohort and identified the activated and suppressed subgroups. Further gene expression variation analysis filtered out 2240 genes and merged them with the results of prognostic analysis. Twenty-three histone modification-associated prognostic genes were enrolled in the LASSO Cox regression analysis. Eventually, we constructed the HMAG signature to predict the prognosis for cervical cancer patients, with nine genes: TNFRSF12A, MYO10, PDK1, FASN, MYLIP, FGFR2, HLF, SOX21, and ZDHHC11. The HMAG signature showed a preferable prognostic value in TCGA-CESC, GSE44001, and GSE52903 cohorts, and the HMAG signature is an independent prognostic factor after adjusting for the influence of other clinicopathological features. In addition, we established a prognostic nomogram combining the HMAG score and tumor status. The nomogram showed better prognostic value than did age, tumor grade, or tumor stage, with a C-index as high as 0.890. Currently, precise treatment raises concerns in the clinic, and we predicted suitable chemotherapy drugs for patients with high or low HMAG scores. The index of SOX21 in the HMAG signature formula is 0.034241118, which means that the higher expression of SOX21 is associated with a favorable prognosis. The elevated expression of histone H3R26-methylase CARM1 can increase the level of Sox21, and CARM1 inhibition can decrease Sox21 levels, which highlights the importance of epigenetic regulation of Sox21 (Goolam et al., 2016). Caglayan et al. (2013) demonstrated that the exogenous expression of Sox21 in tumor cells resulted in a significant decrease in the tumor size. It seems that Sxo21 appears to inhibit stem-like cell properties and initiate the aberrant differentiation of glioma cells. Kurtsdotter et al. also reported that the high levels of SOX5/6/21 in human primary glioblastoma cells enabled the expression of CDK inhibitors and decreased p53 protein turnover, which blocked their tumorigenic capacity through cellular senescence and apoptosis. FASN is a multienzyme protein that serves as the key regulator in lipid metabolism, especially fatty acid synthesis. KDM5C is a histone H3K4-specific demethylase, and the overexpression of KDM5C led to the reduction of H3K4me3 on the promoter and the corresponding downregulation of FASN expression to inhibit FASN-mediated lipid metabolism (Zhang et al., 2020). Du et al. (2022) illustrated that FASN can promote the lymph node metastasis of cervical cancer via cholesterol reprogramming and lymph angiogenesis. The mechanisms of how histone modification influences the tumorigenesis of cervical cancer are complicated and need a further in-depth study.

Conclusion

In summary, we comprehensively evaluated the histone modification status in cervical cancer patients and revealed histone modification-associated prognostic genes to construct the HMAG signature, aiming to provide new insights into prognosis prediction and precise clinical treatment.

Data availability statement

The original contributions presented in the study are included in the article/Supplementary Material; further inquiries can be directed to the corresponding authors.

Author contributions

All authors listed have made a substantial, direct, and intellectual contribution to the work and approved it for publication.

Funding

This work was supported by the 2020–2022 Clinical Research Fund of The Seventh Affiliated Hospital, Sun Yat-sen University.

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.

Supplementary material

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

Supplementary Figure S1 | Data pretreatment of TCGA-CESC, GSE44001 and GSE52903 cohorts.

References

Armenta-Castro, E., Reyes-Vallejo, T., Maximo-Sanchez, D., Herrera-Camacho, I., Lopez-Lopez, G., and Reyes-Carmona, S. (2020). Histone H3K9 and H3K14 acetylation at the promoter of the LGALS9 gene is associated with mRNA levels in cervical cancer cells. FEBS Open Bio 10 (11), 2305–2315. doi:10.1002/2211-5463.12973

PubMed Abstract | CrossRef Full Text | Google Scholar

Audia, J. E., and Campbell, R. M. (2016). Histone modifications and cancer. Cold Spring Harb. Perspect. Biol. 8 (4), a019521. doi:10.1101/cshperspect.a019521

PubMed Abstract | CrossRef Full Text | Google Scholar

Balcacer, P., Shergill, A., and Litkouhi, B. (2019). MRI of cervical cancer with a surgical perspective: Staging, prognostic implications and pitfalls. Abdom. Radiol. 44 (7), 2557–2571. doi:10.1007/s00261-019-01984-7

PubMed Abstract | CrossRef Full Text | Google Scholar

Bannister, A. J., and Kouzarides, T. (2011). Regulation of chromatin by histone modifications. Cell Res. 21 (3), 381–395. doi:10.1038/cr.2011.22

PubMed Abstract | CrossRef Full Text | Google Scholar

Bartke, T., Borgel, J., and DiMaggio, P. A. (2013). Proteomics in epigenetics: New perspectives for cancer research. Brief. Funct. Genomics 12 (3), 205–218. doi:10.1093/bfgp/elt002

PubMed Abstract | CrossRef Full Text | Google Scholar

Berry, W. L., and Janknecht, R. (2013). KDM4/JMJD2 histone demethylases: Epigenetic regulators in cancer cells. Cancer Res. 73 (10), 2936–2942. doi:10.1158/0008-5472.CAN-12-4300

PubMed Abstract | CrossRef Full Text | Google Scholar

Beyer, S., Zhu, J., Mayr, D., Kuhn, C., Schulze, S., and Hofmann, S. (2017). Histone H3 acetyl K9 and histone H3 tri methyl K4 as prognostic markers for patients with cervical cancer. Int. J. Mol. Sci. 18 (3), E477. doi:10.3390/ijms18030477

PubMed Abstract | CrossRef Full Text | Google Scholar

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

PubMed Abstract | CrossRef Full Text | Google Scholar

Caglayan, D., Lundin, E., Kastemar, M., Westermark, B., and Ferletta, M. (2013). Sox21 inhibits glioma progression in vivo by forming complexes with Sox2 and stimulating aberrant differentiation. Int. J. Cancer 133 (6), 1345–1356. doi:10.1002/ijc.28147

PubMed Abstract | CrossRef Full Text | Google Scholar

Chen, H., Deng, Q., Wang, W., Tao, H., and Gao, Y. (2020). Identification of an autophagy-related gene signature for survival prediction in patients with cervical cancer. J. Ovarian Res. 13 (1), 131. doi:10.1186/s13048-020-00730-8

PubMed Abstract | CrossRef Full Text | Google Scholar

Du, Q., Liu, P., Zhang, C., Liu, T., Wang, W., Shang, C., et al. (2022). FASN promotes lymph node metastasis in cervical cancer via cholesterol reprogramming and lymphangiogenesis. Cell Death Dis. 13 (5), 488. doi:10.1038/s41419-022-04926-2

PubMed Abstract | CrossRef Full Text | Google Scholar

Fang, J., Zhang, H., and Jin, S. (2014). Epigenetics and cervical cancer: From pathogenesis to therapy. Tumour Biol. 35 (6), 5083–5093. doi:10.1007/s13277-014-1737-z

PubMed Abstract | CrossRef Full Text | Google Scholar

Geeleher, P., Cox, N. J., and Huang, R. S. (2014). Clinical drug response can be predicted using baseline gene expression levels and in vitro drug sensitivity in cell lines. Genome Biol. 15 (3), R47. doi:10.1186/gb-2014-15-3-r47

PubMed Abstract | CrossRef Full Text | Google Scholar

Goolam, M., Scialdone, A., Graham, S. J. L., Macaulay, I. C., Jedrusik, A., Hupalowska, A., et al. (2016). Heterogeneity in Oct4 and Sox2 targets biases cell fate in 4-cell mouse embryos. Cell 165 (1), 61–74. doi:10.1016/j.cell.2016.01.047

PubMed Abstract | CrossRef Full Text | Google Scholar

Greer, E. L., and Shi, Y. (2012). Histone methylation: A dynamic mark in health, disease and inheritance. Nat. Rev. Genet. 13 (5), 343–357. doi:10.1038/nrg3173

PubMed Abstract | CrossRef Full Text | Google Scholar

Hanzelmann, S., Castelo, R., and Guinney, J. (2013). Gsva: Gene set variation analysis for microarray and RNA-seq data. BMC Bioinforma. 14, 7. doi:10.1186/1471-2105-14-7

PubMed Abstract | CrossRef Full Text | Google Scholar

Johnson, C. A., James, D., Marzan, A., and Armaos, M. (2019). Cervical cancer: An overview of pathophysiology and management. Semin. Oncol. Nurs. 35 (2), 166–174. doi:10.1016/j.soncn.2019.02.003

PubMed Abstract | CrossRef Full Text | Google Scholar

Kong, F., Ma, L., Wang, X., You, H., Zheng, K., and Tang, R. (2022). Regulation of epithelial-mesenchymal transition by protein lysine acetylation. Cell Commun. Signal. 20 (1), 57. doi:10.1186/s12964-022-00870-y

PubMed Abstract | CrossRef Full Text | Google Scholar

Liang, M., Ueno, M., Oomizu, S., Arikawa, T., Shinonaga, R., and Zhang, S. (2008). Galectin-9 expression links to malignant potential of cervical squamous cell carcinoma. J. Cancer Res. Clin. Oncol. 134 (8), 899–907. doi:10.1007/s00432-008-0352-z

PubMed Abstract | CrossRef Full Text | Google Scholar

Liberzon, A., Subramanian, A., Pinchback, R., Thorvaldsdottir, H., Tamayo, P., and Mesirov, J. P. (2011). Molecular signatures database (MSigDB) 3.0. Bioinformatics 27 (12), 1739–1740. doi:10.1093/bioinformatics/btr260

PubMed Abstract | CrossRef Full Text | Google Scholar

Lu, T. Y., Kao, C. F., Lin, C. T., Huang, D. Y., Chiu, C. Y., and Huang, Y. S. (2009). DNA methylation and histone modification regulate silencing of OPG during tumor progression. J. Cell. Biochem. 108 (1), 315–325. doi:10.1002/jcb.22256

PubMed Abstract | CrossRef Full Text | Google Scholar

Lu, X., Meng, J., Zhou, Y., Jiang, L., and Yan, F. (2020). Movics: an R package for multi-omics integration and visualization in cancer subtyping. Bioinformatics 36, 5539–5541. doi:10.1093/bioinformatics/btaa1018

CrossRef Full Text | Google Scholar

Patel, P., Wahan, S. K., Vishakha, S., Kurmi, B. D., Gupta, G. D., and Rajak, H. (2022). Recent progress in histone deacetylase (HDAC) 1 inhibitors as anticancer agent. Curr. Cancer Drug Targets 22. doi:10.2174/1568009622666220624090542

CrossRef Full Text | Google Scholar

Qi, X., Fu, Y., Sheng, J., Zhang, M., Zhang, M., and Wang, Y. (2021). A novel ferroptosis-related gene signature for predicting outcomes in cervical cancer. Bioengineered 12 (1), 1813–1825. doi:10.1080/21655979.2021.1925003

PubMed Abstract | CrossRef Full Text | Google Scholar

Rocha, M. A., de Campos Vidal, B., and Mello, M. L. S. (2022). Sodium valproate modulates the methylation status of lysine residues 4, 9 and 27 in histone H3 of HeLa cells. Curr. Mol. Pharmacol. 15. doi:10.2174/1874467215666220316110405

CrossRef Full Text | Google Scholar

Salvo, G., Odetto, D., Pareja, R., Frumovitz, M., and Ramirez, P. T. (2020). Revised 2018 international federation of Gynecology and Obstetrics (figo) cervical cancer staging: A review of gaps and questions that remain. Int. J. Gynecol. Cancer 30 (6), 873–878. doi:10.1136/ijgc-2020-001257

PubMed Abstract | CrossRef Full Text | Google Scholar

Shrestha, A. D., Neupane, D., Vedsted, P., and Kallestrup, P. (2018). Cervical cancer prevalence, incidence and mortality in low and middle income countries: A systematic review. Asian pac. J. Cancer Prev. 19 (2), 319–324. doi:10.22034/APJCP.2018.19.2.319

PubMed Abstract | CrossRef Full Text | Google Scholar

Shvedunova, M., and Akhtar, A. (2022). Modulation of cellular processes by histone and non-histone protein acetylation. Nat. Rev. Mol. Cell Biol. 23 (5), 329–349. doi:10.1038/s41580-021-00441-y

PubMed Abstract | CrossRef Full Text | Google Scholar

Sun, X., Shu, Y., Ye, G., Wu, C., Xu, M., and Gao, R. (2022). Histone deacetylase inhibitors inhibit cervical cancer growth through Parkin acetylation-mediated mitophagy. Acta Pharm. Sin. B 12 (2), 838–852. doi:10.1016/j.apsb.2021.07.003

PubMed Abstract | CrossRef Full Text | Google Scholar

Taylor, E. M., Koss, B., Davis, L. E., and Tackett, A. J. (2020). Histone modifications as biomarkers for immunotherapy. Methods Mol. Biol. 2055, 213–228. doi:10.1007/978-1-4939-9773-2_10

PubMed Abstract | CrossRef Full Text | Google Scholar

Uckelmann, M., and Sixma, T. K. (2017). Histone ubiquitination in the DNA damage response. DNA Repair (Amst) 56, 92–101. doi:10.1016/j.dnarep.2017.06.011

PubMed Abstract | CrossRef Full Text | Google Scholar

Wu, T., Hu, E., Xu, S., Chen, M., Guo, P., Dai, Z., et al. (2021). clusterProfiler 4.0: A universal enrichment tool for interpreting omics data. Innovation. 2 (3), 100141. doi:10.1016/j.xinn.2021.100141

PubMed Abstract | CrossRef Full Text | Google Scholar

Xiao, Q., Xiao, Y., Li, L. Y., Chen, M. K., and Wu, M. (2022). Multifaceted regulation of enhancers in cancer. Biochim. Biophys. Acta. Gene Regul. Mech. 1865 (6), 194839. doi:10.1016/j.bbagrm.2022.194839

PubMed Abstract | CrossRef Full Text | Google Scholar

Yu, S., Li, X., Zhang, J., and Wu, S. (2021). Development of a novel immune infiltration-based gene signature to predict prognosis and immunotherapy response of patients with cervical cancer. Front. Immunol. 12, 709493. doi:10.3389/fimmu.2021.709493

PubMed Abstract | CrossRef Full Text | Google Scholar

Zhang, B., Zhou, B. H., Xiao, M., Li, H., Guo, L., and Wang, M. X. (2020). KDM5C represses FASN-mediated lipid metabolism to exert tumor suppressor activity in intrahepatic cholangiocarcinoma. Front. Oncol. 10, 1025. doi:10.3389/fonc.2020.01025

PubMed Abstract | CrossRef Full Text | Google Scholar

Zhang, Y., Sun, Z., Jia, J., Du, T., Zhang, N., and Tang, Y. (2021). Overview of histone modification. Adv. Exp. Med. Biol. 1283, 1–16. doi:10.1007/978-981-15-8104-5_1

PubMed Abstract | CrossRef Full Text | Google Scholar

Zhou, H., Wu, L., Yu, L., Yang, Y., Kong, L., and Liu, S. (2022). Identify a DNA damage repair gene signature for predicting prognosis and immunotherapy response in cervical squamous cell carcinoma. J. Oncol. 2022, 8736575. doi:10.1155/2022/8736575

PubMed Abstract | CrossRef Full Text | Google Scholar

Zhou, Y., Zhou, B., Pache, L., Chang, M., Khodabakhshi, A. H., and Tanaseichuk, O. (2019). Metascape provides a biologist-oriented resource for the analysis of systems-level datasets. Nat. Commun. 10 (1), 1523. doi:10.1038/s41467-019-09234-6

PubMed Abstract | CrossRef Full Text | Google Scholar