Epigenomics alternations and dynamic transcriptional changes in responses to 5-fluorouracil stimulation reveal mechanisms of acquired drug resistance of colorectal cancer cells.

A drug-induced resistant cancer cell is different from its parent cell in transcriptional response to drug treatment. The distinct transcriptional response pattern of a drug-induced resistant cancer cell to drug treatment might be introduced by acquired DNA methylation aberration in the cell exposing to sustained drug stimulation. In this study, we performed both transcriptional and DNA methylation profiles of the HCT-8 wild-type cells (HCT-8/WT) for human colorectal cancer (CRC) and the 5-fluorouracil (5-FU)-induced resistant cells (HCT-8/5-FU) after treatment with 5-FU for 0, 24 and 48 h. Integrated analysis of transcriptional and DNA methylation profiles showed that genes with promoter hypermethylation and concordant expression silencing in the HCT-8/5-FU cells are mainly involved in pathways of pyrimidine metabolism and drug metabolism-cytochrome P450. Transcriptional analysis confirmed that genes with transcriptional differences between a drug-induced resistant cell and its parent cell after drug treatment for a certain time, rather than their primary transcriptional differences, are more likely to be involved in drug resistance. Specifically, transcriptional differences between the drug-induced resistant cells and parental cells after drug treatment for 24 h were significantly consistent with the differentially expressed genes (termed as CRG5-FU) between the tissues of nonresponders and responders of CRCs to 5-FU-based therapy and the consistence increased after drug treatment for 48 h (binomial test, P-value=1.88E-06). This study reveals a major epigenetic mechanism inducing the HCT-8/WT cells to acquire resistance to 5-FU and suggests an appropriate time interval (24-48 h) of 5-FU exposure for identifying clinically relevant drug resistance signatures from drug-induced resistant cell models.

Introduction

Colon and rectal cancer (CRC) is a fatal cancer with low overall 5-year survival rate.[1] Until now, 5-fluorouracil (5-FU)-based chemotherapy is the first-line clinical treatment of CRC,[2] but the response rate is only ∼50%,[3] and most of patients acquire resistance during chemotherapy.[4, 5, 6] The process of a tumor acquiring drug resistance during chemotherapy could be mimicked partially by the common process of establishing drug-resistant cells by exposing parental cells to drug treatment in vitro over a period of time. Therefore, it should be helpful to explore drug-induced resistant cell models to understand both innate and acquired drug resistance mechanisms of cancers.

Based on a drug-induced resistant cell model, researchers often first identify differentially expressed genes (DEGs) between parental and resistant cells, termed basally deregulated (BD) genes,[7, 8, 9] and defined these DEGs as drug resistance genes. However, it has been recognized that most of these BD genes might simply represent drug-induced transcriptional changes irrelevant to drug resistance.[10, 11, 12] Recently, taking the DEGs between the nonresponders and responders of CRC patients receiving 5-FU-based therapy (termed CRG5-FU) as reference, we reported that inducible difference (ID) genes, which represent the transcriptional differences between resistant cells and parental cells synchronously exposing to drug for a certain time, are more likely to be involved in drug resistance. This result suggests a novel strategy to identify clinically relevant drug resistance genes from drug-induced resistant cell models.[10] However, the cells analyzed in our previous study were treated with 5-FU for only 6, 12 and 24 h. Thus, it is necessary to prolong the time of drug treatment to investigate the characteristics of sustained responses.[13, 14] Another problem remains to be addressed is to explain the underlying molecular mechanism for the acquired distinct transcriptional response characteristics of a drug-induced resistant cancer cell to drug treatment. It has been reported that aberrant DNA hypermethylation within gene promoters and consequent gene silencing might play a major role in generating drug-resistant phenotypes.[15, 16, 17] Therefore, we could assume that the distinct transcriptional response characteristics of a drug-induced resistant cancer cell to drug treatment could be introduced by acquired DNA methylation aberration in the cell exposing to sustained drug stimulation.

