Low-dose DNA demethylating therapy induces reprogramming of diverse cancer-related pathways at the single-cell level

Epigenetic reprogramming using DNA demethylating drugs is a promising approach for cancer therapy, but its efficacy is highly dependent on the dosing regimen. Low-dose treatment for a prolonged period shows a remarkable therapeutic efficacy, despite its small demethylating effect. Here, we aimed to explore the mechanisms of how such low-dose treatment shows this remarkable efficacy by focusing on epigenetic reprograming at the single-cell level. Expression profiles in HCT116 cells treated with decitabine (DAC) were analyzed by single-cell RNA-sequencing (scRNA-seq). Functional consequences and DNA demethylation at the single-cell level were analyzed using cloned HCT116 cells after DAC treatment. scRNA-seq revealed that DAC-treated cells had highly diverse expression profiles at the single-cell level, and tumor-suppressor genes, endogenous retroviruses, and interferon-stimulated genes were upregulated in random fractions of cells. DNA methylation analysis of cloned HCT116 cells revealed that, while only partial reduction of DNA methylation levels was observed in bulk cells, complete demethylation of specific cancer-related genes, such as cell cycle regulation, WNT pathway, p53 pathway, and TGF-β pathway, was observed, depending upon clones. Functionally, a clone with complete demethylation of CDKN2A (p16) had a larger fraction of cells with tetraploid than parental cells, indicating induction of cellular senescence due to normalization of cell cycle regulation. Epigenetic reprogramming of specific cancer-related pathways at the single-cell level is likely to underlie the remarkable efficacy of low-dose DNA demethylating therapy.

Early clinical studies used DNA demethylating drugs with the maximum tolerated dose (MTD) aiming at cytotoxic effects, and such studies showed only limited therapeutic efficacy with strong toxicity [12]. Afterwards, administration with a dose much lower than the MTD for a prolonged period was shown to have high therapeutic efficacy with limited toxicity [13][14][15]. The therapeutic effect was first explained by re-activation of a specific tumor-suppressor gene, such as CDKN2A [16], and was then shown to be associated with the suppression of tumor-initiating cells by restoration of multiple pathways in tumor cells [17]. In addition, enhancement of antigenicity of tumor cells by activation of endogenous retroviruses [18,19] was found to be an important mode of action. Recently, in addition to the effect on tumor cells, that on tumor cell niche, including cancerassociated fibroblasts and myeloid-derived suppressor cells (MDSCs) has been suggested also to be involved [20][21][22].
Despite the remarkable therapeutic efficacy of lowdose and prolonged treatment with reprograming of multiple target genes, one remaining question is why only partial demethylation of the target genes [15,17] can exert such high therapeutic efficacy. Considering that cells have two alleles for most genes, it is expected that, at the single-cell level, demethylation of a specific gene should be complete, 50%, or none. In this study, we aimed to explore whether complete demethylation of specific genes is really induced at the single-cell level and to analyze the functional consequences of such complete demethylation of specific genes.

DAC-treated single cells had highly diverse expression profiles
Single cell RNA sequencing (scRNA-seq) was conducted using 1783 mock-treated and 1751 DACtreated HCT116 cells (Fig. 1a). On average, expression of 4867 and 5838 genes per cell was detected in mock-and DAC-treated cells, respectively. Uniform Manifold Approximation and Projection (UMAP) analysis was conducted using 14,099 genes that can be induced by DAC treatment (UMI counts ≤ 2 in all the 1783 mock-treated cells). It was shown that expression profiles in DAC-treated cells had high diversity (Fig. 1b). Hierarchical clustering analysis was conducted using highly upregulated genes (top 200 genes with higher mean UMI counts in DAC-treated single cells) selected from the 14,099 genes. It was shown that genes with higher expression levels were different, depending upon DAC-treated clones (Fig.  1c). Among the 1751 DAC-treated single cells, random fractions of cells showed upregulation (UMI counts ≥ 3) of specific established tumor-suppressor genes methylation-silenced in colorectal cancers [23][24][25] (CHFR, FBLN2, ICAM4, and SFRP1), endogenous retroviruses (ERVs) (ERVMER61-1 and ERVW-1), and interferon-stimulated genes (IFI16 and IFI44) (Fig. 1d, e, f). These results showed that DAC-treated single cells had diverse expression profiles.

