Predicting Carcinogenic Mechanisms of Non-Genotoxic Carcinogens via Combined Analysis of Global DNA Methylation and In Vitro Cell Transformation

An in vitro cell transformation assay (CTA) is useful for the detection of non-genotoxic carcinogens (NGTXCs); however, it does not provide information on their modes of action. In this study, to pursue a mechanism-based approach in the risk assessment of NGTXCs, we aimed to develop an integrated strategy comprising an in vitro Bhas 42 CTA and global DNA methylation analysis. For this purpose, 10 NGTXCs, which were also predicted to be negative through Derek/Sarah structure–activity relationship analysis, were first tested for transforming activity in Bhas 42 cells. Methylation profiles using reduced representation bisulfite sequencing were generated for seven NGTXCs that were positive in CTAs. In general, the differentially methylated regions (DMRs) within promoter regions showed slightly more bias toward hypermethylation than the DMRs across the whole genome. We also identified 13 genes associated with overlapping DMRs within the promoter regions in four NGTXCs, of which seven were hypermethylated and six were hypomethylated. Using ingenuity pathway analysis, the genes with DMRs at the CpG sites were found to be enriched in cancer-related categories, including “cell-to-cell signaling and interaction” as well as “cell death and survival”. Moreover, the networks related to “cell death and survival”, which were considered to be associated with carcinogenesis, were identified in six NGTXCs. These results suggest that epigenetic changes supporting cell transformation processes occur during non-genotoxic carcinogenesis. Taken together, our combined system can become an attractive component for an integrated approach for the testing and assessment of NGTXCs.


Transformation Assay
As a first step towards the development of an integrated approach, the in vitro Bhas 42 CTA suggested by the OECD Working Group [8] was performed on 10 NGTXCs. Cell growth assays were conducted to identify the most suitable concentration of each test article for the Bhas 42 CTA (data not shown). Three independent assays were carried out to assess cytotoxicity. Five to seven concentrations of each NGTXC were established based on the results of the cell growth assay. The highest concentration for the cell growth assay was set to 10 mM or 2 mg/mL according to an OECD guidance document [8].
In the case of insoluble substances, the concentration range for the growth assay was set to include one concentration with a visible precipitate due to the limit of solubility. However, in the transformation assay for melamine, rosuvastatin, and D-limonene, no foci were formed even in the concentration range where no cytotoxicity was observed in the cell growth assay (data not shown). Thus, for these three test substances, the CTA was performed at a re-adjusted concentration range, which is low enough to form foci. The promoter activity of 10 NGTXCs was examined, and the results are summarized in Table 3. Among the 10 NGTXCs, three (D-limonene, melamine, and rosuvastatin) showed negative results in the promotion test. Sodium saccharin, methapyrilene, DEHP, diethanolamine (DEA), TCDD, and cholic acid induced a more than two-fold increase in foci formation, while okadaic acid led to a relatively weak increase.