In this work, we performed both transcriptional and DNA methylational profiles of HCT-8 wild type cells (HCT-8/WT) for human CRC and its 5-FU-induced resistant cell line (HCT-8/5-FU). We first uncovered that the genes with both promoter hypermethylation and expression silencing during 5-FU treatment are mainly involved in pathways of pyrimidine metabolism and drug metabolism-cytochrome P450 that might be a major epigenetic mechanism of a cancer cell to acquire 5FU resistance. Then, we confirmed that ID genes detected at the 24 h time point were significantly consistent with CRG5-FU and the consistency became more significant after drug treatment for 48 h (binomial test, P-value=1.88E−06). This result suggested an appropriate time interval (24–48 h) of drug treatment for identifying genes related to drug resistance.

Materials and methods

Cell lines and reagents

HCT-8/WT and HCT-8/5-FU cells were purchased from Keygentec (Nanjing, China). Complete growth medium RPMI-1640 that contains 10% fetal calf serum and 1% penicillin–streptomycin solution was used for cell culture. Both cell lines were maintained at 37 °C in a humidified atmosphere containing 5% CO2. MTT assay was performed to determine the 5-FU resistance in both of HCT-8/WT and HCT-8/5-FU cell lines. 5-FU and penicillin–streptomycin solution (100 ×) were purchased from Sigma-Aldrich Corp. St. Louis, MO, USA. RPMI-1640 and fetal calf serum were purchased from Thermo Fisher Scientific, 168 Third Avenue, Waltham, MA, USA. Cell culture plates and other consumables were purchased from Corning and Axygen, 836 North Street, Building 300 Suite 3401, Tewksbury, MA 01876, USA.

Gene expression and methylation profiling

HCT-8/WT and HCT-8/5-FU cells were cultured in 5-FU free medium for 48 h, and then both cell lines were transferred to the medium for 0, 24 and 48 h with HCT-8/WT half-maximal inhibitory concentration doses of 5-FU (Figure 1). Untreated cells (0 h) and treated cells (24 and 48 h) were subject to extract RNA to be subjected to RNA microarray analysis. Total RNA was isolated from three independent experiments using TRizol (Thermofisher) according to the manufacturer’s instructions. Before the following experiments, total RNA was quality controlled by formaldehyde denatured agarose gel electrophoresis and ultraviolet–visible spectrophotometer. Total RNA (100 ng) from one experiment was sent to commercial kits for complementary DNA synthesis, complementary RNA synthesis, antisense RNA synthesis and its fragmentation (1792, Ambion, Thermo Fisher Scientific, 168 Third Avenue, Waltham, MA, USA). Fragmented antisense RNAs were hybridized to Affymetrix GeneChip PrimeView Human Gene Expression Array (3420 Central Expressway, Santa Clara, CA 95051) following the manufacturer’s protocols.

Figure 1

The experimental design in the study.

Methylation profiling of HCT-8/WT and HCT-8/5-FU cells with three independent experiments, respectively, were performed using Infinium HumanMethylation450 BeadChips (Illumina, 5200 Illumina Way, San Diego, CA 92122, USA). A total amount of 1 mg genomic DNA from one experiment was bisulfite modified using the EZ DNA Methylation Kit (Zymo Research, Irvine, CA, USA). In all, 200 ng of converted DNA was further processed to run BeadArrays according to the manufacturer’s instructions. Each locus is represented by fluorescent signals from two bead types corresponding to the methylated (M) and unmethylated (U) alleles, respectively.

Methylation and expression profiling data were available in NCBI (National Center for Biotechnology Information) Gene Expression Omnibus (GEO accession number: GSE81008).

Data processing

Gene expression profiling data were preprocessed with the Robust Microarray Average algorithm. Each probe set ID was mapped to its Entrez gene ID with the corresponding custom CDF files. If multiple probe sets were mapped to the same gene, the expression value for the gene was defined as the arithmetic mean of the values of the multiple probe sets (on the log2 scale).

The raw signals of unmethylated (U) and methylated (M) bead types were background corrected and computed into β-values using the GenomeStudio software Methylation Module (5200 Illumina Way, San Diego, CA 92122, USA). A total of 25 943 CpG loci located within the proximal promoter regions of the transcription start sites of 14 109 genes were used as background CpG loci in this study. The β-value was used to estimate the methylation level of the CpG locus using the ratio of intensities between methylated and unmethylated alleles:

