GSK-J4-Mediated Transcriptomic Alterations in Differentiating Embryoid Bodies.

Histone-modifying enzymes are key players in the field of cellular differentiation. Here, we used GSK-J4 to profile important target genes that are responsible for neural differentiation. Embryoid bodies were treated with retinoic acid (10 μM) to induce neural differentiation in the presence or absence of GSK-J4. To profile GSKJ4-target genes, we performed RNA sequencing for both normal and demethylase-inhibited cells. A total of 47 and 58 genes were up- and down-regulated, respectively, after GSK-J4 exposure at a log2-fold-change cut-off value of 1.2 (p-value < 0.05). Functional annotations of all of the differentially expressed genes revealed that a significant number of genes were associated with the suppression of cellular proliferation, cell cycle progression and induction of cell death. We also identified an enrichment of potent motifs in selected genes that were differentially expressed. Additionally, we listed upstream transcriptional regulators of all of the differentially expressed genes. Our data indicate that GSK-J4 affects cellular biology by inhibiting cellular proliferation through cell cycle suppression and induction of cell death. These findings will expand the current understanding of the biology of histone-modifying enzymes, thereby promoting further investigations to elucidate the underlying mechanisms.

INTRODUCTION

Chromatin remodeling by DNA methylation and histone modifications (acetylation, methylation, and phosphorylation) plays an important role in the regulation of epigenetic gene expression (Kim and Workman, 2010; Lee et al., 2010). Related genes are either activated or silenced depending on the site of methylation in lysine residues (H3K4, H3K9, H3K27, H3K36, H3K79 and H4K20). Generally, methylation at H3K4 is associated with gene activation, whereas methylation at H3K9 and H3K27 results in transcriptional silencing (Li et al., 2007; Maes et al., 2015). More complexities arise when methylation patterns are also included. For example, the ɛ-amino group of lysine residues can be mono-, di-, or trimethylated, and different patterns of methylation lead to functional diversity. Dimethylation at lysine 4 (H3K4) is associated with both active and inactive genes, whereas trimethylation is only associated with active genes. These posttranslational modifications epigenetically maintain and regulate lineage-specific gene expression during development.

Lysine demethylases (KDMs) are associated with various cellular functions and various disease states. In particular, Jumonji-C domain (JmjC)-containing histone demethylases play important roles in cellular differentiation, development and cancer (Cloos et al., 2006; Rui et al., 2010; Wissmann et al., 2007). The JMJD family has three members, JMJD3 (KDM6B), UTX (KDM6A) and UTY, which are responsible for H3K27me3 demethylation (removal of one methyl residue from H3K27me3 and H3K27me2) (Hong et al., 2007; Xiang et al., 2007). Small molecules that inhibit enzyme activity can be useful tools for observing the molecular mechanisms of different disease states and developing alternative therapies. Additionally, small molecule inhibitors provide excellent opportunities to study different aspects of cellular biology (e.g., cell cycle control, mitosis, and cell signaling) and gene expression (Weiss et al., 2007).

GSK-J4, a small-molecule inhibitor with highly efficient cell permeability (Kruidenier et al., 2012), is the first selective inhibitor of the H3K27 histone demethylases JMJD3 and UTX. Various studies have reported the use of GSK-J4 to observe the roles of JMJD3 in various biological processes (Kang et al., 2015; Lane et al., 2014; Smith and Watters, 2014). It has been recently reported that the inhibition of histone demethylase by GSK-J4 affects cellular viability and cell cycle progression of differentiating cells. Hofstetter et al. (2016) observed that the inhibition of KDM6 by GSK-J4 induced cell death in early differentiating ESCs. It has also been reported that GSKJ4 induces cell death and inhibits the growth of non-small cell lung cancer (NSCLC) cell lines (Watarai et al., 2016). Additionally, induction of cell death and loss of self-renewal capacity in cancer stem cells (CSCs) have also been reported (Sakaki et al., 2015). Despite several reported findings, the mechanisms underlying GSK-J4-mediated cell death remain elusive, particularly with respect to transcriptomics. Here, we performed a transcriptomics analysis to profile GSK-J4-target genes during early neural differentiation, providing an important list of genes that induce cell death with cell cycle suppression and inhibition of cell proliferation. To our knowledge, this is the first report to provide a defined gene list regarding GSK-J4-mediated cell death in differentiating cells.

MATERIALS AND METHODS

Culture conditions and procedures

We cultured NCCIT cells (ATCC, USA) in RPMI 1640 medium (Invitrogen, USA) supplemented with 10% fetal bovine serum (heat inactivated) (Thermo Scientific, USA), 100 U/ml penicillin and 100 μg/ml of streptomycin (Invitrogen, USA). Embryoid body (EB) formation was performed as previously described (Mandal et al., 2016). After 24 h of culturing, EBs were treated with 10 μM all-trans retinoic acid (RA; Sigma-Aldrich, USA) with or without 10 μM GSK-J4 (Tocris Bioscience, United Kingdom). The EBs were maintained in a controlled environment at 5% CO2 and 37°C during culture (48 h).

Isolating total RNAs

Total RNA samples were extracted from EB, EB+RA and EB+RA+GSK-J4 cells using RNAiso Plus (Takara BIO, Japan) according to the manufacturer’s instructions. Briefly, 200 ml of chloroform was mixed with the sample, and the sample was inverted for 5 min. After centrifugation for 15 min at 14,000 × g at 4°C, the upper solution was collected. A total of 600 μl of isopropanol was added, and the solution was subsequently incubated on ice for 1 h. The lysate was centrifuged at 14,000 × g for 15 min at 4°C, and then the isopropanol was decanted. After washing with ice-cold ethanol, the RNA pellets were dried at room temperature for 5 min, and subsequently 20 μl of diethylpyrocarbonate (DEPC)-treated water was added. The quality and quantity of extracted total RNA was measured using an Agilent 2100 Bio-analyzer (Agilent Technologies, Germany) and a spectrophotometer (NanoDrop Technologies, USA), respectively.