Identification and Distribution of DMRs
To detect the critical methylation events associated with NGTXCs, we compared the DNA methylomes from seven Bhas 42 CTA-positive NGTXCs via RRBS analysis, which is an effective and economical method for analyzing the DNA methylation of a single base [18,19]. DMRs, defined as genomic regions exhibiting methylation differences with more than 20% changes between the samples (untransformed cells vs. transformed foci) at a CpG site and achieving a false discovery rate (FDR)-adjusted p-value < 0.05, were used for all the subsequent analyses. The complete list of DMRs associated with each NGTXC is included in Supplementary Tables S1-S7. The average number of DMRs identified in each NGTXC are 730-1574, with the exception of 3796 DMRs for DEHP and 288 DMRs for okadaic acid (OA). DMRs identified in OA were relatively low compared to other CTA-positive NGTXCs probably due to the low transforming activity of OA as shown in Table 3. We calculated the coefficient of correlation to investigate whether there was any relation between the number of DMRs and in vitro cell transformation. The number of DMRs appeared to be positively related to the number of foci formed in CTA, but the coefficient of correlation was low with the value of 0.4782. The scatter diagram and coefficient of correlation are presented in Supplementary Figure  S1. Interestingly, of the 3796 DMRs in DEHP, 3081 DMRs (81%) were hypermethylated, while the remaining 715 (19%) were hypomethylated ( Figure 1A). Conversely, in DEA-induced transformed foci, 730 significant DMRs were detected, of which 146 (20%) were hypermethylated and 583 (80%) were hypomethylated. The DMRs in the promoter region, which were defined as DMRs within 1500 bp upstream and 1500 bp downstream from the transcription start site (TSS), were then segregated on the basis of methylation status ( Figure 1B). With the exception of DEHP, the DMRs in the promoter regions showed slightly more bias toward hypermethylation than DMRs throughout the whole genome. Figure 1C shows the heatmap of the DNA methylation level at DMRs specific for each NGTXC using hierarchical agglomerative clustering with complete linkage.
Then, the distributions of the DMRs from the RRBS data were compared in terms of the following seven genomic features: (i) promoters (DMR within 1500 bp upstream and 1500 bp downstream from TSS); (ii) exon; (iii) intron; (iv) 5 -UTR; (v) 3 -UTR; (vi) downstream (a 3-kb region downstream of transcription termination sites); and (vii) distal intergenic regions. The distribution patterns appeared to be relatively similar among the six NGTXCs, except for OA ( Figure 1D). Almost 40% of the identified DMRs were associated with the gene body (exon + intron). DMRs located within the intergenic regions were almost 40%. We also discovered that, on average, 15% of the identified DMRs were located in promoter regions, which are crucial for gene expression regulation. However, in the case of OA, DMRs were observed in the promoter regions at a frequency of 26.4%.

Overlapping DMRs in Promoter Regions across NGTXCs
We focused on the DMRs in promoter regions that could be directly associated with gene expression. The overlapping DMRs with a methylation difference >20% and adjusted p-value < 0.05 in the promoter region (1500 bp upstream and 1500 bp downstream from the putative TSS) were selected. Of these shared DMRs in NGTXCs, the DMRs exhibiting methylation changes in the same direction were selected as "overlapping" DMRs. A complete list of the overlapping DMRs in promoter regions between each NGTXC is available in Supplemental Table S8. In particular, there were 59 overlapping DMRs in the promoter region between DEHP and CA. This accounts for 31% of the 188 DMRs found in CA. However, the mode of actions (MOAs) of DEHP and CA were peroxisome proliferators and xenobiotic metabolizing enzyme inhibitors, respectively, and there was no similarity between the non-genotoxic carcinogenic mechanisms. Only three DMRs were found between OA and SS. Table 4 lists the genes harboring the overlapping DMRs in four or more NGTXCs. The 13 genes associated with overlapping DMRs in four NGTXCs were identified, of which seven were hypermethylated and six hypomethylated. Although two genes (H2-D1 and Mir705) harboring overlapping DMRs were found to be significantly differentially methylated between the transformed foci and controls across the five NGTXCs, the roles of these genes in cancer have not yet been investigated.

Enriched Functional Annotation and Canonical Pathway
The functions associated with genes harboring DMRs identified for each NGTXC were annotated (Table 5). For enriched function annotation, we used a list of DMRs with more than 30% changes between the samples (controls vs. the transformed foci) at a CpG site and achieving an FDR-adjusted p-value < 0.05. In six NGTXCs, except SA, the genes with DMRs at CpG sites were significantly enriched in cancer-related categories, including "cancer", "cell-to-cell signaling and interaction", and "cell death and survival", although these were not the most enriched categories. On the other hand, the genes harboring DMRs in SA were significantly related to connective tissue development and function (p = 1.93 × 10 −5 ). Additionally, as shown in Table 6, several canonical pathways revealed the significant enrichment of the activation pathway of the retinoid acid receptor (RAR) and retinoid X receptor (RXR) families (in six NGTXCs), gustation pathways (in six NGTXCs), as well as the salvage pathway of pyrimidine ribonucleotides (in five NGTXCs), although these pathways were not observed as top hits. Overall, the annotations represent the trends of the malignant cell transformation of cells with NGTXC treatment.