Concordance score

For DEGs from two independent data sets sharing k DEGs, of which s genes have the same up- or down-regulation directions between the treated cell type and the control cell type, the concordance score was calculated as s/k. If k genes have both methylation and expression changes, among which s genes are hypermethylated (or hypomethylated) and correspondingly downregulated (or upregulated), then the concordance score was calculated as s/k.

The probability of observing a concordance score (s/k) by chance can be evaluated using the cumulative binomial distribution model as following:

in which p is the probability of one gene having the same dysregulation direction in two gene lists by random chance (here, p=0.5).

Identification of reproducible DEGs and differentially methylation genes in cell lines

The fold change and average difference metrics are often used to identify DEGs between two types of cell lines with two or three technical replicates.[10, 18] However, both of the methods depend on an arbitrarily determined cutoff value and lack any statistical control.[18] Here, we applied the reproducibility-based pairwise difference (PD) and pairwise fold change (PFC) methods to identify reproducible DEGs and differentially methylation genes between two types of cell lines based on small-scale cell line data.[18] Briefly, the algorithm treats every paired technical replicates of two cell types as an independent experiment and then identifies genes with top-ranked fold change or average difference between the two cell types that are significantly reproducible in the independent experiments. The DEGs and differentially methylation genes between the two cell types were finally identified according to the reproducibility of dysregulation directions in a treated cell type compared with the control cell type. The two parameters of the algorithm, the initial step (that is, the initial number of the top-ranked genes) and the consistency threshold were set as 300 and 90%, respectively, as suggested in our previous study.[18]

Human protein–protein interaction (PPI) data

The PPI data including 142 583 distinct interactions and 13 693 human proteins were collected as previously described.[19] The types of interactions between proteins included physical interaction, transcriptional regulation and sequential catalysis. Fisher’s exact test was used to test whether the direct PPI links between two gene sets were significant more than what is expected by random chance.

From the Drugbank,[20] KEGG (Kyoto Encyclopedia of Genes and Genomes)[21] and PharmGKB (The Pharmacogenomics Knowledgebase)[22] databases, we collected 291 genes associated with the metabolism, transport or activity of 5-FU, of which 254 could be mapped to the PPI network, denoted as 5-FU nodes.

Results

Acquired methylation aberrations in the drug-resistant cells induced by drug treatment

The half-maximal inhibitory concentration of the HCT-8/5-FU cells determined by MTT assay was 4746.31 μmol l−1, and this was ∼45-fold higher than the HCT-8/WT cells (105.06 μmol l−1). Thus, the resistance level of HCT-8/5-FU cells was large enough for analyzing molecular changes associated with the mechanism of drug resistance.[23] Then, we explored the potential molecular mechanisms determining the resistance of HCT-8/WT cells to 5-FU stimulation. Using the PD and PFC methods (see Materials and methods), we detected a total of 1658 genes with hypermethylated CpG sites located at the promoter regions in the HCT-8/5-FU compared with the HCT-8/WT cells. Among the 1658 hypermethylated genes, 459 genes were also identified as DEGs by the PFC and PD methods, and 67.76% (311) of these genes showed downregulation in the HCT-8/5-FU cells compared with the HCT-8/WT cells. This concordance score between hypermethylation and downregulation was unlikely to occur by chance (binomial test, P-value=1.07E–14; see Materials and Methods).

For the 311 hypermethylation-mediated downregulated genes in the HCT-8/5-FU cells compared with the HCT-8/WT cells, 242 genes were mapped in PPI network. Among the 242 genes, 24.79% (60) had direct PPI links with at least one of the 254 genes associated with the metabolism, transport or activity of 5-FU, documented in Drugbank, KEGG and PharmGKB databases (see Materials and Methods), significantly higher than the corresponding frequency of 17.91% (2416) for the rest 13 492 genes measured in both the gene expression and methylation platforms in the PPI network (Fisher’s exact test, P-value=8.65E–04). As shown in Figure 2, the hypermethylation-mediated downregulated genes participate in pyrimidine metabolism, drug metabolism-cytochrome P450, epidermal growth factor receptor signaling and p53 signaling pathways. The results implied that the hypermethylation-mediated downregulation of genes in the HCT-8/5-FU cells with prolonged time of drug treatment might lead to acquired resistance to 5-FU.