RNA sequencing (RNA-seq)

We used 5 μg of total RNA to remove ribosomal RNAs (rRNAs) with a RiboMinus™ Transcriptome Isolation Kit (Invitrogen). A total of 100 ng of rRNA-depleted RNA was used to construct paired-end transcriptome libraries with the NEBNext® Ultra™ Directional RNA Library Prep Kit for Illumina® (New England Biolabs, USA). Briefly, first-strand cDNA was synthesized, followed by second-strand synthesis with DNA polymerase I and RNase H. The double-stranded cDNAs were subsequently end-repaired and ligated to adaptors. Fragments of 300–400 bp in size were purified using the MinElute Gel Extraction Kit (Qiagen, Germany). For further enrichment, the fragments were amplified and purified using ethanol precipitation, and cDNA fragments of 101 bp in length were sequenced using Illumina HiSeq2500 (National Instrumentation Center for Environmental Management in Seoul National University) with two biological replicates. The gene expression data were submitted to the NCBI Sequence Read Archive (SRA) repository (http://www.ncbi.nlm.nih.gov/sra/) under accession numbers SRX1175001, SRX690524 and SRX690525.

Analysis of the sequencing data

The FASTQ files obtained from the RNA-seq experiments were processed using Trimmomatic software (Bolger et al., 2014). The trimmed FASTQ files were subsequently aligned to the UCSC hg19 reference genome using STAR (version 2.5.1) aligner software (Dobin et al., 2013) with three mismatches. DESeq2 (Love et al., 2014) was used to measure differential gene expression, applying a 1.2 log2-fold difference and P < 0.05. We selected a moderate threshold of a log2-fold change of 1.2 to avoid a low number of DEGs.

Functional annotation of differentially expressed genes (DEGs)

The functional annotation of DEGs was performed using the web accessible program Database for Annotation, Visualization and Integrated Discovery (DAVID, version 6.8) (Huang da et al., 2009). A modified Fisher’s Exact p-value (< 0.05) was used to demonstrate gene ontology (GO). A heat map was constructed to view the relative expression patterns of all DEGs using MultiExperiment Viewer (MeV), enabling advanced analysis through instinctive graphics. Ingenuity pathway analysis (IPA) (Kramer et al., 2014) was also performed to detect cellular and molecular functions and enrich specific pathways.

Obtaining gene expression graphs from the UCSC genome browser

After receiving normalized data from HOMER (Heinz et al., 2010), we used the UCSC genome browser to draw gene expression graphs. The related files were uploaded to the UCSC browser by applying custom track addition. After successful uploading, we used a track setting option to manage the vertical viewing range at a uniform scale to compare the expression of a particular gene. Finally, the expression graphs were exported and edited.

Cell proliferation assay

The proliferation of treated cells was evaluated using a Premix WST-1 Cell Proliferation Assay System (Takara Bio, Inc., Japan). After 3 days of culture, the culture medium was removed, and the cells were washed with phosphate-buffered saline (PBS). WST-1 reagent was then added, and the cells were incubated for 4 h. The results of the WST-1 assay were measured at 450 nm using a Model 680 microplate reader (Bio-Rad, USA).

Cell cycle analysis by flow cytometry

After treating RA-induced EBs with GSK-J4 for 48 h, the cells were harvested to analyze the cell cycle distribution using 0.25% trypsin-EDTA (Invitrogen Life Technologies). After washing with PBS, the cells were probed using a BD CycleTest™ Plus DNA Reagent kit (BD Biosciences, USA) according to the manufacturer’s instructions. The distribution of cycles were analyzed using a FACSCalibur (BD Biosciences, USA). The results are presented as the percentage of cells per phase, which was calculated using ModFit LT 3.0 (Verity Software House, USA).

Apoptosis assays

Cellular apoptosis was measured by FACS using an Annexin V-FITC apoptosis kit from Clontech Laboratories, Inc. (USA). The apoptotic cells in both the early and late stages were identified based on the localization of Annexin V and propidium iodide (PI). The collected EBs were washed twice with ice-cold PBS and subsequently incubated with 200 μl of 1X binding buffer containing 5 μl of Annexin V-FITC and 10 μl of PI. After incubation for 15 min, the cells were analyzed for apoptosis using a flow cytometer.

Western blot analysis

Whole cell extracts were prepared with RIPA buffer supplemented with complete, EDTA-free Protease Inhibitor Cocktail (Roche, Germany). Total proteins were separated via sodium dodecyl sulfate polyacrylamide gel electrophoresis (SDS-PAGE) and transferred to polyvinylidene difluoride (PVDF) membranes (Schleicher & Schuell Bioscience, Inc., USA). After blocking with PBS supplemented with 5% horse serum (Invitrogen) for 2 h, the membranes were blotted with appropriate antibodies at 4°C for 16 h, followed by blotting with HRP-conjugated anti-rabbit IgG antibody for 1 h at room temperature. After washing with TBST (20 mM Tris-Cl, pH 7.6; 500 mM NaCl; 0.1% Tween 20), the reactions were detected using ECL Prime Western Blotting Detection Reagent (GE Healthcare, USA). The antibodies used for western blotting were purchased from Abcam (United Kingdom).

Quantitative real time PCR (qRT-PCR) analysis

We used an ABI 7500 Real-Time PCR System (Applied Bio-systems, Inc., USA) to assay relative mRNA expression levels with SYBR Premix Ex TaqTM II (Otsu-Shi, Japan). The total reaction volume of each sample was 20 μl, and the PCR conditions were as follows: 30 s at 95°C, 40 cycles of 5 s. at 95°C and 34 s. at 60°C, followed by a melting curve analysis step. PCR reactions were considered specific when all amplicons showed a single Tm. Relative quantifications were measured in triplicate and evaluated using the comparative CT (ΔΔCT) method. Expression data were normalized relative to the expression of the GAPDH gene. The primers used for gene validation are listed in Supplementary Table 1.

Statistical analysis

We performed RNA-sequencing using biologically duplicates for both control and GSK-J4-treated samples. To draw a heat map and analyze the data we used a normalized version of dual RNA-seq. For qRT-PCR analysis, we used triplicate biological, and the results are presented as the mean ± standard error of the mean (SEM). For statistical analyses, Student’s t-test was performed using Microsoft Office Excel, 2010 (p < 0.05).

RESULTS

Genes were differentially expressed in GSK-J4-treated EBs

The NCCIT cell-derived EBs were treated with GSK-J4 to analyze the effects during the early stage of neural differentiation. The EBs were exposed to 10 μM of GSK-J4 for 48 h. The detailed schedule of GSK-J4 treatment is shown in Fig. 1A. To observe the effects of GSK-J4 at the gene expression level, we performed RNA sequencing for EB+RA and EB+RA+GSK-J4. Many genes were differentially expressed after GSK-J4 treatment compared with the control. We identified 47 and 58 genes that were up- and down-regulated in EB+RA vs. EB+RA+GSK-J4, respectively, with a cut-off value of 1.2 log2-fold changes (p-value < 0.05). Based on the fold-change criteria used to mine the biological data, a heat map was constructed showing that the expression of genes in GSK-J4-treated EBs was clearly different from that in the controls (Figs. 1B and 1C).

Fig. 1

Functional annotation of differential gene expression

(A) Graphical experimental scheme for the differentiation/treatment protocol. NCCIT cells were stabilized and sub-cultured to form EBs. After stabilization, the EBs were treated with or without GSK-J4 for 48 h. The samples were subsequently collected for further analysis. (B) and (C) show heat maps of up- and down-regulated genes expressed between EB+RA (ER) and EB+RA+GSK-J4 (ERG), respectively. The gene expression level of each gene in the heat map is scaled and represented as relative expression changes. The color key shows relative expression within a range of 0 to 100. (D) and (E) show the bar graphs from the gene ontology analysis (Biological Process) of up- and down-regulated genes, respectively. The top 10 significant categories are plotted here. Related p-values are indicated in front of each bar. The total number of up- and downregulated genes is indicated in parentheses in the title of the bar graphs.

Functional annotation of all DEGs revealed significant enrichments in important biological categories

We used our significant DEGs for gene ontology analysis using the DAVID analysis tool. We identified 9 categories (biological processes) with up-regulated genes (Threshold count: 6, EASE score: 0.05 with Benjamini correction). In contrast, 4 categories were enriched with down-regulated genes (Threshold count: 6, EASE score: 0.05 with Benjamini correction). We observed that the genes up-regulated in GSK-J4-treated cells were largely involved in the regulation of cell proliferation, programmed cell death, and apoptosis, among others. The down-regulated genes were involved in the regulation of transcription and, the RNA metabolic process, among others. The overall GO results showed that GSK-J4 altered several important genes related to basic cellular activities and cellular differentiation. The full chart of enriched categories of both up- and down-regulated genes is presented in Figs. 1D and 1E.

We also conducted pathway analysis using Ingenuity Pathway Analysis (IPA) and observed that the DEGs were significantly involved in a total of 17 pathways (FDR < 0.05, a minimum of 3 genes represented). Important significantly enriched signaling pathways included GADD45, p53, and Wnt/β-catenin signaling, among others. A detailed list of the related genes is presented in Table 1.

Table 1

List of canonical pathways enriched with differentially expressed genes

Ingenuity canonical pathways	−Log (p-value)	Molecules
GADD45 Signaling	4.1	GADD45B, GADD45G, CCND1
p53 Signaling	2.81	JUN, GADD45B, GADD45G, CCND1
ErbB2-ErbB3 Signaling	2.44	JUN, ETV4, CCND1
ATM Signaling	2.26	JUN, GADD45B, GADD45G
Wnt/β-catenin Signaling	2.16	JUN, FZD4, CDH5, CCND1
Glucocorticoid Receptor Signaling	2.04	HMGB1, JUN, CDKN1C, ANXA1, HSPA5
Adipogenesis pathway	1.66	FZD4, TXNIP, FGFRL1
Colorectal Cancer Metastasis Signaling	1.61	JUN, ARRB1, FZD4, CCND1
Hereditary Breast Cancer Signaling	1.59	GADD45B, GADD45G, CCND1
Human Embryonic Stem Cell Pluripotency	1.58	FZD4, FGFRL1, LEFTY2
Ovarian Cancer Signaling	1.58	ARRB1, FZD4, CCND1
Acute Phase Response Signaling	1.4	HMOX1, FN1, JUN
Role of Macrophages, Fibroblasts and Endothelial Cells in Rheumatoid Arthritis	1.31	FN1, JUN, FZD4, CCND1
Regulation of the Epithelial-Mesenchymal Transition Pathway	1.28	FZD4, FGFRL1, FGF19
ILK Signaling	1.24	FN1, JUN, CCND1
IL-8 Signaling	1.24	HMOX1, JUN, CCND1
Molecular Mechanisms of Cancer	0.633	JUN, FZD4, CCND1

DEGs were involved in inhibiting cell proliferation in GSK-J4-treated EBs

For a more in-depth analysis, we performed IPA to determine the cellular and molecular functions of all DEGs and identified 19 enriched categories (FDR < 0.05, a minimum of 5 genes represented). The top 10 categories are plotted in Fig. 2A. The top ranked functional category was “cellular growth and proliferation”, in which 60 DEGs are listed. The enrichment of 57.14% (60 out of 105) of all DEGs suggested that GSK-J4 was subjected to cell proliferation regulation. A total of 57 out of 60 genes were directly associated with cell proliferation. The IPA provided a network for these 57 DEGs showing that they were associated with the inhibition of cell proliferation (Fig. 2B). To elucidate the related expression patterns of these DEGs, we used the UCSC genome browser. We randomly selected 6 genes, 3 up-regulated genes and 3 down-regulated genes, and the exported graphs are presented in Fig. 2C. The expression of the selected genes was differentially altered by GSK-J4 as clearly plotted in Fig. 2C. Thus, it was suggested that a majority of the altered genes were involved in the inhibition of cell proliferation. Thus, we assessed whether GSK-J4 also inhibited the proliferation of EBs in vitro. We performed a cell proliferation assay (WST-1 assay) of EB+RA and EB+RA+GSK-J4 samples and observed that GSK-J4 significantly inhibited the overall proliferation of EBs (Fig. 2D). The results suggested that GSK-J4 inhibited cell proliferation by altering the expression level of related genes in differentiating EBs.

Fig. 2

Inhibition of cell proliferation by GSK-J4 in differentiating EBs

(A) Functional annotation of all DEGs with respect to cellular and molecular functions. The top 10 significant categories are plotted. Related p-values are indicated beside each block. Numbers indicated inside each block represent enriched gene numbers in each category. (B) The network for the “proliferation of cells” category adopted from IPA. The relationship is indicated as dotted lines. The blue-colored dotted lines represent the inhibition of cellular proliferation. Molecule legends are presented at the bottom of the network. The up- and down-regulated genes are indicated in orange and green, respectively. (C) Validation of related gene expression using the UCSC genome browser. We plotted the relative gene expression of randomly selected genes (BTG2, TXNIP, CDKN1C, GBX2, GPR3 and FZD4), where different colors represent different samples. d) Effect of 10 μM GSK-J4 treatment on the proliferating cell number. Proliferating cells were significantly reduced after GSK-J4 treatment. The values are represented as the average proliferating cells ± SEM bars, n = 3 replicates. Asterisks indicate statistically significant changes based on adjusted p-values < 0.05.

