Identification of potential pathways and microRNA-mRNA networks associated with benzene metabolite hydroquinone-induced hematotoxicity in human leukemia K562 cells

Hydroquinone (HQ) is a phenolic metabolite of benzene with a potential risk for hematological disorders and hematotoxicity in humans. In the present study, an integrative analysis of microRNA (miRNA) and mRNA expressions was performed to identify potential pathways and miRNA-mRNA network associated with benzene metabolite hydroquinone-induced hematotoxicity. K562 cells were treated with 40 μM HQ for 72 h, mRNA and miRNA expression changes were examined using transcriptomic profiles and miRNA microarray, and then bioinformatics analysis was performed. Out of all the differentially expressed genes (DEGs) and differentially expressed miRNAs (DEMs) induced by HQ, 1482 DEGs and 10 DEMs were up-regulated, and 1594 DEGs and 42 DEMs were down-regulated. HQ-induced DEGs were involved in oxidative stress, apoptosis, DNA methylation, histone acetylation and cellular response to leukemia inhibitory factor GO terms, as well as metabolic, Wnt/β-catenin, NF-κB, and leukemia-related pathways. The regulatory network of mRNAs and miRNAs includes 23 miRNAs, 1108 target genes, and 2304 potential miRNAs-mRNAs pairs. MiR-1246 and miR-224 had the potential to be major regulators in HQ-exposed K562 cells based on the miRNAs-mRNAs network. This study reinforces the use of in vitro model of HQ exposure and bioinformatic approaches to advance our knowledge on molecular mechanisms of benzene hematotoxicity at the RNA level.


Open Access
*Correspondence: yizc@buaa.edu.cn 1 School of Biological Science and Medical Engineering, Beihang University, Beijing 100191, China Full list of author information is available at the end of the article [14,15]. The underlying molecular mechanisms involved in HQ toxicity are not fully understood.
MicroRNAs (miRNAs) are endogenous, highly conserved small non-coding RNAs that modulate various biological processes like cell survival, proliferation, metabolism, differentiation, and apoptosis [16][17][18][19]. The miRNA-mediated gene silencing can involve translational repression, co-translational protein degradation, competition for the cap structure, inhibition of ribosomal subunit joining, inhibition of mRNA circularization through deadenylation, or deadenylation-dependent decay [20,21]. But recently, miRNAs have been found to mediate gene activation through the bidirectional transcription of the human genome, binding to the enhancer, or recruiting a protein complex with transcriptional activators to the gene promoter [22,23]. Repression is more common in eukaryotes whereas posttranscriptional upregulation has been observed in specific cell types with distinct transcripts or conditions [23]. MiRNAs have been identified as exerting a powerful effect upon acute myeloid leukemia (AML) development and can be used as potential biomarkers for leukemia diagnosis and prognosis [24][25][26]. Various studies have assessed the role of miRNAs in particulate matter-exposed human pulmonary epithelial cells, sulfur mustard-resistance of the keratinocyte cell line, alternariol and altertoxin II-treated HepG2 cells, and environmental toxicants-exposed aquatic organisms to explore the corresponding toxicity mechanisms [27][28][29][30]. As Liang et al. point out, miRNA-451a and miRNA-486-5p expression are notably lower in HQ-treated CD34+ hematopoietic progenitor cells and K562 cells [13]. HQ may activate apoptotic signals via inhibiting the tumor-suppressive effects of miR-7-5p in TK6 lymphoblastoid cells [31]. However, the role of miRNAs and their target mRNAs expression in benzene hematotoxicity has not been fully addressed yet and needs further study.
Transcriptome sequencing with high sensitivity and good reproducibility is a bargain for detecting lowexpression genes [32]. The integration of genomic tools contributes to investigating the molecular mechanism of toxicity [33]. Human leukemia K562 cells were derived from a patient with chronic myeloid leukemia [34]. Here we performed transcriptomic profiles and miRNA microarray to identify mRNAs and miRNAs changes, constructed the mRNAs and miRNAs regulatory network, and performed an in-depth bioinformatic analysis in HQ-induced K562 cells. The differentially expressed mRNAs and miRNAs participated in the metabolic pathways, Wnt/β-catenin pathway, NF-κB pathway, etc. by regulating oxidative stress, apoptosis, DNA methylation, histone acetylation, resulting in hematotoxicity including leukemia. These findings provided a theoretical basis for understanding the molecular mechanisms of benzene hematotoxicity, laying the foundation for future validation of in vivo models as well as therapeutic targets and prognostic factors.