Figure 2

The hypermethylation-mediated downregulations of genes in the HCT-8/5-fluorouracil (5-FU) cells have frequent protein–protein interaction (PPI) links with 5-FU activity-related genes. 5-FU activity-related genes: genes involved in 5-FU transport, metabolism and other downstream effects (such as DNA repair, apoptosis and cell cycle regulation) in the public databases. The green nodes denote 5-FU activity-related related genes. The yellow nodes denote the hypermethylation-mediated downregulations of genes in the HCT-8/wild-type (WT) compared with the HCT-8/5-FU cells. The red nodes were overlapped between the two kinds of genes.

ID48 genes were most significantly consistent with CRG5-FU

The 131 CRG5-FU genes, which represent DEGs between CRC tissue samples of nonresponders and responders for 5-FU-based therapy extracted from three independent data sets,[10] were used to evaluate the clinical relevance of a list of candidate drug resistance genes at the transcriptional level.

First, we further confirmed the previous report[10] that BD genes, which represent DEGs between parental and resistant cells, mainly represent drug-induced changes. We combined BD genes detected by either the PFC or PD methods and obtained 6405 BD genes in total. The concordance score between BD genes and CRG5-FU was 66.67%, suggesting significant but still weak consistency (Table 1). In addition, the concordance score of BD genes with IP24 and IP48 were 97.84% and 98.85%, respectively (binomial test, all P-values <1.11E−16, Table 2). Here, the IP24 and IP48 genes represented DEGs between parental and drug-treated parental cells that underwent drug treatment for 24 and 48 h, respectively. The results confirmed that most of the BD genes detected in the HCT-8 cells mainly represent drug-induced changes.[10]

Table 1

The consistency scores between BD genes, ID genes and CRG5-FU

Gene set	Overlapped DEGa	Consistent DEGb	Consistency score (%)c	Binomial P-valued
BD	54	36	66.67	9.92E−03
ID24	21	15	71.43	3.92E−02
ID48	71	55	77.46	1.88E−06

Abbreviations: BD, basally deregulated; DEG, differentially expressed gene; 5-FU, 5-fluorouracil; ID, inducible difference.

CRG5-FU indicates DEGs between the pre-chemotherapy clinical specimens of responders and nonresponders receiving 5-FU-based chemotherapy; ID genes are genes detected between drug-treated parental cell line and drug-treated resistant cell line; BD genes are basally deregulated genes detected between parental cell line and resistant cell line.

The number of BD genes or ID genes overlapped with CRG5-FU.

The number of overlapped DEGs with the same deregulation direction.

The consistency score of BD genes or ID genes with CRG5-FU.

The binominal distribution P-value.

Table 2

The consistency scores between BD genes and IP genes

Gene set	Overlapped DEGa	Consistent DEG	Consistency score (%)b	Binomial P-value
IP24	1667	1631	97.84	<1.11E−16
IP48	2176	2151	98.85	<1.11E−16

Abbreviations: BD, basally deregulated gene; DEG, differentially expressed gene; IP, inducible parental gene.

The number of BD genes overlapped with IP genes.

The consistency score of BD genes and IP genes.

The PD and PFC methods were further used to identify ID genes that represent the transcriptional difference between resistant cells and parental cells in response to drug treatment for a certain time. The ID genes were denoted as ID24 and ID48 for the conditions where cells underwent drug treatment for 24 and 48 h, respectively. In total, 1922 ID24 genes and 6177 ID48 genes were identified. The ID24 genes overlapped with 21 genes of the CRG5-FU and the concordance scores of the overlapped genes was 71.43% (binomial test, P-value=3.92E−02, Table 1), whereas the ID48 genes overlapped with 71 genes of the CRG5-FU and the corresponding score increased to 77.47%. The results showed that the ID48 genes were most significantly consistent with CRG5-FU (binomial test, P-value=1.88E−06, Table 1). Especially, many of the 55 ID48 genes that were consistent with CRG5-FU (Supplementary Table S1) are known to be correlated with drug resistance or tumor cell proliferation, as exemplified in Supplementary Table S2 for the top 20 ID48 genes ranked by the PFC method. For example, upregulation of CD24 in CRC cell lines exhibits cancer stem cell-like characteristics and enhances 5-FU resistance.[24] SLC35A1 is a member of solute carriers and its overexpression can activate the process of absorption and transport of cell inhibitors.[25] DLGAP5, DTL, FGFBP1 and CDKN3 promote tumor cell proliferation and their overexpression could stimulate drug resistance.[26, 27, 28]