Inhibition of cell cycle progression was observed in GSK-J4-treated EBs

As shown in Fig. 2A, a significant number of genes were enriched in the “cell cycle” category. A total of 27 out of 105 DEGs (25.71%) were listed in this category. Thus, we wanted to identify the relationship between the cell cycle and inhibition of cell proliferation. We analyzed the 27 DEGs and observed that 24 of the 27 genes were associated with cell cycle progression. The IPA provided a network for these 24 DEGs, showing an association with the overall progression of the cell cycle (Fig. 3A). To view the related expression patterns of these DEGs, we used the UCSC genome browser. We randomly selected 6 genes, 3 up-regulated genes and 3 down-regulated genes, and the exported graphs are presented in Fig. 3B. Expression of the selected genes was differentially altered by GSK-J4, as clearly plotted in Fig. 3B. At this stage, we hypothesized that the inhibition of cell proliferation by GSK-J4 resulted from the alteration of cell cycle dynamics. To verify whether GSK-J4-mediated alterations in cell cycle-related gene expression altered the overall cell cycle dynamics, we performed a cell cycle assay for both control and GSK-J4-treated EBs. We observed that GSK-J4 significantly altered the overall cell counts in each phase of the cell cycle (Fig. 3C). We observed that the cells accumulated in the S and G2 phases with a reduced fraction in the G1 phase in GSK-J4-treated EBs. These results suggested that GSK-J4 mediated alterations in cell cycle-related gene expression, thereby inducing alterations in overall cell cycle dynamics. Furthermore, this suppression of the cell cycle might have consequences for GSK-J4-mediated inhibition of cell proliferation.

