Identification of HIST1H2BH as the hub gene associated with multiple myeloma using integrated bioinformatics analysis

ABSTRACT

Objectives

Identifying the specific biomarkers and molecular signatures of MM might provide novel evidence for MM prognosis and targeted therapy.

Methods

Bioinformatic analyses were performed through GEO and TCGA datasets. The differential expression of HIST1H2BH in MM sample was validated by the qRT-PCR. And the CCK-8 assay was performed to detect the proliferation activity of HIST1H2BH on MM cell lines.

Results

A total of 793 DEGs were identified between bone marrow plasma cells from newly diagnosed myeloma and normal donors in GSE6477. Among them, four vital genes (HIST1H2AC, HIST1H2BH, CCND1 and TCF7L2) modeling were constructed. The increased HIST1H2BH expression was correlated with worse survival of MM based on TCGA datasets. The transcriptional expression of HIST1H2BH was significantly up-regulated in primary MM patients. And knockdown HIST1H2BH decreased the proliferation of MM cell lines.

Conclusions

We have identified up-regulated HIST1H2BH in MM patients associated with poor prognosis using integrated bioinformatical methods.

KEYWORDS:

• Multiple myeloma
• bioinformatics analysis
• independent prognostic factor‌
• differentially expressed gene
• hub genes
• epigenetic regulation
• H2B protein family
• histone modification

Introduction

As an incurable, biologically heterogeneous disease of the plasma cells, MM is one of the hematologic cancers that account for approximately 10% of all hematologic malignancies and is the second most common hematologic malignancy [1–5]. The incidence of MM in China is 0.9/100,000 per year, much lower than that in the United States (6.56/100,000 per year) [6–10]. Overproduction of immunoglobulins or immunoglobulin chains that are unfunctional due to abnormal monoclonal plasma cell growth in bone marrow. Bone lesions, anemia, hypercalcemia and other symptoms are caused by the accumulation of intact immunoglobulins and aberrant monoclonal plasma cells in the bone marrow [11–16]. Recent studies have found that MM cells induce changes in the bone marrow micro-environment through complex genetic and epigenetic abnormalities, affecting the tumor's clinical course and directly impacting patient prognosis [17]. Therefore, it is crucial to explore the pathogenesis of MM at the molecular level to identify new tumor markers and therapeutic targets and improve early diagnosis, treatment outcomes and prognosis.

There is evidence that changes in histone modifications can specifically activate or suppress the transcriptional activity of specific genes and that histone modification-related proteins are closely associated with MM epigenetic regulation [18, 19]. H2B protein plays an essential role in histone modification [20–22]. Previous studies have demonstrated that mutations in H2B may represent a class of oncogenic drivers. HIST1H2BH, as a critical member of the H2B histone family, has also been shown to serve as a predictive factor for cervical cancer prognosis [23].