Network Identification
We identified networks with scores over 5 for each NGTXC. Figure 2 depicts the top two highest-scoring IPA-generated networks. One network (control vs. DEHP) with a score of 32, including genes for ERK1/2, NFκB (complex), NOS1, NOTCH4, PITX2, and PTPN11, was associated with cell-to-cell signaling and interaction, gene expression, nervous system development, and function. The other network (control vs. OA) with a score of 28 included genes for Ccl2, IFNG, IRS1, RUNX2, TNF, and TP53, was associated with cancer, cell death and survival, and organismal development. Networks related to "cell death and survival", which were considered to be linked to carcinogenesis, were identified in the other six NGTXCs, except for CA. Supplementary Table S9 features detailed information, including the molecules, score, and focus molecules, for each network.  function. The other network (control vs. OA) with a score of 28 included genes for Ccl2, IFNG, IRS1, RUNX2, TNF, and TP53, was associated with cancer, cell death and survival, and organismal development. Networks related to "cell death and survival", which were considered to be linked to carcinogenesis, were identified in the other six NGTXCs, except for CA. Supplementary Table S9 features detailed information, including the molecules, score, and focus molecules, for each network.

Discussion
Genotoxic carcinogens (GTXCs) can be detected with a battery of genotoxicity tests, such as the measurement of primary DNA damage. However, NGTXCs act through a large and diverse variety of different and specific mechanisms that do not involve direct interaction with DNA [1]. Thus, in contrast to GTXCs, a panel of assays addressing particular biological endpoints will be needed for the detection of NGTXCs [27]. There is cumulative evidence that the in vitro Bhas 42 CTA is suitable for detecting NGTXCs acting via non-genotoxic mechanisms (Validation Study Report [4,5]). However, in vitro CTA approaches do not provide information on the non-genotoxic mechanisms of carcinogenicity. Recently, it has been suggested that the development of an IATA could be the

Discussion
Genotoxic carcinogens (GTXCs) can be detected with a battery of genotoxicity tests, such as the measurement of primary DNA damage. However, NGTXCs act through a large and diverse variety of different and specific mechanisms that do not involve direct interaction with DNA [1]. Thus, in contrast to GTXCs, a panel of assays addressing particular biological endpoints will be needed for the detection of NGTXCs [27]. There is cumulative evidence that the in vitro Bhas 42 CTA is suitable for detecting NGTXCs acting via non-genotoxic mechanisms (Validation Study Report [4,5]). However, in vitro CTA approaches do not provide information on the non-genotoxic mechanisms of carcinogenicity. Recently, it has been suggested that the development of an IATA could be the mechanistic basis for the carcinogenicity of NGTXCs [27]. IATA can include a combination of methods that includes in vitro testing, QSAR, and -omics technologies.
Here, in order to identify the epigenetic modifications associated with NGTXC-induced carcinogenesis, we provided an integrated strategy consisting of (QSAR), in vitro Bhas 42 CTA, and DNA methylomes employing RRBS analysis. DNA methylation is one of the most important epigenetic modifications because aberrant methylation patterns are associated with the development diseases like cancer [45,46]. RRBS is an efficient and high-throughput technique for analyzing genome-wide methylation profiles at a single nucleotide level [25]. In addition, RRBS has been used to study several human diseases, including cancer [47], neurodegenerative diseases [48], aging [49], and immunological diseases [50]. Ten test articles selected as NGTXCs in this study were also predicted to be negative through Derek/Sarah analysis. Among the 10 NGTXCs, seven had positive results in the promotion testing of the in vitro Bhas 42 CTA. The average number of DMRs identified in the CTA-positive NGTXCs was 730-1574, with the exception of 3796 DMRs for DEHP and 288 DMRs for OA. We also discovered that, on average, 15% of the identified DMRs were located in promoter regions, which are crucial for gene expression regulation. With the exception of DEHP, the DMRs in the promoter regions showed slightly more bias toward hypermethylation than DMRs across the whole genome. Moreover, the 13 genes associated with overlapping DMRs in four NGTXCs were identified, of which seven were hypermethylated and six hypomethylated. In particular, although two genes (H2-D1 and Mir705) harboring overlapping DMRs were found to be significantly differentially methylated between transformed foci and controls across the five NGTXCs, and the roles of these genes in cancer have not yet been investigated.
To identify the major genetic events and processes occurring during malignant cell transformation by NGTXCs, we examined enriched functional annotations. The IPA analysis showed that, in six NGTXCs, except SA, the genes with DMRs at CpG sites were highly enriched in cancer-related categories, including "cancer", "cell-to-cell signaling and interaction", and "cell death and survival", although these were not the most enriched categories. Overall, the annotations represent the trends of the malignant cell transformation of the cells with NGTXC treatment. Interestingly, the potential contribution of the synaptic transmission and development throughout cell transformation by NGTXCs is evident from the enrichment of genes with annotations related to "cell-to-cell signaling and interaction". Although immunological synapses have been described previously in the context of the immune escape of glioma [51], their role in cell transformation has not been determined as of yet.
Several canonical pathways were established to have a significant enrichment of the activation pathway of the RAR and RXR families (in six NGTXCs), gustation pathways (in six NGTXCs), as well as the salvage pathway of pyrimidine ribonucleotides (in five NGTXCs), though these pathways were not observed as top hits. With this, the cell growth and differentiation of malignant cells by retinoids has been reported to be regulated by RARs and RXRs [52]. In addition, our results are in line with previous reports indicating that gustation pathways and the salvage pathway of pyrimidine ribonucleotides were identified via IPA in cancer cells [53,54]. IPA also determined the functional networks relevant to NGTXC-induced cell transformation. The most significant network was defined as "cell-to-cell signaling and interaction" (with a score of 32), where genes including ERK1/2, NFκB (complex), NOS1, NOTCH4, PITX2, and PTPN11 are clustered. The second top-ranked network was "cell death and survival" (with a score of 28).
Taken together, we posit that the genes with DMRs at CpG sites in the transformed foci are linked to malignant cell transformation that the cells undergo during NGTXC treatment. Thus, the framework in this study demonstrated the possibility of using a combination of several tools to take advantage of an integrated approach. However, the functional annotations and canonical pathways identified through our IPA analysis did not exhibit a robust correlation with the modes of action previously reported for each NGTXC. Therefore, further study will be required for the development of a more efficient integrated approach that includes a methylation profile analysis.   was included as a positive control for the promotion assay. Test articles were dissolved in sterile distilled water or DMSO. The final concentration of distilled water was 1% in the culture medium, and that of DMSO was 0.1%.