Fig. 3

Suppression of cell cycle progression by GSK-J4 in differentiating EBs

(A) Network for the category “cell cycle progression” adopted from IPA. The relationship is indicated as dotted lines. The blue-colored dotted lines represent the inhibition of cell cycle progression. Molecule legends are presented at the bottom of the network. The up- and down-regulated genes are indicated in orange and green, respectively. (B)Validation of related gene expression using the UCSC genome browser. We plotted the relative gene expression of randomly selected genes (GADD45B, GADD45G, HMOX1, ID3, CCND1 and LIN28A), where different colors represent different samples. (C) Cell cycle analysis of control and demethylase-inhibited samples. i) and ii) show the results of the cell cycle analysis using flow cytometry in EB+RA and EB+RA+GSK-J4, respectively. The peaks in the illustration correspond to the G1, S and G2 phases of the cell cycle. iii) Bar graphs showing the percentages of cells in each phase of the cell cycle. The asterisks indicate statistically significant changes based on adjusted p-values < 0.05.

Cell death was induced in EBs after GSK-J4 exposure

We also observed that 45.71% of all DEGs (48 out of 105 genes) were listed in the category of “cell death and survival”, as shown in Fig. 2A. We had expected to observe a relationship between all DEGs. After analyzing these 48 genes, we observed that the induction of cell death was mediated by apoptosis and necrosis. Among the 48 DEGs enriched in “cell death and survival” 41 and 40 genes were listed under “apoptosis” and “necrosis”, respectively. These two categories showed similar genes (almost 90% in common). The IPA provided a network for these DEGs demonstrating that GSK-J4-treated EBs showed induction of apoptosis and necrosis (Fig. 4A). To determine the related expression patterns of these DEGs, we used the UCSC genome browser. We randomly selected 8 genes, 4 up-regulated genes and 4 downregulated genes, and the exported graphs are presented in Fig. 4B. The expression of the selected genes was differentially altered by GSK-J4, as clearly plotted in Fig. 4B. To observe whether GSK-J4-mediated alterations of apoptosis and necrosis-related genes lead to overall cell death, we performed a cellular apoptosis assay. Using an Annexin V-FITC apoptosis kit, we observed that GSK-J4 treatment resulted in a 2.23% increase in cellular necrosis, whereas cellular apoptosis was achieved in 0.39% cells with overall death of 3.13% of the cells, and these results were statistically significant (Fig. 4C). These results suggested that GSK-J4-mediated alterations of apoptosis and necrosis-related genes induced cell death.