In the present study, bioinformatics methods were used to investigate the differential gene expression profiles of MM samples in both the Gene Expression Omnibus (GEO, https://www.ncbi.nlm.nih.gov/geo/) database and The Cancer Genome Atlas Program (TCGA, https://portal.gdc.cancer.gov/). A critical hub gene analysis was then conducted and HIST1H2BH level was found to be used as independent prognostic factors for MM. Additionally, higher expression levels of MM patients were more likely to have a reduced survival time. Furthermore, the expression level of HIST1H2BH was then validated in primary MM patients and its’ function was explored in vitro experiments. Our research may be helpful to identify new molecular signature and therapeutic target of MM patients.

Methods

Data collection

GSE6477, related to MM dataset, was downloaded from the GEO database. Moreover, this dataset, with data platform GPL96, has CD138-selected bone marrow plasma cells from 69 new MM and 15 normal samples, consisting of array data corresponding to Homo sapiens.

Differential expression analysis

Differentially expressed genes (DEGs) between 69 new MM and 15 healthy samples of the GSE6477 dataset was identified by the ‘limma’ R package (https://www.bioconductor.org/packages/ release/bioc/html/limma.html) [24] based on the GEO database [25]. The selection criteria for DEGs were: adjusted P-value (adj. P) < 0.05 and |log FC| > 1, indicating that the gene is dis-regulated.

GO and KEGG pathway enrichment analyses

An online website (bioinformatics.com.cn) was used to analyze GO and KEGG pathway enrichment. A total of 793 DEGs genes were analyzed to detect the GO terms that displayed enrichment in three categories (biological processes, cellular components and molecular functions) and the KEGG pathways.

Protein−protein interaction analysis

The online website STRING (https://string-db.org/) was used to analyze PPI [26]. Firstly, Protein IDs were translated into gene symbols, and those PPIs without corresponding gene names were excluded. Then, Cytoscape 3.9.1 was used to construct the PPI network (https://cytoscape.org/). The cytoHubba plugin was used to identify hub genes.

Patients

In this study, from November 2021 to June 2022, the primary diagnosed MM patients with complete data were enrolled, and the patients with greater than 20% MM cells of bone marrow. The institutional ethics committee approved the study with informed consent (The ethics serial number is QYFY WZLL 28115.).

Real-time PCR validation

Real-time PCR was performed to validate the expression of HIST1H2BH in both MM patients and healthy controls. Primer Premier 5.0 software was used to design the HIST1H2BH primer. The primer sequences were as following: Fwd primer: GATCCAGCTAAGTCCGCTC, Rev primer: ATGGCTTTGGAGGAGATGCC.

The PCR reactions were performed in a total volume of 20 µL containing 10 µL of 2X ChamQ Universal SYBR qPCR Master Mix. The amplification program was as follows: initial hot start (30s, 95 °C) and 40 cycles with denaturation (10s, 95 °C), annealing (30s, 60 °C), extension (60s, 72 °C) with a final melting curve (15s, 95 °C) (60s, 60 °C) (15s, 95 °C).

Gene set enrichment analysis and correlation analysis

Gene set enrichment analysis (GSEA) was used as a computational method to determine if a pre-defined set of genes shows statistically significant differences between two biological states. The program multiple GSEA was used to analyze the gene expression data from TCGA. The selected criteria were P < 0.05, and only pathways related to tumors were retained. The differentially expressed genes, up- and down-regulated, were then subjected to a correlation analysis.

Survival curve, independent prognosis analysis and clinical correlation analysis

Differential gene expression analysis, pathway enrichment analysis (GO and KEGG) and protein network analysis were performed on the data from GEO. A critical hub gene analysis was conducted and HIST1H2BH was selected as the target gene for further study. Independent prognostic analysis and clinical correlation analysis were conducted based on the data from the TCGA (https://portal.gdc.cancer.gov).

Cell viability

The proliferation activity of HIST1H2BH was tested on U266 and RPMI8226 cell lines. Stable cell lines were established using the lentivirus system, and puromycin was used to establish the HIST1H2BH knockdown (siHIST1H2BH) stable lines (1.5μg/ml for RPMI-8266 and 2.5μg/ml for U266). These cell lines were seeded in 96-well plates (1 × 10⁴ cells/well), and cell viability was evaluated with CCK-8 assay. The assay was performed according to the manufacturer’s protocol.

Results

Differential genes analysis of MM dataset

Data downloaded from the GEO database was annotated to convert the gene probe matrix into a gene matrix. DEG analysis was then performed, and a total of 793 DEGs were obtained, including 285 up-regulated genes and 508 down-regulated genes (Figure 1A). The volcano plot was conducted by the ggplot2 package (Figure 1B). The X-axis represents the log2-transformed FC and the Y-axis represents the adjusted P-value. The abscissa-above red dots represent up-regulated genes, while the abscissa-below red dots represent down-regulated genes in MM patients. The black dots represent genes with no significant differences. Correlation analysis was performed on the top 40 differentially expressed genes regarding up- and down-regulation. The green circles indicate a negative correlation and red circles indicate positive correlation. The circle size has a direct ratio to the correlation results (Figure 1C).

Figure 1. The expression of DEGs based on GEO database datasets. (A) Heatmap showed 793 DEGs, including 285 upregulated genes and 508 down-regulated genes and (B) Violin diagrams, the abscissa-above red dots represent upregulated genes in MM patients, the abscissa-below red dots represent down-regulated genes in MM patients and the black dots represent genes with no significant differences. (C) Correlation plot of 40 DEGs. The positive correlation was marked with red and negative correlation was marked with blue. The size of circle represents the absolute value of correlation coefficients.

GO and KEGG analysis

GO analysis was performed on the DEGs. The results of the biological process analysis showed that the differentially expressed genes associated with tumor were mainly enriched in subsets such as leukocyte-mediated immunity, leukocyte migration and activation of the immune response etc. Cellular component analysis showed that the differentially expressed genes were mainly enriched in subsets such as collagen-containing extracellular matrix and the external side of the plasma membrane etc. Molecular function analysis showed that the differentially expressed genes were mainly enriched in subsets such as serine-type peptidase activity and serine hydrolase activity (Figure 2A) etc. KEGG pathway enrichment analysis of the DEGs showed enrichment in pathways such as the PI3K-Akt signaling pathway, hematopoietic cell lineage and cell adhesion molecules (Figure 2B) etc. GSEA analysis of the entire gene sets in the dataset showed enrichment of pathways mainly related to cell cycle, DNA replication, homologous recombination, P53 signal pathway, RNA degradation and spliceosome (Figure 2C).

Figure 2. The GO, KEGG enrichment and GSEA analyses using the 793 DEGs. (A) GO functional analysis shows enriched items in the DEGs. (B) KEGG analysis shows the enriched items in the DEGs. BP, biological process; CC, cellular component; MF, molecular function. (C) GSEA plot of genes involved in cell cycle, DNA replication, homologous recombination, P53 signaling pathway, RNA degradation, and spliceosome signaling pathway.

The DEGs obtained from the screening were analyzed on the String website to generate a protein interaction network. The interaction score with the highest confidence (0.9) was chosen as the minimum required interaction score (Figure 3A). Cytoscape was used to select hub genes, and the cytoHubba plugin was used to screen for functional modules in the PPI network to identify sub-networks. The sub-net scored four and showed interactions between HIST1H2AC, HIST1H2BH, CCND1 and TCF7L2 (Figure 3B).

Figure 3. The protein–protein interaction (PPI) network and 4 hub gene interaction network analyzed using Cytoscape software (A) PPI network. The interaction score with the highest confidence (0.9) was chosen as the minimum required interaction score. (B) The sub-net scored four and showed interactions between HIST1H2AC, HIST1H2BH, CCND1 and TCF7L2.

Survival curves, independent prognosis analysis and clinical correlation analysis

The survival package was used to perform multivariate Cox regression analysis according to patient HIST1H2BH expression levels, age, disease stage and gender in the TCGA datasets (Figure 4A). The results showed that HIST1H2BH levels, age and stage could be used as independent prognostic factors for MM, while gender was not an independent factor associated with patient survival. The expression of HIST1H2BH showed an increasing trend with the disease progression, although it was not statistically significant (Figure 4B). The expression of the HIST1H2BH (Figure 5A, P < 0.05), CCND1 (Figure 5B, P < 0.05), and HIST1H2AC (Figure 5C, P < 0.001) genes in MM patients was up-regulated compared to normal individuals. While the expression of the TCF7L2 gene in MM patients was down-regulated compared with normal individuals (Figure 5D, P < 0.001). The HIST1H2BH gene of MM patients in the high-risk group was more likely to have a reduced survival time compared with those in the low-risk group (Figure 5E, P < 0.001), while TCF7L2 gene was direct opposite (Figure 5H, P < 0.05). There were no significant difference on HIST1H2AC and CCND1 (Figure 5F and G, P > 0.05).

Figure 4. (A) Multivariate Cox regression analysis according to HIST1H2BH expression levels, patient age, gender and disease stage based on DEGs in the TCGA dataset. (B) HIST1H2BH expression levels on different clinical stage based on TCGA dataset.

Figure 5. The expression of HIST1H2BH increased in MM samples and siHIST1H2BH can decrease the proliforation of MM cell lines. Relative HIST1H2BH expression was determined in MM samples and healthy people by quantitative real-time PCR analysis. (B and C) Cell growth of siHIST1H2BH and siControl transfected U266 and RPMI8226 cell lines was analyzed by CCK-8 assays.

The expression of HIST1H2BH in MM patients and HIST1H2BH effects on MM cells proliferation

Sample clinical features have been shown in Table 1. Real-time PCR was used to detect the expression of HIST1H2BH in MM primary patients and the healthy controls. And the correlation between HIST1H2BH expression and clinical characteristics of MM patients was shown in Table 2, it was found that there was a statistically significant difference in bone lesion and renal impairment groups. Additionally, RT–PCR result showed that the significant increased expression of HIST1H2BH in MM patients compared to that of the healthy controls (Figure 5A, P < 0.05).

Table 1. Sample clinical features of MM.

Characteristics	MM	HD	P-value
Age (years)
≤60	11	5	0.737
>60	20	7
Sex
Male	18	8	0.735
Female	13	4
R-ISS stage
I/II	14	NA
III	15	NA
Light chain type
Κ	17	NA
Λ	12	NA
Heavy chain type
IgG	18	NA
IgA	6	NA
IgD	1	NA
Bone lesion
yes	27	NA
no	4	NA
Renal impairment
yes	8	NA
no	23	NA
LDH, U/L
≥250	8	NA
<250	23	NA
ALB, g/L
≥35	11	NA
<35	20	NA
β2-MG
<5.5	16	NA
≥5.5	15	NA

Note: MM: multiple myeloma; HD: healthy donors; NA: not applicable.

Table 2. Correlation between HIST1H2BH expression and clinical characteristics of MM patients.

Characteristics	MM	HIST1H2BH expression	P-value
Age (years)			0.508
≤60	11	12.38 ± 15.72
>60	20	14.11 ± 14.96
Sex			0.945
Male	18	13.06 ± 13.93
Female	13	14.09 ± 16.91
R-ISS stage			0.383
I/II	14	12.12 ± 14.69
III	15	15.53 ± 16.14
Light chain type			0.578
Κ	17	13.39 ± 10
Λ	12	15.7 ± 14.87
Heavy chain type			0.219
IgG	18	11.9 ± 9.92
IgA	6	7.07 ± 6.84
IgD	1	15.73
Bone lesion			0.038
yes	27	14.86 ± 15.69
no	4	4.23 ± 4.64
Renal impairment			<0.0001
yes	8	32.67 ± 16.95
no	23	6.83 ± 6.44
LDH, U/L			0.586
≥250	8	13.81 ± 10.42
<250	23	13.14 ± 16.6
ALB, g/L			0.115
≥35	11	17.56 ± 16.44
<35	20	11.26 ± 14.08
β2-MG			0.087
<5.5	16	8.02 ± 8.16
≥5.5	15	19.33 ± 18.55

Proliferation activity of MM cell lines was measured by CCK-8 assay to assess short-term cell viability of siHIST1H2BH transfected MM cells. Absorbance values were measured every 24 h. After 96 h of continuous incubation, the proliferation curve showed that the HIST1H2BH-knockdown cells (defined as the siHIST1H2BH group) exhibited significantly slower proliferation when compared with the control group at 24, 48, 72 and 96 h for both U266 and RPMI8226 cell lines (p < 0.05, Figure 5B and C). Our results suggested that down-regulation of HIST1H2BH decreased MM cell proliferation.

Discussion

Until now, MM still remains incurable because novel treatments and plenty of biomarkers that reflect host- or tumor-related factors have failed to explain the heterogeneity in clinical outcomes [27–31]. Recently, evaluation systems based on prognostic parameters have been established for stratifying risk [32–35]. Therefore, it is crucial to identify novel clinical marker to predict prognosis and tailor treatment plans for MM patients.

An epigenetic trait is information beyond DNA sequence that can be inherited from one generation to the next. Extensive research over the years has revealed that histone modifications, DNA methylation and RNA-mediated mechanisms are the primary methods of regulating epigenetics. During transcription, histone modifications like ubiquitination, acetylation and methylation happen simultaneously and are crucial for successful transcription [20]. Evidence suggests that epigenetic regulation of MM is closely linked to histone modification-related proteins [18, 19]. Post-translational modification of histones refers to the addition of different acylation modifications to the amino acid residues at the tail of histones, which can regulate gene transcription and expression by altering chromatin structure without altering DNA sequence [36, 37].

In histone modification, the H2B protein plays an important role [20–22], Tathe et al. found that SHP-1-mediated dephosphorylation of Y121 may be a prerequisite for efficient ubiquitination of H2B in the absence of Y121 phosphorylation. SHP-1-mediated H2B dephosphorylation contributes to maintaining basal autophagic flux in cells through the efficient transcription of genes associated with autophagy and lysosomes [20]. H2B is an essential post-translational modification, Wang’s study provides new mechanistic insights into chromatin remodeling in reprogrammable cells by showing that RNF20-mediated H2B ubiquitination promotes chromatin relaxation [21]. In addition, H2B involved in differentiation [38, 39], transcription initiation and elongation [40, 41], meiotic recombination [42] and other fundamental physiological processes. Monoubiquitination occurs at Lysine 120 (K120) at which histone H2B (H2Bub1) is monoubiquitinated, causing chromatin strands to be physically disrupted, leading to open and accessible fibers [43]. Studies have demonstrated that mutations in H2B may play an important role in the development of cancer. As a key member of the H2B histone family, HIST1H2BH also serves as a predictor of cervical cancer outcomes [23]. At present, only a few publications focus on the role of HIST1H2BH in tunors, such as head and neck squamous cell carcinoma [44] and endometrial cancers [45]. However, there is no literature about the effect of HIST1H2BH on MM (Figure 6).

Figure 6. The expression of four different genes in tumor and normal tissue from TCGA cohort (A.HIST1H2BH gene, B.CCND1 gene, C.HIST1H2AC gene, D.TCF7L2 gene) and Kaplan-Meier curves for the OS of patients in the high-risk group and low-risk group in the TCGA cohort for four genes (E.HIST1H2BH gene, F.CCND1 gene, G.HIST1H2AC gene, H.TCF7L2 gene).

In the present work, we compared HIST1H2BH gene expression datasets from the GEO database with MM samples and identified 793 DEGs. Among them, the expression level of 285 was up-regulated, while the expression level of 508 was down-regulated. According to GO analysis, the results of the biological process analysis showed that the differentially expressed genes associated with tumor were mainly enriched in subsets such as leukocyte-mediated immunity, leukocyte migration and activation of the immune response etc. Cellular component analysis showed that the differentially expressed genes were mainly enriched in subsets such as collagen-containing extracellular matrix and the external side of the plasma membrane etc. Molecular function analysis showed that the differentially expressed genes were mainly enriched in subsets such as serine-type peptidase activity and serine hydrolase activity etc. KEGG pathway enrichment analysis of the DEGs showed enrichment in pathways such as the PI3K-Akt signaling pathway, hematopoietic cell lineage and cell adhesion molecules etc. Additionally, we built a PPI network containing 127 nodes and 413 edges. This PPI network's top module was then screened with the cytoHubba plugin. In the meantime, 4 DEGs were discovered from the as-built PPI network. Furthermore, survival analysis revealed that the selected hub gene HIST1H2BH was strongly associated with poor OS in MM patients. Finally, we confirmed that HIST1H2BH was up-regulated in MM primary samples via qPCR. And CCK-8 assay showed knockdown HIST1H2BH could decrease the proliferation of MM cells. In conclusion, the present study had identified HIST1H2BH as the hub gene associated with multiple myeloma using integrated bioinformatics analysis and validated by in vitro experiments. According to the present research, HIST1H2BH level was expressed significant higher in the MM patients group compared to the healthy donors, and the survival curve results showed a negative correlation between HIST1H2BH expression and MM prognosis. Our results suggested that HIST1H2BH may functioned as an oncogene in MM, that could be therapeutic target for personalized and precision treatment of MM patients. As a key member of the H2B histone family, HIST1H2BH may function as a regulator through epigenetic regulation. Its’ biological function requires further exploration in the future.

Authors’ contributions

Conceived and designed the experiments: W.X.Z., X.C. and Y.Y.Z. Performed the experiments: W.X.Z., X.C., Y.Y.Z., Q.W., J.F., Z.W., D.X.W., Z.C.W. Analyzed the data: J.X.T., S.Y.Q. Contributed reagents/materials/analysis tools: W.X.Z., X.C. and Y.Y.Z. Wrote the paper: W.X.Z., X.C. and Y.Y.Z. The authors read and approved the final manuscript.

Limitation of the present study

The major limitation of the present study is the small sample size. The results would be more reliable with more enrolled MM patients.

Ethics approval and consent to participate

The present study was approved by the Ethics Committee of the Affiliated Hospital of Qingdao University. The ethics serial number is QYFY WZLL 28115.

Availability of data and materials

The datasets used and/or analyzed during the current study are available from the corresponding author on reasonable request.

Disclosure statement

No potential conflict of interest was reported by the author(s).

Additional information

Funding

The present study was funded by the Natural Science Foundation of Shandong Province, ZR2020MH314 (to Xian Chen), the Natural Science Foundation of Shandong Province, ZR2021MH311 (to Yunyuan Zhang).