Taken together, the above results suggested an appropriate time point for identifying ID genes from the drug-induced resistant cell model.

Discussion

5-FU-based chemotherapy is widely used in the treatment of CRC and other solid tumors like gastrointestinal, head and neck and breast cancers.[29] However, response rates of patients with solid tumors to 5-FU-based chemotherapy are still very low.[30] Our analyses suggested that aberrant promoter hypermethylation might be a major epigenetic mechanism contributing to the acquired resistance of cancer cells to 5-FU, supporting the previous reports that CRC patients with a CpG island hypermethylation phenotype could not benefit from 5-FU-based adjuvant chemotherapy.[31, 32]

Furthermore, the PPI network analysis showed that many downregulated genes mediated by promoter hypermethylation had intensive direct interactions with genes participating in pyrimidine metabolism, drug metabolism-cytochrome P450, epidermal growth factor receptor signaling and p53 signaling pathways. Especially, three hub nodes (CYP3A5, GSTP1 and GPX2) in the PPI network are in the drug metabolism-cytochrome P450 pathway (Figure 2) and their suppression by promoter hypermethylation could induce drug inactivation within CRC cells[33, 34] and/or promote CRC cell survival and migration.[35] Notably, the suppression of some of these genes, such as GSTP1,[34] might mediate the multidrug resistance in CRC cells. More generally, some aberrant pathways induced by the epigenetically altered genes might increase multidrug resistance levels regardless of canonical drug mechanisms.[36] Further study is needed to investigate whether the acquired epigenetic aberrations in a resistant cancer cell induced by a particular drug could also increase the cell’s resistance to other drugs that might be an epigenetic mechanism of clinical multidrug resistance.

The DNA methylation aberrations acquired in the process of inducing drug-resistant cells make the drug-induced resistant cancer cells different from the parental cells in transcriptional response to drug treatment. Our analysis confirmed that genes with transcriptional differences between a drug-induced resistant cell and its parent cell after drug treatment for 24–48 h are more likely to be involved in drug resistance. Compared with the ID24 genes,[10] the ID48 genes became more consistent with the clinically relevant drug resistance genes (CRG5-FU). Obviously, if the drug treatment time is too short, the cancer cells could not receive the maximal benefit from drug dosing[23] and the transient changes in expression levels might be unstable.[10, 37] On the other hand, if the drug treatment time is too long, the expression changes of genes might reflect secondary and tertiary responses[37] and could not mirror the cycles of chemotherapy that patients receive in clinical treatment.[23] In fact, we have treated both HCT-8/WT and HCT-8/5-FU cells with 5-FU for 72 hours. The log2FC values and AD values of the 55 ID48 genes which were consistent with CRG5-FU were shown in Figure 3a–b. We found that the transcriptional changes of these genes at 72 h time point were different from those at 24–48 h. Moreover, we found the cellular morphology feature and density changed largely at 72 h, which would be unsuitable for experimental operations to obtain signatures related to mechanisms of drug action.[37]

Figure 3

The log2 fold change (FC) values (a) and average difference (AD) values (b) of 55 genes of 5-fluorouracil (5-FU) resistance.

In summary, this study reveals a major epigenetic mechanism of the HCT-8/WT cells to acquire drug resistance to 5-FU. Furthermore, the results suggest an appropriate time interval (24–48 h) to extract clinically relevant drug resistance signatures from drug-induced resistant cell models. Further study should be done to explore whether we could extract predictive signatures from the drug resistance genes identified from drug-induced resistant cell models to predict responses of patients to 5-FU-based chemotherapy.[10, 38, 39]