Fig. 4

GSK-J4 treatment induces cell death through apoptosis and necrosis in differentiating EBs

(A) Network for the categories “apoptosis” and “necrosis” adopted from IPA. The relationship is indicated using dotted lines. The mastered colored dotted lines represent the activation of cellular death either through apoptosis or necrosis. Molecule legends are presented at the bottom of the networks. The up-and down-regulated genes are indicated in orange and green, respectively. (B) Validation of related gene expression using the UCSC genome browser. The relative gene expression of randomly selected genes (JUN, MSX2, ULK1, PPP1R15A, HOXA1, HOXB1, MiR17HG and VASH2) is plotted, where different colors represent different samples. c) Cellular apoptosis was measured using FACS and an Annexin V-FITC apoptosis kit. A two-parameter histogram dot plot displays FITC on the x-axis and PI on the y-axis. i) and ii) represent the percentage of cells that were positive for Annexin V-FITC and/or propidium iodide inside the quadrants (Q1 = necrotic cells, Q2 = late-stage apoptotic cells, Q3 = living cells, and Q4 = early-stage apoptotic cells). iii) shows the percentage of total cell death in the EB+RA and EB+RA+GSK-J4 samples. The asterisks indicate statistically significant changes based on adjusted p-values < 0.05.

Profiling of upstream transcriptional regulators of all DEGs altered by GSK-J4 in differentiating EBs

Next, we searched for a set of master regulators responsible for inhibiting cell proliferation. We identified four common genes in all three constructed networks: cyclin D1 (CCND1), growth arrest and DNA damage-inducible beta (GADD45B), growth arrest and DNA damage inducible gamma (GADD45G) and miR-17-92a-1 cluster host gene (MIR17HG). Consistent with previously published reports, we proposed that these 4 genes play an important role in the GSK-J4-mediated inhibition of cell proliferation through cell cycle suppression and activation of cell death. A network was exported from the IPA of these four genes and their upstream transcriptional regulators to determine whether v-myc avian myelocytomatosis viral oncogene homolog (MYC) was a master regulator in GSK-J4-treated EBs (Fig. 5A). We also listed the upstream transcriptional regulators of all 105 DEGs. We exported a list from IPA, identifying 49 candidates as the upstream regulator. The detailed list is presented in Table 2.

Fig. 5

Motif enrichment analysis

(A) Network showing up-streaming transcriptional regulators of CCND1, GADD45B, GADD45G and MIR17HG. (B) and (C) represent transcription factor binding motifs of all up- and down-regulated genes, respectively. We performed enrichment analysis using the known TF motifs in the JASPAR database. The top 5 motifs are plotted for up- and down-regulated genes as indicated.

Table 2

List of upstream transcriptional regulators of all differentially expressed genes