DAC-treated single cells exhibited complete DNA demethylation of specific CGIs
To analyze the functional consequences of the diverse expression profiles and DNA methylation at the single-cell level, DAC-treated HCT116 cells were cloned by the limiting dilution method, and nine clones were obtained (Fig. 2a). Using these nine clones and their parental cells in bulk, DNA methylation was analyzed in a genome-wide manner, and methylation levels of 270,249 genomic blocks in DAC-treated bulk cells and clones were compared with those in mock-treated bulk cells. In the bulk cells, only partial reduction of DNA methylation levels was induced by DAC treatment (Fig. 2b, Fig.  S1), and none of the 1039 TSS200CGIs (200 bp upstream regions from TSSs with CGIs) fully methylated in HCT116 cells were completely demethylated (β < 0.2). In contrast, in all of the clones, methylation levels were markedly reduced by DAC treatment (Fig.  2c, Fig. S2), and 126-293 of the 1039 TSS200CGIs were completely demethylated. These results showed that DAC-treated single cells exhibited complete demethylation of a large number of CGIs.
Demethylated genes were highly diverse among the DAC-treated clones To examine whether similar or different sets of genes were demethylated in different clones, DNA methylation levels of all the 270,249 genomic blocks in individual clones were compared. A large fraction of blocks showed different methylation levels among clones (Fig. 3a). DNA methylation status was then compared among clones for 1039 TSS200CGIs fully methylated in HCT116. Completely demethylated TSS200CGIs were highly diverse among the nine clones (Fig. 3b). Only 15.2-54.1% of completely demethylated TSS200CGIs overlapped (Table S1), and 700 (67.4%) of the 1039 TSS200CGIs were completely demethylated in at least one of the nine clones. These results showed that completely demethylated genes were highly diverse among the DAC-treated clones.
We next focused on genes involved in colon cancerrelated pathways, namely cell cycle regulation, WNT pathway, p53 pathway, and TGF-β pathway. Regarding the cell cycle regulation, CDKN2A and CHFR were partially and fully methylated in HCT116, respectively. CDKN2A and CHFR were completely demethylated in one and three of the nine clones, respectively (Fig. 3c). Regarding the WNT pathway, six of its negative regulators, DKK3, SFRP1, SFRP2, SFRP5, SOX17, and WIF1, were fully methylated in HCT116. None of the six genes were demethylated in the DAC-treated bulk cells, but five of the six genes were completely demethylated in 1-3 of the nine clones (Fig. 3c). Genes involved in the p53 pathway, IGFBP7 and MIR34B, and the TGF-β pathway, BMP2, BMP3, BMP6, and BMPR1B, were also completely demethylated, depending upon clones. It was confirmed that complete demethylation led to the upregulation of methylation-silenced genes at the protein level (Fig. 3d). These results showed that different sets of genes in different cancer-related pathways were completely demethylated at the single-cell level. Changes of cell characteristics were first analyzed by examining cell proliferation of the nine clones along with mock-treated and DAC-treated bulk cells. The growth rates of the nine clones were slower than that of their parental cells, but the decreases were different among the clones (51.4-80.2%) (Fig. 4a). The doubling time of the individual clone was not associated with the number of completely demethylated genes (Fig. S3). These results suggested the diverse sets of re-activated genes, rather than the number of completely demethylated genes, led to diverse functional consequences.
To confirm that complete demethylation of a specific gene in a clone actually resulted in normalization of a specific cancer-related pathway, cell cycles in the DAC-treated clones, along with the DAC-treated bulk cells, were compared to that of mock-treated bulk cells. The H3-32 clone, which exhibited complete demethylation of CDKN2A, had a fraction of tetraploid (> 4 N) cells much larger than the mock-treated bulk cells (p = 0.0047) (Fig. 4b, c). In contrast, the seven of the other eight clones that did not exhibit complete demethylation, along with the DAC-treated bulk cells, did not have an increased fraction of tetraploid cells compared to the mock-treated bulk cells. Since expression of p16, encoded by CDKN2A, is known to play a key role in the induction of cellular senescence [26,27], the increased fraction of tetraploid cells was considered to be due to the induction of cellular senescence by CDKN2A re-activation. This result showed that complete demethylation of a specific gene actually induced normalization of a specific cancer-related pathway. Fig. 3 Difference of demethylated genes between DAC-treated clones. a Comparison of DNA methylation levels between DAC-treated clones. A large fraction of genomic blocks showed different methylation levels between the two clones compared. b Unsupervised hierarchical clustering analysis using DNA methylation profiles of 1039 TSS200CGIs fully methylated in HCT116. Demethylated genes were highly diverse among the DAC-treated clones. c DNA methylation status of genes involved in colon cancer-related pathways. Genes involved in cell cycle regulation, the WNT pathway, the p53 pathway, and the TGF-β pathway were completely demethylated, depending upon clones. d Protein expression levels of DKK3 and SFRP1 in DAC-treated clones. Complete demethylation led to the upregulation of methylation-silenced genes at the protein level in DAC-treated clones