Cell Lines and Cell Culture
Bhas 42 cells, which were established from the BALB/c 3T3 cells through transfection with the v-Ha-ras gene [55], were obtained from the Japan Health Sciences Foundation Health Science Research Resources Bank (Osaka, Japan). As previously described [24], the cells were grown in MEM supplemented with 10% FBS (M10F). Only 60-70% of the confluent cells were sub-cultured. The cells were cultured in a humidified incubator under 5% CO 2 and 95% O 2 . Three days prior to treatment with test articles, Bhas 42 cells were maintained in DMEM/F12 supplemented with 5% FBS (DF5F).

Cell Growth Assay
The cell growth assays were performed prior to the transformation assays to determine the appropriate treatment concentrations of the test articles, as well as concurrently with every transformation assay to estimate the effect of the test articles on cell growth and survival. The cells were seeded, treated with a test article, and cultured until day 7. Three wells were prepared for each treatment group. On day 7, the cells were fixed with 10% formalin and stained with a 0.1% CV solution in 5% ethanol. CV was extracted from the stained cells with a solution containing 0.02 mol/L hydrochloric acid in 50% ethanol. The optical density of the extract was measured at a wavelength between 540 and 570 nm.

In Vitro Bhas 42 CTA for the Promotion Test
The in vitro transformation assays were conducted according to the procedures reported by Sakai et al. [7] with slight modifications. The Bhas 42 cells were seeded into each well of 6-well microplates at 2 mL volumes (7000 cells/well, day 0). Six wells were prepared for each treatment group. The cells were cultured for four days without medium exchange. On day 4, day 7, and day 11, the culture medium was replaced with fresh medium containing a test article or vehicle alone, and the treatment in the promotion phase was continued until day 14. The media were not changed during the last 7 days of culture. The cells were then fixed with methanol and stained with Giemsa solution.