Upstream transcriptional regulator	p-value of overlap	Target genes in present dataset
CREB1	8.27E-10	BTG2, CCND1, FGF19, FN1, GADD45B, GADD45G, GAP43, GPR3, HECA, HMOX1, HSPA5, JUN, MVD, PPP1R15A, UPP1, YPEL4
CTNNB1	2.42E-08	ANXA1, CCND1, ETV1, ETV4, FN1, GADD45B, GADD45G, GAP43, HILPDA, ID3, JUN, LBH, MSX2, RASSF4, SERPINA5, TNC
SMAD3	2.95E-08	CCND1, CDH5, FN1, GADD45B, HMOX1, ID3, JUN, LEFTY2, PMEPA1, TNC
FOXO1	1.94E-07	CCND1, CDKN1C, FN1, GADD45B, HMOX1, HSPA5, JUN, LEFTY2, MVD, PPP1R15A, TXNIP
NFKBIA	2.5E-07	ACKR3, BCL11A, BTG2, CCND1, FN1, GADD45B, GADD45G, HMOX1, HSPA5, JUN, PMEPA1, TNC
SOX2	7.45E-07	CCND1, DPPA4, ETV4, ETV5, GADD45B, GBX2, HOXB1, JUN, LIN28A, MSX2
KLF4	8.99E-07	CCND1, CDKN1C, ETV1, FN1, GBX2, ID3, LIN28A, MSX2, TNC
NFYA	1.21E-06	GADD45B, HSPA5, JUN, PPP1R15A, PTTG1, TXNIP
TP63	2.75E-06	CCND1, CDKN1C, FN1, HAS3, ID3, LIN28A, MIR205HG, TNC, ULK1, UPK1A
TCF4	9.56E-06	CCND1, CDKN1C, FN1, JUN, MSX2
JUN	1.15E-05	ANXA1, CCND1, FN1, GAP43, HMOX1, JUN, MVD, NBPF10, PTBP2, TNC
RBPJ	2.46E-05	CCND1, ECT2, JUN, LEFTY2, MIB2, TNC
SMAD2	2.53E-05	CDH5, FN1, GADD45B, HMOX1, LEFTY2
SMAD4	3.19E-05	CCND1, CDH5, FN1, GADD45B, HMOX1, MSX2, TNC
TP53	0.000059	ANXA1, BTG2, CCND1, FN1, GADD45G, GBX2, HMGB1, HMOX1, ID3, JUN, MIR17HG, MVD, PMEPA1, PPP1R15A, PTTG1, ULK1, UPP1
EP300	6.64E-05	CCND1, FN1, HMOX1, ID3, JUN, LEFTY2, TNC, TXNIP, YPEL2
JUNB	0.000165	CCND1, FN1, HMOX1, MVD, PTBP2
HIF1A	0.000177	CCND1, FN1, GADD45B, HILPDA, HMOX1, HSPA5, JUN, P4HA2
VHL	0.000185	CCND1, CDKN1C, FN1, HMOX1, NEDD9
TCF3	0.00032	CDH5, ECT2, GADD45B, JUN, RASSF4, SCN4A
CEBPA	0.000355	ANXA1, BTG2, EEF1A2, GBX2, HMOX1, HSPA5, JUN, MSX2
POU5F1	0.000365	CCND1, GBX2, HOXB1, ID3, LEFTY2, LIN28A, PMEPA1
MYC	0.000381	CCND1, FBN2, FN1, GADD45B, GADD45G, GBX2, HMOX1, ID3, JUN, MIR17HG, MSX2, TNC, TXNIP
TWIST1	0.000402	FBN2, FN1, LEFTY2, MRC2, OLFML3
REL	0.000406	CCND1, ETV5, FN1, GADD45B, JUN, PPP1R15A
FOS	0.000413	CCND1, FN1, HMOX1, HOXA1, HSPA5, JUN, NBPF10, PMEPA1, PTTG1
YY1	0.00063	BTG2, GADD45G, HMOX1, HSPA5, MSX2, PPP1R15A
FOXO3	0.000772	CCND1, CDKN1C, GADD45B, PPP1R15A, TXNIP
EGR1	0.000834	CCND1, FN1, GADD45B, HMOX1, JUN
PDX1	0.000834	CCND1, GAP43, ID3, JUN, TXNIP
MYOD1	0.0011	CDKN1C, ETV4, GADD45G, PMEPA1, SCN4A
CDKN2A	0.00117	BTG2, CCND1, GADD45G, JUN, P4HA2, PMEPA1
EPAS1	0.00135	CCND1, FN1, GADD45B, HILPDA, HSPA5
PAX3	0.00141	ANXA1, ID3, MSX2, NEDD9, TNC
SP3	0.00195	CCND1, HSPA5, JUN, PTTG1, TNC
ETS1	0.00199	BCL11A, CCND1, CDH5, FN1, HMOX1
E2F4	0.00247	CCND1, ECT2, HMGB1, ID3, PTTG1
TP73	0.00255	BTG2, CDKN1C, LBH, SERPINA5, SPINT1, ULK1
HTT	0.00258	ATP5G1, CCND1, ETV1, FN1, GADD45G, GAP43, HOXB1, HSPA5, JUN
NUPR1	0.00519	ETV1, HILPDA, HIST2H3D, P4HA2, PPP1R15A, UPP1, ZNF483
MYCN	0.00555	CCND1, FN1, MIR17HG, MRC2, MSX2
SP1	0.0074	CCND1, FN1, HMOX1, HSPA5, JUN, PTTG1, TNC
TCF7L2	0.0109	ANXA1, CCND1, ID3, LBH, YPEL2
E2F1	0.0122	ATP5G1, CCND1, CDKN1C, ECT2, HMGB1, ID3
SRF	0.0124	ETV1, GADD45G, JUN, LBH, TNC
STAT3	0.0212	CCND1, CDH5, ECT2, FN1, GADD45G, HMOX1
CREBBP	0.0223	CCND1, CDH5, JUN, TNC, YPEL2
RELA	0.0223	BTG2, CCND1, FN1, HMOX1, JUN
SMARCA4	0.0345	CCND1, FN1, GADD45G, ID3, JUN, LRAT

In Fig. 2A, we observed an enrichment of the “gene expression” category, prompting us to profile the TF-encoding genes among all 105 DEGs in the dataset. To identify GSKJ4-targeted downstream TFs, we compared the DEGs identified in the present study with a list of previously identified human TFs (Vaquerizas et al., 2009). We identified a total of 16 TF-encoding genes among all DEGs, which are listed in Table 3. We also applied a TF motif analysis to observe the functional association of TFs in GSK-J4-mediated gene expression in RA-treated EBs. The Pscan software tool was used to perform the in silico computational analysis of over-represented cis-regulatory elements within the 5′-promoter regions of coordinately regulated genes using the JASPAR database. Applying this score to the promoters (from −950 bp to +50 bp) of the differentially expressed genes revealed that the putative binding sites for different motifs were significantly enriched. We identified a total of 13 TFs with binding sites that were significantly over-represented in the promoter region of the up-regulated genes (p < 0.001). Only 9 TFs were down-regulated genes with binding sites that were significantly over-represented in the promoter regions (p < 0.01). The top 5 enriched motifs of both up- and down-regulated genes are presented in Figs. 5B and 5C, respectively. Among the up-regulated genes, the top ranked motif was E2F3, and among down-regulated genes, it was ZNF740. Furthermore, we found that 44 of 47 up-regulated genes contained an E2F3 binding motif in the promoter sequence (from −950 bp to +50 bp) (Supplementary Table 2). In the case of ZNF740, we observed that 55 out of 58 down-regulated gene promoters contained binding motifs (Supplementary Table 3). These findings suggested that E2F3 and ZNF740 had distinct role in GSK-J4-mediated gene transcription. Similar findings were observed for other enriched TF binding motifs that need further functional experiments for verification. We also attempted to identify common TFs between the 16 TFs among all DEGs and 22 TFs resulting from promoter motif analysis. We observed no matches as expected. The expression of 16 TFs was determined at in the mRNA level. These TFs will likely alter the cell fate when expressed as proteins. In contrast, the 22 TFs resulting from promoter motif analysis showed an important role in regulating the expression of the listed DEGs.

Table 3

List of transcription factor-encoding genes in the data set