Discussion
Single-cell analysis revealed that diverse expression profiles were produced by a low-dose treatment with a DNA demethylating agent, here DAC. Different sets of genes in different cancer-related pathways were completely demethylated in individual cells, and such complete demethylation induced normalization of a specific cancer-related pathway, such as regulation of cell cycle. The finding solved the long-standing question, at least in part, of why low-dose DAC treatment can have strong therapeutic effects despite the fact that only partial demethylation is induced in bulk cells. Epigenetic reprogramming of different cancer-related pathways at the single-cell level, namely heterogeneous responses of individual cells to DAC, is likely to underlie the remarkable efficacy of low-dose DNA demethylating therapy (Fig. 4d). The heterogeneous response is considered to be due to the heterogeneous allelic demethylation. Based Fig. 4 Analysis of a specific cancer-related pathway in DAC-treated clones. a Cell growth rates of DAC-treated clones. The growth rates of the nine clones were slower than those of the mock-treated and DAC-treated bulk cells, but the decreases were different depending upon the clones. b Cell cycle analysis of DAC-treated H3-32 clone. Cells with tetraploid were detected in the H3-32 clone, which had complete demethylation of CDKN2A, but not in the DAC-treated bulk cells. c Fraction of cells in specific phases of cell cycle. The H3-32 clone had a fraction of tetraploid cells much larger than that of mock-treated bulk cells. Data represent mean ± SD of triplicate experiments. d The model of the mechanism underlying therapeutic efficacy of low-dose demethylating therapy. Complete demethylation of a specific gene is induced at the single-cell level by demethylating therapy, and this leads to normalization of a specific cancer-related pathway, such as cell cycle regulation, the WNT pathway, and the p53 pathway on these findings, it is considered that detection of only partial demethylation, but not complete demethylation, as a whole is sufficient to obtain therapeutic efficacy.
Usual single-cell analysis can provide vivid pictures of cellular heterogeneity even among cells that appear to be homogeneous. Here, to analyze functional consequences of the diverse expression profiles among DACtreated cells, we combined the method of single-cell cloning with single-cell RNA-seq analysis. DNA methylation analysis and phenotypic analysis using cloned cells clearly showed normalization of cancer-related pathways at the single-cell level due to complete DNA demethylation and diversity in functional consequences. Acceleration of the step of single-cell cloning by additional innovation would further advance research addressing complex groups of cells without morphological differences.
Epigenetic reprograming of diverse pathways in individual cells could theoretically lead to production of more aggressive cells. However, in a clinical setting, good response is often observed, and cells in bulk show slower growth and less aggressive phenotypes [17]. In this study, the growth rates of all the nine DAC-treated clones were indeed slower than that of mock-treated and DAC-treated bulk cells, and clones with increased growth rates did not appear. As a possible mechanism of why only clones with less aggressive phenotypes are produced despite complete demethylation of various sets of genes at the single-cell level, we can note that cancer cells have already been highly selected as the fittest for multiple environments, and it is very difficult to reprogram so that the cancer cells would become even more aggressive.
The normalization of cellular functions in diverse directions by a DNA demethylating drug suggests that epigenetic therapy would be better combined with drugs with broad actions, such as cytotoxic drugs, but not with molecular targeted drugs. Indeed, DNA demethylating drugs were frequently combined with various types of cytotoxic drugs, such as cisplatin, temozolomide, and docetaxel, in many clinical trials [6,[28][29][30], but with molecular targeted drugs, such as panitumumab and erlotinib, in only a limited number of trials [31,32]. Higher overall response has been achieved in combination therapy with cytotoxic drugs [6].
Surviving cells after DAC treatment, namely DACtreated clones, were used for analysis of normalization of specific cancer-related pathways. The surviving cells were only partially demethylated as the dead cells ( Fig.  S1) but had 7.8% less demethylation based upon the mean of the changes of the β values in all genomic blocks. The 7.8% increased demethylation in dead cells was considered to lead to an increase of completely demethylated pathways to a similar degree. However, although some bias for the normalized pathways and their numbers may have existed, the concept of normalization of different sets of pathways in different cells is likely to be valid even in the dead cells. Unfortunately, dead cells cannot be used for single cell RNA-seq analysis or analysis of cellular functions, and this could be a limitation of this study.

Conclusions
Different sets of genes in different cancer-related pathways were completely demethylated at the single-cell level, and this is likely to underlie the remarkable efficacy of low-dose DNA demethylating therapy.

Single cell RNA sequencing
A human colon cancer cell line, HCT116, was purchased from the American Type Culture Collection (Rockville, MD). The absence of Mycoplasma infection was confirmed using the MycoAlert mycoplasma detection kit (Lonza, Basel, Switzerland). The cells (1 × 10 5 ) were seeded on day 0 and treated with 0.5 μM of DAC (Sigma-Aldrich, St. Louis, MO) on days 1 and 3. The cells were placed in fresh medium without DAC on day 5 and harvested on day 7. With DAC treatment between 0.005 to 5 μM, the strongest DNA demethylation was observed with 0.5 μM of treatment (Fig. S4).
Using 1000 live cells of DAC-treated or mock-treated HCT116, a library was prepared using the Chromium Controller (10x Genomics, Pleasanton, CA) and Chromium Single Cell 3′ GEM, Library & Gel Bead Kit v3 (10x Genomics). The prepared scRNA-seq libraries were sequenced using an Illumina HiSeq platform (Illumina, San Diego, CA), and then, obtained sequence reads were de-multiplexed and mapped to the human genome (build hg38) using the Cell Ranger software (10x Genomics, version 3.0.2).

Cluster analysis of scRNA-seq data
Dimension reduction and visualization of scRNA-seq data were conducted for 14,099 genes with UMI counts of 2 or less in all the mock-treated single cells by UMAP [33] (URL; https://arxiv.org/abs/1802.03426) using the R version 3.6.2 with the umap package. Hierarchical clustering analysis was conducted using the top 200 highly upregulated genes (top 200 genes with higher mean UMI counts in DAC-treated single cells) selected from the 14,099 genes using R function hclust (https://cran.rproject.org/).

Cloning of HCT116 cells after DAC treatment
HCT116 cells (3 × 10 5 ) were seeded on day 0 and treated with 0.5 μM of DAC on days 1 and 3, and the cells were placed in a fresh medium without DAC on day 5. The cells were cloned by the limiting dilution method on day 17 (DAC clones), and clones were harvested on days 48-58.

Genome-wide DNA methylation analysis by BeadChip array
Genome-wide analysis of DNA methylation status was performed using an Infinium HumanMethylation450 BeadChip array as described [34]. The methylation level of each CpG site was represented by a β value that ranged from 0 (unmethylated) to 1 (fully methylated). All the 485,512 probes on the array were assembled into 279,922 genomic blocks based on their location against a transcription start site (TSS) and relative location against a CGI [35]. Among them, 7839 were located within TSS200CGIs. The DNA methylation level of a genomic block was evaluated using the mean β value of all the probes within the genomic block, and the methylation status of the genomic block was classified into unmethylated (β value < 0.2), partially methylated (0.2 ≤ β value < 0.8), and fully methylated (0.8 ≤ β value). Hierarchical clustering analysis was conducted for 1039 TSS200CGIs fully methylated in HCT116 using R function hclust.
DNA methylation status of the genes involved in cancer-related pathways, namely cell cycle regulation, WNT pathway, p53 pathway, and TGF-β pathway, were evaluated in their TSS200CGIs, except for CDKN2A. This lacked probes in its TSS200CGI, and its methylation was evaluated by a genomic block most 5′ of its exon 1, as described [36]. DNA methylation data was submitted to the Gene Expression Omnibus (GEO) database (accession no. GSE149255).

Analysis of cell cycle
Cultured cells were washed with 1× PBS (-) and fixed with 70% ethanol. The fixed cells were treated with RNase A (Thermo Fisher Scientific, Waltham, MA) for 20 min, stained with 50 μg/ml propidium iodide (Sigma-Aldrich) for 10 min, and analyzed using an Attune NxT Flow Cytometer (Thermo Fisher Scientific).