Isolation of Transformed Foci
For the methylation profiling analysis, the foci that formed in the culture dish were cloned with cloning cylinders (6 mm internal diameter, 8 mm external diameter, 8 mm length). The foci were inspected, and the desired foci were circled. Growth medium was removed from the dish and cells washed with PBS. Then, a sterile cloning cylinder with one edge coated with silicone grease was placed around the selected focus, and 40 µL of trypsin was added to the cloning cylinder well. Cell lines were established from each focus, which were placed in 6-well plates containing M10F. After cellular growth for 10 to 15 days with frequent changes of medium, the genomic DNA from one established cell line was analyzed for methylation with RRBS.

Genomic DNA Extraction and Reduced-Representation Bisulfite Sequencing (RRBS)
For the whole-genome DNA, RRBS was performed by C&K Genomics, Inc. (Seoul, Korea) Briefly, genomic DNA was quantified in a Victor2 spectrophotometer (PerkinElmer, Boston, MA, USA) using the Quant-iT™ PiboGreen™ RNA Assay Kit (Invitrogen, Carlsbad, CA, USA) according to the manufacturer's protocol. To construct the Msp1 and Apek1 digested RRBS library, 500 ng of genomic DNA was incubated at 37 • C for 24 h. The digested products were purified with a MinElute PCR Purification Kit (Qiagen, Valencia, SC, USA). After purification, the digested products were blunt-ended, and then the dA was added, followed by methylated-adapter ligation. A range of 160-420 adapter-ligated fractions was excised from a 2% agarose gel. Bisulfite conversion was conducted using a ZYMO EZ DNA Methylation-Gold Kit™ (ZYMO research, Irvine, CA, USA) following the manufacturer's instructions. The final libraries were generated by PCR amplification using PfuTurbo Cx Hotstart DNA Polymerase (Agilent Technologies, Santa Clara, CA, USA). RRBS libraries were evaluated by an Agilent 2100 Bioanalyzer (Agilent Technologies, Palo Alto, CA, USA). Libraries were quantified by qPCR using the CFX96 Real-Time System (Bio-Rad, Hercules, CA, USA). After normalization, the sequencing of the prepared library was conducted on the Nextseq system (Illumina, San Diego, CA, USA) with 75-bp paired-end reads.

RRBS Data Analysis
RRBS data were preprocessed by the following steps. Briefly, low-quality adapter sequences in the raw reads were trimmed using TrimGalore (version 0.6.5) with Cutadapt (version 1.18) [56]. Quality checks of the clean data were conducted with FastQC (version 0.11.5), then high-quality reads were mapped to the reference genome (Mouse genome assembly GRCm38.p6) using Bismark (version 0.22.3) with bowtie2 (version 2.3.5) [57]. Data quality was shown in Supplementary Table S10. During the NGTXC-induced cell transformation, the DMRs were identified using DMRcaller (version 1.18.0), which uses the Bismark output data as input [58]. We considered 100-bp tilling bins for computing the differentially methylated cytosines, and performed a score test between the methylated and total reads in a bin for the samples (untransformed cells vs. transformed foci). We selected the bins where the FDR-adjusted p-value was less than 0.05 and the difference in level of methylation was at least 20% in the CG context, containing at least four cytosines (with at least four reads per cytosine). The identified DMRs were annotated using ChipSeeker [59] with gene annotation information from the R package TxDb.Mmusculus.UCSC.mm10.knownGene and org.Mm.eg.db [60]. Promoters were defined as the DMRs within 1500 bp upstream and 1500 bp downstream from the TSS.

Molecular Network, Pathway, and Functional Analysis
Enriched functional annotation, gene network identification, and canonical pathway analysis (generalized pathways that represent the common properties of a particular signaling module or pathway) were performed using Ingenuity Pathway Analysis software (Ingenuity Systems, Redwood City, USA). We utilized a list of DMRs with more than 30% changes between the samples (controls vs. transformed foci) at a CpG site and achieving an FDR-adjusted p-value < 0.05. The networks were identified using the scores calculated for each network according to the fit of the network to the set of the focus genes of the ingenuity analysis. The significance of the biofunctions and canonical pathways were tested by p-value. In addition, the canonical pathways were ranked by the ratio (number of genes from the input DMR data that maps to the pathway divided by the total number of molecules that exist in the canonical pathway).

Conflicts of Interest:
The authors declare no conflict of interest.