Gene symbol	Description	Log2 FC
BCL11A	B-cell CLL/lymphoma 11a	1.633008456
JUN	Jun proto-oncogene, ap-1 transcription factor subunit	1.585211074
MSX2	Msh homeobox 2	1.409957904
ZNF483	Zinc finger protein 483	−2.003457678
HOXB1	Homeobox b1	−1.894954643
ID3	Inhibitor of DNA binding 3, HLH protein	−1.760864137
GBX2	Gastrulation brain homeobox 2	−1.707942582
ETV5	ETS variant 5	−1.704617738
ATOH8	Atonal BHLH transcription factor 8	−1.697759027
OLIG3	Oligodendrocyte transcription factor 3	−1.565186008
ETV4	ETS variant 4	−1.554508961
NFE2L3	Nuclear factor, erythroid 2 like 3	−1.425514822
HOXA1	Homeobox a1	−1.411087323
HMGB1	High mobility group box 1	−1.296979403
LARP7	La ribonucleoprotein domain family member 7	−1.270947331
ETV1	ETS variant 1	−1.206732073

The DEGs found in GSK-J4-mtreated EBs were unique and are not shared with the RA-treated data set

Our experimental model is an early differentiating model of neural differentiation. We exposed EBs to RA to induce differentiation and observed GSK-J4-mediated effects during differentiation. The results showed that GSK-J4-altered genes that were significantly enriched in cell proliferation and cell death. These two cellular functions are basic features of a differentiating cell. To ascertain whether RA treatment also shared similar gene expression patterns during differentiation, we performed RNA-seq analysis for EB vs EB+RA. Genes that were differentially expressed included 51 and 9 up- and down-regulated genes in RA-treated EBs respectively, with a cut-off value of 1.2 log2-fold changes (p-value < 0.05) (Fig. 6A). We compared RA-mediated and GSK-J4-mediated DEGs to select shared genes. The results showed no shared genes between these two data sets (Fig. 6A). We also performed gene ontology analysis for the DEGs found in RA-treated EBs using the DAVID analysis tool. We identified 5 categories (biological processes) with up-regulated genes (threshold count: 6, EASE score: 0.05 with Benjamini correction). In contrast, no category was enriched with down-regulated genes. We observed that the genes that were up-regulated in GSK-J4-treated cells were largely involved in multicellular organism development and regulation of transcription, among others (Fig. 6B). Additionally, the pathway analysis using IPA revealed 4 signaling pathways (FDR < 0.05, a minimum of 3 genes represented). The top-ranked significantly enriched signaling pathway was RAR Activation. A complete list of the enriched pathways and related genes is presented in Supplementary Table 4. We also determined the cellular and molecular functions of all DEGs using IPA and identified 13 enriched categories (FDR < 0.05, a minimum of 5 genes represented). The top 10 categories are plotted in Fig. 6C. The enriched functions also contained categories such as cellular proliferation, gene expression and cell survival, similar to Fig. 2A but with a different list of genes. These results suggest that the DEGs found in GSK-J4-mtreated EBs were unique and were not share with the RA-treated data set.

Fig. 6

Functional annotation of differential gene expression in RA-treated EBs

(A) The intersection between (EB vs EB+RA) and (EB+RA vs EB+RA+GSK-J4) showed that no genes were shared. (B) Bar graphs from the gene ontology analysis (Biological Process) of up- regulated genes. Five enriched categories are plotted with their p-values indicated in front of each bar. (C) Functional annotation of all DEGs with respect to cellular and molecular functions. The top 10 significant categories are plotted. Related p-values are indicated beside each block. Numbers indicated inside each block represent enriched gene numbers in each category.

GSK-J4 treatment facilitated the neural differentiation process

To assess the methylation status in GSK-J4-treated differentiating EBs, we performed western blotting for H3K27me3 and observed elevated expression (Fig. 7A). This result suggested potential consequences of the GSK-J4-mediated change in methylation status and DEGs in the RNA-seq analysis. Further functional analyses are needed to validate this conjecture. We also performed qRT-PCR analysis for 10 DEGs to observe their relative expression. The results showed that GSK-J4 treatment of differentiating EBs significantly altered the expression of important genes (Fig. 7B).

Fig. 7

qRT-PCR analysis to validate the RNA-seq results and neural marker expression

(A) Changes in H3K27me3 levels in response to GSK-J4 during EB differentiation. Whole cell extracts were collected from cells treated for 48 h, isolated using RIPA buffer, and immunoblotted with H3K27me3 and H3 antibodies. ER and ERG represents EB+RA and EB+RA+GSK-J4, respectively. (B) qRT-PCR results. We randomly selected 10 genes to observe their expression levels in GSK-J4-treated differentiating EBs. The expression value was normalized to the GAPDH expression level. Values are represented as the average mRNA expression ± SEM bars, n = 3 replicates. Asterisks indicate statistically significant changes based on adjusted p-values < 0.05. (C) qRT-PCR result showing the relative expression of neural marker genes in GSK-J4-treated differentiating EBs. (D) Schematic presentation of GSK-J4-mediated inhibition of cell proliferation through the suppression of cell cycle and activation of cell death.

Cell proliferation and differentiation show a remarkable inverse relationship (Ruijtenberg and van den Heuvel, 2016). Precursor cells continue to divide before acquiring a fully differentiated state, while terminal differentiation usually coincides with proliferation arrest and permanent exit from the division cycle. In our experimental model, we observed a rapid reduction of cell proliferation status due to the exposure to GSK-J4. To check whether GSK-J4 favored differentiation, we performed qRT-PCR analysis to assess the known neural precursor markers NESTIN and PAX6. The results showed that both markers were significantly induced compared with the RA-treated control (Fig. 7C).

DISCUSSION