Cell culture
Human leukemia K562 cells were purchased from Cell Resource Center, PekingUnion Medical College (CRC/ PUMC, China), and were cultured as described previously [15,34]. After K562 cells were treated with 40 μM HQ (Sigma-Aldrich) for 72 h, the cells were harvested for further study.

Transcriptome analysis
The genome-wide transcriptome analysis was analyzed as described previously [35]. Differentially expressed genes (DEGs) analysis was indicated with the absolute log 2 (fold change of HQ/C) values and the adjusted P-value <0.05 was considered differentially expressed. Data are representative of three independent experiments. The volcano plot of DEGs was performed by creating scatter plots in Excel software. Select the data of downregulated DEGs including log 2 (fold change of HQ/C) values and the adjusted P-value; choose the scatter plot to represent the relationship between the data sets; Add green color for the dots. The data of upregulated DEGs and genes without differential changes were performed in the same plot and added red and blue colors, respectively. Gene Ontology (GO) enrichment analysis of DEGs was conducted by over-representation analysis (ORA) using the online database WebGestalt (http:// www. webge stalt. org). Significantly enriched GO terms in DEGs compared to the genome background were defined by Wallenius' non-central hypergeometric distribution adjusting for gene length bias [36]. The pathway analysis of DEGs was conducted using the Kyoto Encyclopedia of Genes and Genomes (KEGG) database (http:// www. kegg. jp/) [37].

miRNA microarray analysis
The miRNA microarray was analyzed as described previously [38]. Total RNA was extracted with Ribozol ™ RNA Extraction Reagent (Invitrogen, USA) according to the manufacturer's protocol. The expression of microRNAs contained in the miRBase was analyzed by microarray using the μParaflo ™ Microfluidic Biochip Technology of LC sciences (Huston, TX, USA). A total of 1 μg of miRNA-enriched RNA was labeled with Cy3 using the ULS ™ microRNA labeling kit (Kreatech, USA) and hybridized on the microarray. Hybridization images were collected using a laser scanner and digitized using Array-Pro image analysis software. The signals were normalized using a locally-weighted regression (LOWESS) filter.
Differentially expressed miRNAs (DEMs) analysis was indicated with the absolute log 2 (fold change of HQ/C) values and P-value <0.05 was considered differentially expressed. Data are representative of three independent experiments. Hierarchical clustering analysis of the DEMs was performed and graphs were generated using the ggplot2 package in R software.

The effects of HQ on DEGs and DEMs in K562 cells
The volcano plot of overall gene expression showed that 3076 genes were differentially expressed in K562 cells treated with 40 μM HQ for 72 h, including 1482 upregulated and 1594 downregulated DEGs (Fig. 1a). Moreover, there were 84 upregulated and 307 downregulated DEGs over a 2-fold change after HQ exposure. The top 10 upregulated DEGs and top 10 downregulated DEGs based on the absolute log 2 (fold change of HQ/C) values after HQ exposure for 72 h were set out in Table 1. What stands out in the table is that LILRA6 (log 2 (HQ/C) = 4.17) was the most upregulated DEG and HBB (log 2 (HQ/C) = −4.95) was the most downregulated DEG.
The heatmap and hierarchical clustering analysis showed that 52 miRNAs were differentially expressed after HQ exposure, including 10 upregulated and 42 downregulated DEMs (Fig. 1b). In addition, 7 upregulated and 23 downregulated DEMs were over 2-fold change after HQ exposure. These results indicated that majority of genes and miRNAs were downregulated after HQ exposure. The top 10 upregulated DEMs and top 10

DEGs involved in the oxidative stress, apoptosis, DNA methylation, histone acetylation and cellular response to leukemia inhibitory factor GO terms in HQ-induced K562 cells
The top 30 enriched GO terms of target DEGs were presented in Fig. 2. GO analysis showed that HQ-upregulated DEGs enriched in biological processes for RNA processing, cellular nitrogen compound metabolic process, and ribonucleoprotein complex biogenesis, whereas HQ-downregulated DEGs enriched in biological processes for single-organism cellular process, localization, and response to stimulus. Our previous studies have demonstrated that DNA methylation, histone acetylation, and ROS have major influences on the transcription of erythroid-specific genes [15,[39][40][41]. Here, we focused on the GO terms of oxidative stress, apoptosis, DNA methylation, and histone acetylation in HQ-induced K562 cells. HQ treatment considerably increases the relative ROS levels in HD3 chicken erythroblast cells, HL-60 promyelocytic leukemia cells (HL-60 cells), Jurkat T-lymphoblastic leukemia cells (Jurkat cells), and K562 cells [4,40,42,43]. In the present study, there were 38 upregulated DEGs and 37 downregulated DEGs in the response to oxidative stress term (GO:0006979), among which there were 27 upregulated DEGs and 20 downregulated DEGs in the cellular response to oxidative stress term (GO:0034599) in HQ-induced K562 cells. Taken together, these indicated the crucial role of oxidative stress in HQ exposure in K562 cells.
Activation of CASP3 and CASP8 both have pivotal roles in the execution phase of cell apoptosis. As mentioned, HQ induces apoptosis in HL-60 cells, Jurkat cells, human bone marrow mononuclear cells (BMMNC), and K562 cells accompanied by CASP3 or CASP8 activation [7,44,45]. Consistently, HQ exposure markedly upregulated the mRNA levels of CASP3 and CASP8 to 1.53-fold and 1.38-fold of that in the control K562 cells in the present study. Moreover, there were 17 upregulated DEGs and 14 downregulated DEGs in the positive regulation of apoptotic signaling pathway (GO:2001235), contributing to HQ-induced apoptosis in K562 cells.
Exposure to HQ can markedly change the mRNA levels and DNA methylation levels of erythroid-specific genes, as well as reactive oxygen species (ROS) levels in K562 cells [14,[39][40][41]. Transcriptomic analysis contributed to high-throughput data for HQ-induced DEGs in K562 cells, which further support the role of oxidative stress, apoptosis, DNA methylation, histone acetylation in benzene metabolite HQ-induced hematotoxicity. In addition, there were 20 upregulated DEGs and 10 downregulated DEGs in the cellular response to leukemia inhibitory factor term (GO:1990830) in HQ-induced K562 cells (Table 3). Different biological processes in response to HQ exposure might result in these genes expression changes and take part in the mechanism of benzeneinduced leukemia.

DEGs involved in the metabolic, Wnt/β-catenin, NF-κB, and leukemia-related pathways in HQ-induced K562 cells
The HQ-upregulated DEGs were annotated with 243 KEGG pathway whereas the downregulated DEGs Table 3 DEGs in the GO term of cellular response to leukemia inhibitory factor * represents P value < 0.05; ** represents P value < 0.01; *** represents P value < 0.005  Fig. 3. KEGG pathway analysis demonstrated that HQ-upregulated DEGs were the most significantly enriched in the ribosome biogenesis in the eukaryotes pathway (annotated 50 DEGs, P = 4.03E-24), spliceosome, RNA transport, and proteasome whereas HQ-downregulated DEGs enriched in the axon guidance pathway (annotated 28 DEGs, P = 4.96E-05), biosynthesis of amino acids, ribosome and carbon metabolism (Fig. 3). In addition, most genes were annotated as the metabolic pathways, attaining 121genes in both the upregulated and the downregulated DEGs. The cytochrome P450 family 1 subfamily A member 1 (CYP1A1) participates in metabolic pathways. The frizzled class receptor 2 (FZD2) gene encodes a protein in the beta-catenin canonical signaling pathway. HQ has been reported to considerably upregulate the mRNA level of CYP1A1, whereas downregulate the mRNA level of FZD2 in human epidermal keratinocytes cells (HEK cells) [49]. In our results, the mRNA levels of CYP1A1 and FZD2 in K562 cells showed similar trends to the HEK cells.
RNA-seq has identified that hydroquinone promotes DNA homologous recombination repair via activating the NF-κB pathway in the human osteosarcoma cell line (U2OS/DR-GFP) [50]. In our study, the mRNA level of NFKB1 was upregulated after HQ exposure, indicating the NF-κB pathway might be activated in HQ-induced K562 cells.
Furthermore, there were 10, 12, 13, and 39 DEGs annotated as the KEGG terms of chronic myeloid leukemia, acute myeloid leukemia, hematopoietic cell lineage, and human T-cell leukemia virus 1 infection in HQ-induced K562 cells, respectively (Table 4) [37]. The related DEGs were highlighted in the KEGG maps in Supplementary  Fig. 2-5. HQ exposure might change these genes' expression to activate the leukemia-related pathways, resulting in potential benzene hematotoxicity.

Identification of miR-1246 and miR-224 as major regulators in the miRNA-mRNA network in HQ-induced K562 cells
The mutual regulation of DEMs and DEGs was analyzed to construct an interaction network after HQ exposure in K562 cells. It is a widely held view that miRNAs can induce the degradation of their target mRNAs, so we focused on those miRNAs and mRNAs with opposite expression trends. There were 10 and 2304 potential miRNAs-mRNAs pairs in the network. Moreover, the miRNA-mRNA network of DEGs and DEMs over a 2-fold change was constructed by a cystoscope (Fig. 4).  The miRNA-mRNA network was constructed with 109 potential pairs containing 7 upregulated DEMs and 75 downregulated target DEGs in Fig. 4a; the network had 30 potential pairs containing 8 upregulated DEMs and 15 downregulated target DEGs in Fig. 4b. Per miRNA had 2-29 target genes and per gene was connected to 1-5 miRNAs. The lists of top miRNAs and target genes in the network were shown in Table 5 and Table 6. On the other hand, we found that majority of DEGs and DEMs were downregulated after HQ exposure.
HMDD (the Human microRNA Disease Database) is a database for experiment-based miRNA and human disease associations (http:// www. cuilab. cn/ hmdd) [51]. When the DEMs over a 2-fold change were screened from the HMDD 3.2, 5 miRNAs (miR-27a, miR-224, miR-1246, miR-18a, and miR-18b) were associated with leukemia, indicating potential roles in benzene-induced leukemia. According to the miRNA-mRNA regulatory network, there were more interactions between upregulated miRNAs and downregulated genes than the other way around. In our study, miR-1246 and miR-224 were the most upregulated DEMs with 17 target DEGs whereas miR-18a and miR-18b were downregulated in HQ-induced K562 cells. However, it remains unknown why MiR-27a was upregulated in HQ-induced K562 cells but downregulated in acute leukemia cell lines and primary samples compared to hematopoietic stem-progenitor cells [52]. Therefore, miR-1246 and miR-224 were identified to be major regulators for HQ exposure in K562 cells. These results indicated the crucial role of miR-1246 and miR-224 in benzene hematotoxicity.

Discussion
Our previous studies have demonstrated that hemininduced erythroid differentiation is concentrationdependently and time-dependently inhibited by benzene metabolites exposure (phenol, 1,2,4-benzenetriol, and hydroquinone) in K562 cells [14,15,39,46]. In the present study, the concentration of 40 μM HQ was selected to correspond to no obvious cytotoxicity but markedly inhibiting erythroid differentiation in K562 cells when exposed for 72 h [14]. Hemoglobin subunit beta (HBB) loci are associated with beta-thalassemia, sickle cell anemia, and heinz body anemias [53][54][55] . HBB was the most dramatically downregulated gene after HQ exposure, which might play important roles in benzene hematotoxicity.

WDR43
WD repeat-containing protein 43 3 up It has been reported that the NF-κB pathway is activated in the development of chronic myeloid leukemia and acute myeloid leukemia [56][57][58][59]. In our study, HQ exposure upregulated the mRNA level of NFKB1 in chronic myeloid leukemia and acute myeloid leukemia pathway, which partly supported that HQ might promote leukemia development by activating the NF-κB pathway. HQ exposure downregulated the mRNA levels of MAP 2 K2 and MAPK3 in chronic myeloid leukemia and acute myeloid leukemia pathway (Table 4). This may inhibit the MAPK signaling pathway and thus inhibit cell proliferation ( Supplementary Fig. 2, 3). These results were consistent with our previous study that HQ induced a concentration-dependent decrease in the viabilities in K562 cells [44]. The roles of metabolic pathways, Wnt/βcatenin pathway, and NF-κB pathway in benzene hematotoxicity need further study.
The downregulated miR-1246/1248 are key nodes that reveal the possible relapse mechanisms for pediatric T cell acute lymphoblastic leukemia. As reported, miR-1246 is one of the most highly enriched miRNAs in AML derived extracellular vesicles [16,60,61]. MiR-1246, a hundreds-fold alteration in microvesicles from three different leukemia cell lines (K562, Nalm-6, and Jurkat), is upregulated to activate the expression of C6orf211 and C19orf10 to promote tumor progression in patients diagnosed as chronic myeloid leukemia [62]. MiR-224 expression is considerably upregulated in the bone marrow of pediatric AML patients and can be used as noninvasive biomarkers for the early prediction of hepatocellular carcinoma development [63]. These results indicated the crucial role of miR-1246 and miR-224 in hematotoxicity. MiR-29a is highly expressed in arsenicinduced peripheral neuropathy, which is consistent with higher expression after HQ exposure in this study [64]. Low expression of miR-18a distinguishes pediatric and adult acute lymphoblastic leukemia from each other [65]. PML/RARα-regulated miR-181a/b cluster targets the tumor suppressor RASSF1A in acute promyelocytic leukemia [66].
In further research, the expression, target genes, and biological function of miR-1246 and miR-224 will be confirmed by more techniques in HQ-induced K562 cells, CD34+ hematopoietic progenitor cells, U937 human leukemia cells, and human peripheral blood mononuclear cells (PBMCs). Inhibition of miR-1246 and miR-224 will contribute to their regulation for HQ exposure. On the other hand, hydroquinone has been reported to inhibit PRV infection in mouse neuroblastoma N2a cells, protect neurons from transient cerebral ischemia, and reduce gliosis in a gerbil model of transient cerebral ischemia [67,68]. The study on the effect of HQ on neural systems will contribute to a full understanding of benzene toxicity. Further assessment in suitable animal models or any data of miRNAs expression of benzene-exposed patients will be immensely beneficial to further studies.

Conclusion
In summary, the present study identified differentially expressed genes and miRNAs in HQ-induced K562 cells using transcriptomic profiles and miRNA microarray. The miRNA-mRNA network can help us better understand the molecular mechanisms between miRNAs and their target genes. MiR-1246 and miR-224 had the potential to be major regulators for HQ exposure in K562 cells based on the miRNAs-mRNAs network and were reported to be associated with leukemia, suggesting potential biomarkers for the evaluation of benzene hematotoxicity. Apart from the miRNAs-mRNAs regulation, how miRNAs regulate protein expression remains elusive. Further proteomics study is to be performed to elucidate the underlying mechanism. The GO and KEGG pathways provide a framework for further studies in suitable in vitro and animal models, which will contribute to developing new strategies for the prevention and rapid diagnosis of benzene hematotoxicity.