Epigenetic modifications and gene expression regulation are key events during cellular differentiation. Different types of histone modifications have been previously implicated in cell fate determination. Histone modification by lysine methylation has been directly associated with epigenetic inheritance. The critical roles of histone modification have not yet been well characterized. Common methods to reveal the molecular mechanisms of histone modifications are the use of knock-down/knock-out techniques and/or of an inhibitor that disables the catalytic activities of related enzymes. RNA interference and site-specific mutagenesis are readily used, but these techniques may affect the integrity of the demethylases (Burgold et al., 2008; Morales Torres et al., 2013). From this aspect, a chemical inhibitor that binds to the active site and prevents the interaction between the enzyme and the substrate without changing the conformation of the enzymes may be the best solution. GSK-J4 is a widely used inhibitor of JMJD3/UTX-selective histone demethylases (Ferreri et al., 2010). Kruidenier et al. (2012) first reported the use of GSK-J4 in a study observing the role of histone demethylase in macrophages. Subsequently, several reports using different cell types were published in which the functions of JMJD3/UTX were the main focus (Ferreri et al., 2010; Heinemann et al., 2014; Kamikawa and Donohoe, 2015; Kang et al., 2015). A few studies have reported the association of GSK-J4 with the suppression of embryonic development, (Xie et al., 2016) cellular differentiation (Liu et al., 2015b) and increased apoptosis (Tennant et al., 2015). A recent study reported that the inhibition of histone demethylase by GSK-J4 affected cellular proliferation and cell cycle progression in differentiating cells (Hofstetter et al., 2016; Sakaki et al., 2015; Watarai et al., 2016), but no reports have described the genes responsible for these functions. Here, we performed a transcriptomics analysis to profile GSK-J4-targeted genes during early neural differentiation and provided a list of important genes that induce cell death with cell cycle suppression and inhibition of cell proliferation. Thus, the present study is novel in this respect.

After a critical literature review of all DEGs, we identified four candidate genes CCND1, GADD45B, GDD45G and MIR17HG, which play critical roles in the inhibition of cell proliferation, cell cycle suppression and cell death. CCND1 is a regulatory subunit of the cyclin-dependent kinases CDK4/6, which are positive regulators of cell proliferation (Dai et al., 2016). The down-regulation of CCND1 gene expression induced the inhibition of cell proliferation. A recent study showed that the inhibition of CCND1 suppresses cell proliferation in NSCLC cell lines (Li et al., 2016). Several studies have reported the inhibition of cell proliferation, (Chen et al., 2013; Zhu et al., 2010) suppression of cell cycle progression (Chen et al., 2016; Kapral et al., 2011) and increased apoptosis (Dai et al., 2016; Lai et al., 2010; Wang et al., 2014). In the present study, we also observed that GSK-J4 induced the suppression of CCND1 gene expression and might be responsible for the suppression of cell cycle progression, followed by the inhibition of cell proliferation and increased cell death.

Among the potent growth-arrest and DNA-damage-inducible genes, GADD45 family genes are well known. GADD45B contributes to the induction of apoptosis and cell cycle alterations in response to several stimuli. Previous studies have reported that the inhibition of histone deacetylase induces GADD45B and suppresses cell cycle progression (Greenberg et al., 2001; Nettersheim et al., 2016). Although the distinct function of GADD45B remains unclear, this gene has been shown to inhibit cell proliferation, (Kageyama et al., 2015; Liu et al., 2005) alteration of the cell cycle (Mak and Kultz, 2004; Vairapandi et al., 2002; Wang et al., 1999) and apoptosis (Cho et al., 2010; Liu et al., 2015a). However, GADD45G also induces similar functions regarding cell cycle regulation and cell death. Studies have implicated the induction of GADD45G in the inhibition of cell proliferation, (Flores and Burnstein, 2010; Ishida et al., 2013) suppression of cell cycle progression (Guo et al., 2013; Vairapandi et al., 2002; Zhao et al., 2014) and activation of apoptosis (Salvador et al., 2013; Zhao et al., 2014). GSK-J4 treatment also induces GADD45B/G expression in differentiating EBs, and this elevated expression might be involved in regulating cell proliferation, cell cycle and apoptosis.

MiRNA17HG targets multiple cellular pathways, favoring the proliferation of cells and inhibiting apoptosis. A previous study suggested that miRNA17HG increases cell proliferation, (Sandhu et al., 2013; Zhang et al., 2013) by targeting p21 and sequentially activating the cyclinD1-CDK4 complex (Jiang et al., 2010). miRNA-17-92 also minimizes MyC-induced apoptosis by inducing BCL2 (Xiao et al., 2008). We observed that miRNA17HG is also involved in cell cycle regulation (Yang et al., 2014). The down-regulation of miRNA17HG gene regulation might induce cell proliferation and activate apoptosis. In the present study, we observed similar incidents with GSK-J4 treatment. The complex functions of miRNA-17-92 suggested the presence of additional targets for regulating the cellular fate in GSK-J4-treated EBs, warranting further attention in future studies.

The mechanisms that maintain the balance between cell proliferation and differentiation are often compromised. Generally, the proliferation rate is reduced with the progression of differentiation. In our model, GSK-J4 caused a rapid reduction of cell proliferation. The induction of neural markers suggested that there might be a correlation of GSK-J4-mediated suppression of the cell proliferation rate and the differentiation progress; however, further functional analyses are required to validate this relationship.

In summary, the transcriptomics analysis and subsequent functional study revealed a significant number of DEGs involved in the inhibition of cell proliferation, suppression of cell cycle progression and activation of cellular death in GSK-J4-treated EBs. We also validated the predicted impact on differentiating EBs. However, we did not perform any additional experiments to determine whether the listed genes were epigenetically altered (directly or indirectly). It would be interesting to examine this topic in future studies. Here, we provided a platform to determine the molecular mechanisms by which GSK-J4 mediates the inhibition of cell proliferation, suppression of cell cycle progression and activation of cellular apoptosis and/or necrosis (Fig. 7D).

CONCLUSION

In conclusion, we treated differentiating cells with GSK-J4 to profile target genes and observed that GSK-J4 inhibited cell proliferation, suppressed cell cycle progression and activated cell death by apoptosis/necrosis. Transcriptomics expression of the listed genes was highlighted in the present study. These findings will expand the current understanding of the biology of histone-modifying enzymes. Furthermore, we propose that elucidating the mechanism underlying the alteration of GSK-J4-target genes may provide useful information for the development of drugs, particularly for cancer treatment.