DNA methylation changes from primary cultures through senescence-bypass in Syrian hamster fetal cells initially exposed to benzo [ a

Current chemical testing strategies are limited in their ability to detect non-genotoxic carcinogens (NGTxC). Epigenetic anomalies develop during carcinogenesis regardless of whether the molecular initiating event is associated with genotoxic (GTxC) or NGTxC events; therefore, epigenetic markers may be harnessed to develop new approach methodologies that improve the detection of both types of carcinogens. This study used Syrian hamster fetal cells to establish the chronology of carcinogen-induced DNA methylation changes from primary cells until senescence-bypass as an essential carcinogenic step. Cells exposed to solvent control for 7 days were compared to naïve primary cultures, to cells exposed for 7 days to benzo[ a ]pyrene, and to cells at the subsequent transformation stages: normal colonies, morphologically transformed colonies, senescence, senescence-bypass, and sustained proliferation in vitro. DNA methylation changes identified by reduced representation bisulphite sequencing were minimal at day-7. Profound DNA methylation changes arose during cellular senescence and some of these early differentially methylated regions (DMRs) were preserved through the final sustained proliferation stage. A set of these DMRs (e

1 Matthew J Meier, Cummings-Lorbetskie, C., Rowan-Carroll, A., and Desaulniers, D. (2023).Dataset on DNA methylation and gene expression changes induced by conducted.To date, there are no validated assays recommended by regulatory agencies to detect NGTxC (Jacobs et al. 2020;Luijten et al., 2020).Induction of epigenetic alterations (e.g., global genome DNA hypomethylation, hypermethylation of tumour suppressor genes and hypomethylation of oncogenes (Allis and Jenuwein, 2016;Feinberg et al. 2016;Gama-Sosa et al. 1983)) is a key characteristic of carcinogens (Smith et al., 2020).Such alterations of the epigenetic system occur in cancers regardless of the initial carcinogenic events (whether GTx or NGTx).Therefore, measurements of epigenetic changes are being considered to improve testing strategies for the identification of chemical carcinogens, including NGTxC (Desaulniers et al. 2021;Jacobs et al. 2020).
The current in vitro work uses the SH fetal cell model and aims at identifying DNA methylation changes that could be used to predict chemical-induced carcinogenic events.It provides empirical data to test the hypotheses that DNA methylation changes can be detected early during cellular transformation to predict senescence-bypass as an essential forward carcinogenic step and potential AOP key event.SH fetal cells were exposed to dimethylsulfoxide (DMSO) for 7 days and used as control group.The chronology of DNA methylation changes was established by comparing the control group first to naïve primary cultures, then to cells exposed for 7 days to B[a]P, and to cells at the subsequent transformation stages: normal colonies (Nc), morphologically transformed colonies (MTc), senescence (SEN), senescence-bypass (SENbp), and sustained proliferation in vitro (Fig. 1).The endpoints that were monitored include differentially methylated regions (DMRs) detected by reduced representation bisulphite sequencing (RRBS), global genome changes in DNA methylation, and DNA methylation of a DNA interspersed repeated sequence (Yamada et al. 2006) with sequence homology to the murine and human long interspersed nuclear element-1 open reading frame 2 (L1-ORF2), which is hereafter referred to L1-like (L1l 5 ).In a final experiment, to ensure the reproducibility of the identified DMRs, a set of these DMRs in 6 genes was confirmed by pyrosequencing assays using a series of clones obtained from different laboratories and animals.Overall, important cancer-related genes were found with DMRs appearing prior to and during senescence and that persisted during sustained proliferation at 213 days of culture.Such DMRs might be useful markers to develop assays to predict senescence-bypass as an essential carcinogenic step and to improve chemical carcinogen testing strategies.

Animal procedures and stock of primary fetal cells
Animal procedures and housing conditions were approved by the Animal Care Committee of Health Canada and conformed to the Guidelines of the Canadian Council on Animal Care.Pregnant SH (Charles River Laboratories, St. Constant, Quebec, Canada) arrived in our laboratory at gestation-day-8 (GD8) and were allowed to acclimatise until GD13.5 when the fetuses were collected.The methodologies to process the fetuses to generate cell suspension from the trunk and to cryopreserve the cells were based on previous protocols (OECD, 2015;Pickles et al. 2016) and adapted to our lab as described in the supplementary Word document.The sex of the fetus can impact the DNA methylation profile.Given that the fetal sex could not be visually determined even under a magnifying lens or stereoscope, a quantitative real-time PCR (qPCR) method was developed for sex determination of the fetus and later from the clones emerging from cultures of pools of cells from different fetuses (companion paper (Meier et al. 2023)).

Chronology of DNA methylation changes
The cells were cultured in complete growth media (CGM) composed of Leboeuf's DMEM (Quality Biochemicals, Gaithersburg, MD, USA) at pH 6.7, including 20% FBS (ThermoFisher), 1% GlutaMAX™ Supplement (ThermoFisher), and penicillin-streptomycin (50 U/mL, Thermo-Fisher), and incubated at 37 • C with 10% CO 2 .In contrast to the classical assay methodology that requires growing the SH fetal cells on a feeder layer of irradiated cells (OECD, 2015), here the cells were seeded in 50% preconditioned media as previously described (Pant et al. 2010;Pant et al. 2012;Pickles et al. 2016).To prepare the preconditioned media, 2.5 × 10 6 frozen cells were seeded in CGM in a 6 cm dish (Nunc, ThermoFisher), and 24 h later the cells were passaged into 3 dishes (10 cm) with 12 mL CGM each.The media was harvested 72 h later, cleaned through a 0.45 µm filter (Millipore, Mississauga, ON, Canada) and stored at 4 • C (for a maximum of two weeks).
Fig. 1 describes the experimental protocol we used that led to the collection of 8 cell groups to establish the chronology of DNA methylation changes toward senescence-bypass and sustained proliferation as essential key carcinogenic events.Freshly dissected cells were seeded in a 6 cm dish in CGM, then 24 h later the cells were lifted using TrypLE™ (ThermoFisher) and harvested.Viable cells were counted (Scepter™ handheld cell counter, SigmaAldrich) and used to either prepare aliquots of frozen cells for molecular analyses as primary cells (group-1), or were seeded in 6 well plates (Nunc, ThermoFisher) at a clonal density of 45 cells per well in 2 mL of 50% preconditioned CGM.After 24 h of culture, 1 mL of CGM containing 3x concentrated test chemical and/or vehicle was added to reach a final concentration of 5 µg/mL (20 µM) B [a]P (a concentration repeatedly demonstrated to induce SH cell transformation (OECD, 2007;OECD, 2015;Vasseur and Lasne, 2012)) and/or 0.2% DMSO.After 7 days of exposure, DMSO and B[a]P treated cells (colonies picked across wells to create pools of 60 colonies per sample) were collected and frozen, representing group-2 and -3, respectively.Additionally, four Nc and four MTc were picked from B[a]P treated wells and reseeded separately in 100% CGM containing no chemical treatment.The Nc and MTc were visually identified based on the disorganized growth pattern and changes in nuclear/cytoplasmic ratio, according to previously described SH fetal cell transformation photo catalogues (Bohnenberger et al. 2012;Maire et al. 2012).To obtain sufficient cells for DNA extraction and RRBS analyses, MTc and Nc were expanded until they reach a total of 19-25 days of culture for an average of 20 more population doublings.Nc and MTc represent group-4 and -5, respectively.
Samples at SEN, SENbp, and sustained proliferation were created as follow.From the remaining colonies that were exposed to B[a]P for 7 Fig. 1.Experimental and cellular events relative to days in culture and cumulative population doubling levels (PDL).The PDL was calculated as described in the supplementary document.The blue insert represents the lower graph.The cells were collected from fetuses at Day-0 and allowed to proliferate in culture for 24-48 h and then frozen in aliquots of 2.5 × 10 6 to establish a bank of primary cells.To perform an experiment the frozen primary cells were cultured for 24 h, then collected and counted.Viable cells were distributed into treatment groups and exposed to DMSO or B[a] P for 7 days.At that time (Day 9), separate cultures of Nc and MTc were initiated and terminated between Day 19-25 when a sufficient number of cells were collected for molecular analyses.Cultures eventually senesced, however one male clone HC26d1 emerged (upper panel).The dots represent passage time at which cells were counted and reseeded, while the arrow indicates the samples used for RRBS analyses representing senescence (growth arrest between Day 37-58), senescence-bypass (Day 94), and sustained proliferation (Day 213).
days, approximately one-half of each colony was collected by scraping with a sterile plastic pipet tip then cultured in 24-well dishes.The remaining half colonies were fixed and stained with premade Giemsa stain solution (Fisher Scientific, Ottawa, ON) to confirm the morphology as Nc or MTc.Colony-derived cultures were expanded to successively larger culture dishes until senescence.Out of 20 colonies that reached senescence, one (which we named clone HC26d1), showed growth arrest following 30 population doublings for a period of over 20 days prior to bypassing senescence and sustained growth until we stopped the culture.Four aliquots of the HC26d1 clone were collected during the growth arrest/senescence at 48 days of culture (passage 10, group-6), four other aliquots after senescence-bypass (Day 94, passage 15, group-7), and finally samples were collected at 213 days of culture (passage 31, group-8) after displaying sustained proliferation.
A sheet in the Excel supplementary file provides the sample details including fetal identification, sex, and dams of fetuses.The experiment used cells from male fetuses, except for one group of four sibling female fetuses representing the primary cell group.In summary, the chronology of DNA methylation changes (Fig. 1) was investigated using the following pairwise comparisons: 1. Dispersed primary cells from four female fetuses (expanded in culture for 24-48 h after dissection, no colony was formed at this early stage) vs a male fetus (id#10) primary cells exposed for 7 days to DMSO (0.2%) from which a pool of 60 colonies were collected and divided in 4 aliquots as the control group (CG).2. The CG vs a male (id#10) primary cells exposed for 7 days to B[a]P (20 µM) from which a pool of 60 colonies were collected and divided in 4 aliquots.3. The CG vs Nc (from id#10) picked after exposure to B[a]P and expanded in chemical-free medium until 19-25 days of culture (4 samples, 1 colony per sample).4. The CG vs MTc (from id#10) picked after exposure to B[a]P and expanded in chemical-free medium until 19-25 days of culture (4 samples, 1 colony per sample).5.The CG vs senescent cells (SEN, four aliquots) at 48 days of culture (originate from a pool of primary cells exposed to B[a]P from 5 male and 3 female fetuses identified as pool#1).6.The CG vs the male clone HC26d1 that emerged from pool#1 and bypassed senescence (SENbp).A sample was collected at day-94 of culture and divided in four aliquots.7. The CG vs the clone HC26d1 during sustained proliferation collected at day-213 of culture (four aliquots of one sample).8. Nc vs MTc (as previously described).9. SEN vs SENbp (as previously described).

Reduced representation bisulphite sequencing (RRBS) and bioinformatics
All RRBS analyses were from cells harvested using TryplE (Ther-moFisher), washed twice with dPBS, and frozen as cell pellet.Genomic DNA was extracted with DNEasy kits (Qiagen, Mississauga, ON, Canada).DNA libraries were prepared from 100 ng genomic DNA using the Premium RRBS kit from Diagenode (Denville, NJ, USA).Briefly, MspI digested DNA samples were individually barcoded, pooled, then bisulfite treated to convert unmethylated cytosine bases to uracil that are then amplified to thymine.Libraries were sequenced on an Illumina NextSeq 550.Fastq files were trimmed using TrimGalore-0.5.0 with parameters -Illumina -RRBS.Base calls were generated by Bcl2fastq v. 2.20.0.42.Alignments were generated with Bismark v0.22.1 (Krueger and Andrews, 2011) and Bowtie2 v2-2.3.5.1.The software MethylKit (v1.10.0, running in R v 3.6.1)(Akalin et al. 2012) was used to generate the count files through pairwise comparisons.The statistical test used by MethylKit for pairwise comparisons activated by the DiffMeth() function is a chi-squared test with overdispersion correction.DMRs were clusters of 200 base pairs (bp) with read depth > 20, q< 0.05, and methylation difference > |25%|.A cluster length of 200 bp is considered appropriate for CpG island and promoter analyses (Carmona et al. 2017), and the transition from an inactive to active transcription start site (TSS) can involve a chromatin segment of only a dinucleosomal length (Wolff et al. 2010) and a nucleosome is surrounded by a DNA segment of 114 bp.The RRBS data were deposited in the NCBI GEO super series GSE220238.

Ingenuity Pathway Analyses (IPA)
The IPA knowledgebase software (Qiagen) does not recognize the SH genome; consequently, the mouse orthologue genes downloaded from Ensembl Biomart were used instead.Percent methylation changes for DMRs located in gene promoters ( ± 3000 bp from the TSS) were imported into IPA to identify diseases and biofunctions (cut-off at Benjamini-Hochberg adjusted (B-H) p < 10 − 8.5 ) and canonical pathways (cut-off at B-H adjusted p < 0.01) enriched with affected genes.

Pyrosequencing
Bisulfite pyrosequencing assays were designed using Assay Design Software v2 (Qiagen) to measure methylation of individual CpGs.Post design, the PCR primer sequences were analysed for unintended alternative PCR products and mispriming sites using the online ePCR tool, Bisearch (bisearch.enzim.hu,(Aranyi and Tusnady, 2007)), with the bisulfite-converted SH genome selected.Oligonucleotide primers were purified by desalting except for the 5 ′ -biotin-labelled PCR primers, which were HPLC-purified (Integrated DNA technologies, Coralville, Iowa, USA).
Genomic DNA was isolated using Qiagen's mini DNEasy blood and tissue kit, and then 400 ng was bisulfite-converted using the EZ DNA Methylation-Gold kit (Zymo Research, Orange, CA, USA).The PCRs were carried out using a C1000 thermocycler (BioRad) in 1X Pyromark PCR mix (Qiagen) with 0.4 µM PCR primers in a total volume of 25 µL and amplified according to Qiagen's instructions.The amplification protocol started with polymerase activation at 95 • C for 15 min, then 45 cycles of denaturation (94 • C, 30 s), annealing (56 • C, 30 s) and extension (72 • C, 30 s), followed by a final extension at 72 • C for 10 min.CpG methylation was measured using a PyroMark Q96 MD instrument (Qiagen) operated by PyroMark CpG software (version 1.011, Qiagen).The specificity of the PCR was verified by a single band resulting from electrophoresis of amplicons through 2% agarose (Fig. S1).
To validate the pyrosequencing assays, highly methylated and unmethylated SH DNA was generated by enzymatic treatment of genomic DNA isolated from excess fetal tissue (dissected heads, see section 2.1.1).Highly methylated standard DNA was created by treating 1 µg of DNA with 25 units Sss1 methylase and 0.16 mM of S-adenosylmethionine in 1X buffer-2 (New England Biolabs, Ipswich, MA, USA) for 16 h at 37 • C followed by a 20 min enzyme denaturation step at 65 • C. Low methylated standard DNA was prepared by whole genome amplification of 25 ng DNA using the Repli-G kit (Qiagen) following the supplied protocol.Both DNA standards were cleaned on a DNEasy purification column (Qiagen).To further validate the pyrosequencing assays and test for PCR bias mixtures of each methylation standard were assayed at 0%, 25%, 50%, 75%, and 100% of the high methylated standard.
Pyrosequencing assays were designed for specific DMRs (Bhlhe22, Aifm3, Pou4f1, Bgalnt2, Gja8, Klf17) and against a repetitive element sequence cloned from SH genomic DNA (GenBank accession # AB185085, (Yamada et al. 2006)).The repetitive element assay targets methylation of cytosines at positions 354 and 430 within the AT rich constitutive heterochromatin sequence AB185085, which is widely distributed across the SH genome and similar to the rodent Long Interspersed Nuclear Element-1 (Line-1 or L1) (Yamada et al. 2006).The DNA sequence of AB185085 was used to search the database of repetitive DNA families (https://dfam.org.(Storer et al. 2021)), which yielded the highest similarity with "L1_Mur3_orf2: L1 Non-LTR Retrotransposon from Muridae" and "ORF2 from L1 retrotransposon, L1M2_orf2 subfamily" from Homo sapiens.Therefore, this assay is referred to here as Long interspersed nuclear element-1-like (L1l) assay.The Table S1 lists the pyrosequencing assay characteristics including the primer sequences and genome coordinates targeted for the DMRs in L1l and in 6 other genes Bhlhe22, Aifm3, Pou4f1, Bgalnt2, Gja8, Klf17, selected for validation across experiments (see next section).

Additional validation experiments
To gain confidence in the reproducibility of the DMR data across studies and laboratories, selected DMRs were measured at Health Canada (HC) in B[a]P treated SH fetal cell samples obtained from Brunel University (Pickles et al. 2016) and originating from two laboratories.These analyses were compared to cultures obtained at HC.The selected DMR-bearing genes were Bhlhe22, Aifm3, Pou4f1, Bgalnt2, Gja8, and Klf17, that were chosen based on their different methylation abundance, their genomic sequence permitting the design of high quality pyrosequencing assays, and their biological functions (supplementary Word document).To highlight concordance between pyrosequencing and RRBS data, correlation analyses between both data types are presented from an extended series of DMR-bearing genes (supplementary Word document, Fig. S4).Finally, to identify gene expression changes associated with DMRs, a companion paper (Meier et al. 2023) summarize correlations between DMRs identified by RRBS with the related RNAseq data from cells exposed to the DNA methylation inhibitor 5-aza-2 ′ deoxycytidine (5aCdR).

General observations
The RRBS analyses indicate relatively stable global genome DNA methylation at approximately 39% prior to SEN, then rises to 43% during SEN and SENbp, and descends to 41.5% at 213 days in culture (Fig. 2).These data suggest greater DNA methylation reprograming beginning after 19-25 days of culture, which is supported by hierarchical clustering of the RRBS data indicating separate groupings starting at senescence (Fig. S2).
Table 1 summarizes the number of DMRs based on cell transformation stages relative to the DMSO control group at 9 days of culture.For example, this table shows 198 DMRs (q < 0.05) between the vehicle and the B[a]P treated samples (7 days of B[a]P exposure with a total of 9 days of culture since fetal dissection).Among these 198 DMRs, 39 or 20% were located in promoter regions, while the other DMRs were in introns, exons, but mostly in intergenic sequences (Table S2).An important increase in the number of DMRs was occurring as the cells reached senescence, it increased from 337 to 1188 in the MTc and SEN groups, respectively.The Excel supplementary file provides the data for all contrasts described in Table 1, including primary cells vs the DMSO control group.The complete list includes 9698 DMRs with p < 0.05, or 5544 DMRs with q< 0.05 over nine comparisons/contrasts.

Chronology of gene specific methylation changes
IPA were performed to reveal potential links between DMR-bearing genes and cancer.Based on DMRs located in gene promoters ( ± 3000 bp from the TSS), IPA identified 29 diseases and biofunctions that were significantly enriched with DMRs across the transformation stages (Fig. 3).A cut-off of B-H adjusted p < 10 − 8.5 was used to ensure readable text in Fig. 3.The number of statistically significant diseases and biofunctions increased mostly at senescence.Of the 29 diseases and biofunctions, 26 were cancer related.The number of enriched canonical pathways (cut-off at B-H adjusted p < 0.01) increased with the Nc, but included pathways that differed from those in the MTc.Such difference is due to different DMR-bearing genes in both groups, as illustrated through the Venn's diagram analyses described in the next paragraph.
Venn's diagram analyses (Oliveros, 2007) were used to further explore differences in DMRs between Nc, MTc, SEN and SENbp (Figs. and 6).To create the Venn's diagrams the data set included the genome wide DMRs (i.e., not limited to the promoters) from the supplementary Excel File, filtered with q < 0.05 and methylation difference >|25%|.Entries were limited to those associated with a gene name.Only one entry was kept when the same gene included multiple DMRs so that each gene is represented only once in the Venn's diagrams.The filtered list of genes for each contrast (Fig. 4 The data from the contrasts (Ctrl vs B[a]P, Ctrl vs Nc, Ctrl vs MTc, Ctrl vs SEN, Ctrl vs SENbp, Ctrl vs Sustained proliferation) were aligned based on the common DMRs to identify those that persistent across the various events.Table 2 provides typical examples of persistent and transient DMRs with the percentage of DNA methylation change across time relative to the DMSO control group at 9 doc.For example, DMRs in Gja8 and Bhlhe22 were persistent, detected as early as in the contrast with MTc and were heritable across cellular generations until at least 213 doc, sometimes with progressive changes in methylation abundance (e.g.: Bhlhe22).The progressive changes in methylation were not only in magnitude, but also in length extending over 400 bp (e.g.: Pitx2,

Pou4f1
). Occasionally, DMRs were detected early but were missed in later contrasts (e.g.: Gja8, Phactr3), others like Klf17, were detected in only one contrast, whether these are truly transient methylation changes requires further investigation.Nevertheless, the series of contrasts permitted the identification of a list of 24 DMR-bearing genes (like Gja8 and Bhlhe22) that persisted across 4 or 5 contrasts with at least one of these contrast with > |25%| difference in methylation (q<0.05),these are: Adgrd1, Arhgef12, Arhgef3,Bhlhe22,Cacng4,Cfap126,Cul7,Dhx57,Eid3,Elavl2,Fam25a,Fxyd5,Gata6,Gja8,Gne,Krr1,Myog,Nppc,Palm2,Phc2,Slc36a4,Synm,Trpc6,Zfp64.IPA of this small list of 24 DMRbearing genes revealed their potential involvement in 253 categories of diseases and functions (B-H adjusted p < 0.05), and with cancer as a recurrent dominant descriptor of these categories (Table 3).The DMRbearing genes involved in the identification of diseases and functions are presented in two sheets of the supplementary Excel file (one sheet for all diseases, the other for cancers only).Furthermore, with these 24 DMR-bearing genes IPA identified connections with molecular networks that include major cancer related genes.A first molecular network of 35 molecules with the IPA title "Cell Cycle, Cellular Development, Organ Development", included 11 of the 24 DMR-bearing genes (IPA score 25, Fig. 6A).A second molecular network of 35 molecules entitled "Cellular Development, Cellular Growth and Proliferation, Hematological System Development and Function", included 10 of the 24 DMR-bearing genes (IPA score 22, Fig. 6B).These IPA Figures suggest that many of the 24 DMR-bearing genes can interact directly and/or indirectly with highly important cancer related genes in the first (e.g., MYCN, TP53, EZH2, RASSF, CDKN2A, P38MAPK, in Fig. 6A) and second network (e.g., CREBBP, NFkB, Akt, ERK1/2, in Fig. 6B).

Reproducibility of DMRs across clones and laboratories
The DMR-bearing genes in bold characters in Table 2 (Bhlhe22, Aifm3, Pou4f1, Bgalnt2, Gja8, and Klf17) were used to examine by pyrosequencing the reproducibility of DMRs across clones obtained in different laboratories and animals.The data are presented in Fig. 7.The regression lines show DNA methylation changes across the Health Canada (HC) samples in blue, including clone HC26d1 (data points represented by the acronym HC).The HC cultures started to senesce after approximately 32 population doublings or 35 days in culture (senescence/proliferation arrested from passage 8-10).Brunel University (Pickles et al. 2016) provided samples from two laboratories, these include 12 clones from B[a]P exposed SH fetal cells shown in red and numbered 1-12.A red line was drawn to connect the data when two samples from the same clone were collected at different passage number.A second set include DMSO treated samples (black open circles) from which, unexpectedly, a MTc (shown as black filled circles) developed but did not reach the stage of sustained proliferation.Generally, after passage 5 these graphs demonstrate reproducibility of trends across passage numbers and clones.Promoter DMRs in Bhlhe22, Aifm3, Pou4f1, and Bgalnt2, became hypermethylated through time.For example, a plateau was rapidly reached at senescence for Bhlhe22 but the increase in methylation was delayed after passage 15 for Pou4f1.Gja8 and Klf17 became hypomethylated through time, with Klf17 showing a longer delay prior to hypomethylation in the HC samples.Interestingly, Table 2 indicates that Klf17 hypomethylation was detected by RRBS in the contrast Ctrl vs SENbp but not in the later contrast Ctrl vs sustained proliferation; this observation is supported by some pyrosequencing data showing lower methylation abundance from passage 14-17 relative to later passages in the HC samples and DMSO-MTc samples (Fig. 7), both data sets (Fig. 7 and Table 2) indicate partial remethylation in late passages.Overall, Fig. 7 generally demonstrate the reproducibility of DNA methylation data across clones from different laboratories and animals.
The L1l data in Fig. 8, as well as from other HC cultures indicate a gradual decline in DNA methylation (from 33% down to 14%) from as early as 8 days in culture (Fig. S3), which is long before senescence and emergence of clones.Collectively with the other six genes (Fig. 7), the L1l assay pyrosequencing data reinforce the reproducibility of DNA methylation data across clones from different laboratories and animals.

Discussion
Our work established the chronology of DNA methylation changes from the primary cell stage through senescence and sustained proliferation.It identified persistent DMRs that could be tested for predictability of chemical-induced carcinogenic transformation steps.Here, exposure of SH fetal cells to B[a]P for 7 days induced morphologically transformed colonies from which a clone emerged and sustained proliferation, as previously showed by others (Pickles et al. 2016;Vasseur and Lasne, 2012).It has been reported that the injection of such morphologically transformed cells in syngeneic animals creates tumours frequently identified as sarcomas (Barrett et al. 1979;Elias et al. 1989;Pienta et al. 1977).We found that the cells at the end of the 7-day exposure to B[a]P had mostly transient DMRs no longer detectable at 19-25 days of culture in Nc or MTc, however the MTc colonies could be distinguished from the Nc by a series of DMRs.The MTc showed persistent hypo and hypermethylated DMRs that were preserved through senescence, senescence-bypass, and sustained proliferation at 213 days of culture.A list of 24 genes was established with DMRs detected mostly in MTc and that persisted across cellular generations, through cell transformation, SEN, and SENbp.IPA identified these genes in molecular network with well known cancer related genes, but also

Table 1
Number of differentially methylated regions (DMRs: clusters of 200 bp, read depth >20) based on cell transformation stages across the global genome (Total) or located in promoter areas ( ± 3000 bp from the transcription start site (TSS)) based on p value or q value adjusted for the false discovery rate.The first 6 groups were compared to the control DMSO-treated group, while the last 2 groups have their comparator listed below.with EZH2.The epigenetic enzyme EZH2 has the active site responsible for histone 3 lysine 27 methylation (H3K27), which play major roles in mediating the poised and silenced states of genes (Chammas et al. 2020;Schuettengruber et al. 2017).Deregulation of a network that involve EZH2 may promote further epigenetic changes and carcinogenesis (Chammas et al. 2020).The 24 DMR-bearing genes are important and may contribute to further carcinogenic steps or perhaps were the targets of previous cellular dysfunctions.Some of these DMRs (Bhlhe22,Gja8,Aifm3,Pou4f1,Bgalnt2,Klf17,L1l) were further investigated by pyrosequencing and were showed to be reproducible across clones and animals from other laboratories.Such reproducible and persistent DMRs deserve further investigation as potential early markers of SEN, SENbp and carcinogenicity.

B[a]P as carcinogenic and epigenotoxic agent
The International Agency for Research on Cancer classified B[a]P as a human carcinogen (IARC Working Group, 2012).Toxic and epigenetic effects of B[a]P exposures were recently reviewed (Bukowska et al. 2022;Bukowska and Sicinska, 2021).B[a]P binds and activates the aryl hydrocarbon receptor that increases the expression of cytochromes p450 involved in the transformation of B[a]P into mutagenic B[a]P diol-epoxide (BPDE) metabolites forming depurinating DNA adducts (Bausinger et al. 2016;IARC Working Group, 2012).BPDE preferentially creates DNA adducts on the guanine in methylated 5 ′ -CpG-3 ′ (CpG) sequences (Weisenberger and Romano, 1999).The methyl group creates hydrophobicity that stabilizes the BPDE-DNA adducts and thus CpGs are hot spot for mutations, for example on p53 (Zhang et al. 2005) and codon-14 on the K-ras oncogene (Hu et al. 2003).Similarly, molecular alterations in B[a]P-treated Syrian hamster cell cultures include transversion point mutations in the DNA-binding domain of p53 coupled with Ink4 alterations, loss of expression of p15 (Yasaei et al. 2013), Bmi1 upregulation, monoallelic deletion of the Cdkn2A/B locus, and p16 silencing through promoter methylation (Pickles et al. 2016).Here, the p16 promoter methylation was not affected.
The SH fetal cells were exposed to B[a]P for seven days (as recommended by (OECD, 2015)), which produced DMRs that were transient and others that persisted through time.Similarly, B[a]P induced transient and persistent hyper and hypomethylation events over 96 h in breast cancer cell lines (Sadikovic and Rodenhiser, 2006).The occurrence of transient and persistent DMRs may be linked to the sequence of toxic and transformation events following B[a]P exposure.The following kinetic analyses suggest mutagenesis as an early event.The maximum amount of DNA adducts occurred 1 h after the addition of BPDE to cultures of normal human skin fibroblasts, followed by a decline in abundance of BPDE-DNA adducts associated with excision repair and BPDE aqueous instability (Kootstra, 1982).This indicates that the B[a] P-mutagenic impact on gene expression and epigenetic patterns can occur early in the 7-day exposure period.Once bound to DNA, BPDE-DNA adducts change the chromatin configuration, create new transcription factor binding sites and inhibit others (MacLeod, 1996).BPDE-DNA adducts were shown to inhibit cytosine methylation at CpG sites (Subach et al. 2007;Wojciechowski and Meehan, 1984), even though they strengthen bonding of DNMT3a to the DNA, they reduce the methylase activity (Lukashevich et al. 2011).In contrast, the oxidative stress adduct 7,8-dihydro-8-oxoguanine (8-oxo-G) weakened DNMT3a bonding to the DNA, to also reduce DNA methylation (Lukashevich et al. 2011;Maltseva et al. 2009).Reactive oxygen species (ROS) created during transformation of B[a]P are important contributors to the oxidative stress, which is intensified by B[a]P-induced decreases in abundance of antioxidant molecules and enzymes (Bukowska et al. 2022).B[a]P-induced epigenotoxic effects are delayed relative to mitochondrial dysfunction in human peripheral lymphocyte cultures (Bhargava et al. 2020).In human bronchial epithelial cell cultures, B[a] P increases DNMT1 abundance (2.5 µM), but decreases the abundance of DNA methylation and DNMT3b methylase at 20 µM and 40 µM, respectively (Xia et al. 2016).Overall, dynamic changes in methylation may reflect sequences of genetic, toxic, and epigenetic events through time toward a resetting of DNA methylation and gene expression pattern conducive to carcinogenesis.The sequence of events between genotoxic and non-genotoxic carcinogens differ, but ultimately senescence bypass and sustained proliferation must be achieved as essential carcinogenic steps, which may be predictable with identified persistent DMRs.

Senescence
The cells in our study went through a period of growth arrest at passages 8-10 after approximately 30 population doublings, suggesting a senescent phase (Fig. 1, 37-58 days in culture).This observation is in line with reports indicating that cells from small rodents have limited life-spans and enter senescence under normal conditions after 20-30 population doublings (Russo et al. 1998).Senescence can be triggered by telomere shortening, activation of oncogenes, DNA damage,  mitochondrial dysfunction (Crouch et al. 2022;Wiley and Campisi, 2021), and/or activation of retrotransposons (De Cecco et al. 2019).Large increases in the total number of DMRs and many in promoters were observed in cells reaching senescence (Table 1).Widespread epigenetic changes are known to be associated with the induction of senescence, with the activation of the senescence-associated secretory phenotype, and the formation of senescence-associated heterochromatin foci (Crouch et al. 2022;Cruickshanks et al. 2013).Such senescence-associated epigenetic changes complicate the discovery of epigenetic markers responsible for SENbp and further carcinogenic steps.

Gene specific DMRs
Numerous factors can induce DMRs; indeed, replication stress (Nikolov and Taddei, 2016;Raurell-Vila et al. 2017), hypoxic stress (Lu et al. 2014), and cell culture itself is known to induce differential methylation in genes independent of the tumorigenic state (De Carvalho et al. 2012;Nestor et al. 2015).However, the reproducibility of identified DMRs across multiple clones obtained from different laboratories was demonstrated by pyrosequencing, highlighting hypermethylation in Bhlhe22, Aifm3, Pou4f1, Bgalnt2, and hypomethylation in Gja8 and Klf17 (the functions of these six genes are discussed in the supplementary Word file).It would be speculative to suggest that these DMRs are associated with increases or decreases in gene expression, or that the induced DMRs would be the same regardless of the initiating carcinogenic mechanisms.Indeed, the origin of oncogenic pathways (CYCLIN E1, WNT1, or HRasv12) influence tumour characteristics and their molecular subtype (Fonti et al. 2019), consequently, it can be speculated that other initiating mechanisms than B[a]P exposure may induce cellular transformation with different molecular subtypes.Nevertheless, these DMRs that are reproducible across clones, are among the list of 24 DMR-bearing genes detected early and that persisted across cellular generations through SEN and SENbp, all together these genes deserve further investigation as potential early epigenetic biomarkers of carcinogenesis.
The results suggest that Klf17 displayed DNA hypomethylation and partial remethylation after SENbp.Klf family members regulate numerous carcinogenic functions (Lewis et al. 2022;Mas et al. 2022;Zhu et al. 2022).Induced KLF17 expression suppresses growth of colorectal cancer cells (Jiang et al. 2019); hypothetically, perhaps a reduction in Klf17 DNA methylation around the time of SEN and early SENbp was required to increase Klf17 expression and to suppress SH cell proliferation during SEN.Then a remethylation of Klf17 would reduce its expression and favor cell proliferation after SENbp.Remethylation events are not unusual, for example, there are transient changes in the opposing DNA methylation and H3K4 methylation marks in regulating PD-1 expression during infections (Bally et al. 2020).Klf17 DNA methylation dynamic and roles during SEN and SENbp deserve further investigation.

Demethylation of L1l vs global genome DNA methylation
The L1l data revealed a gradual decline in methylation after 8 days of culture, which occurred whether the cells were treated or not with B[a] P. DNA methylation patterns as well as histone post-translational modifications, are known to change rapidly under in vitro conditions (Nestor DMR: genomic position of a 200 base pair region differentially methylated located at a distance less than 3000 bp from the transcription start site (TSS).This table provides examples of the same DMR within a gene detected in multiple contrasts with at least one with ≥ |25%| differences in methylation relative to the DMSO controls (q < 0.05).DMRs highlighted in bold were selected for further pyrosequencing validation (Fig. 7).B[a]P: colonies collected at 9 days of culture including a 7day benzo[a]pyrene exposure period.Nc: normal colonies.MTc: morphologically transformed colonies.SEN: senescence.SENbp: senescence-bypass.genome (Beck et al. 2022).More than 90% of all methylated cytosines reside within repetitive elements of constitutive heterochromatin despite their low CpG density (Beisel and Paro, 2011).Factors associated with DNA hypomethylation include numbers of cell divisions and late replicating loci (Zhou et al. 2018).Progressive loss in DNA methylation occurs predominantly in late replicating heterochromatin with DNA associated with the nuclear lamina characterized by low gene and GC density (Zhou et al. 2018).DNA demethylation proximal to transposable elements was greater than global demethylation (Kong et al. 2019), which is consistent with the global and L1l demethylation observations in our study.
The RRBS data at SEN and SENbp displayed a slight but statistically significant elevation in global genome DNA methylation after 48 days of culture, followed by a small reduction in methylation in cells that sustained proliferation.Hypothetically, the small increase in DNA methylation at SEN might be linked to a transient increase in Dnmt3l expression (Yu et al. 2020), and/or redistribution of repressive histone marks (H3K9me3, H3K27me3) by Rb-E2f-dependent methyltransferases (Yu et al. 2018) with associated Dnmt complex (Rose and Klose, 2014;Yang et al. 2022).The subsequent small decrease in DNA methylation may be associated with the SENbp upregulation of H3K9me3 demethylases (LSD1, JMJD2C; (Yu et al. 2018)), and/or elevated oxidative stress (Zhou et al. 2016).Oxidative stress induces the formation of 8-oxo-7,8-dihydroguanosine (8-oxo-dG), which is recognized by 8-oxo-dG DNA glycosylase (OGG1) to initiate DNA repair.In addition to being mutagenic, 8-oxo-dG was suggested as an epigenetic transcription marker to regulate gene expression and to stimulate demethylation of cytosines adjacent to 8-oxo-dG (Giorgio et al. 2020;Zarakowska et al. 2014).A model was proposed by which OGG1 bound to 8-oxo-dG recruits the DNA demethylase TET1, which then initiates cytosine demethylation (Zhou et al. 2016).Alternative mechanisms such as the availability of cofactors for methylases and demethylases may contribute to methylation changes.The specific mechanisms responsible for these changes in indices of global genome DNA methylation at SEN and SENbp need further investigation.Overall, a decrease in heterochromatin (L1l assay) DNA methylation appears to occur spontaneously (in B[a]P exposed and non-exposed cultures) with cumulative doubling in vitro.The constant decline in L1l methylation observed in our study is of sufficient magnitude to propose that future experiments should test if a hypomethylation threshold can be used to predict key events of cancer AOPs, such as enhanced mutagenic L1 expression, DNA repair deficiency, or genomic instability.However, such DNA methylation declines in vitro may be dependent on cell types and cell culture conditions.L1 expression differ Fig. 7. Relative DNA methylation abundance in six differentially methylated region-bearing genes measured by pyrosequencing (Bhlhe22, Aifm3, Pou4f1, Bgalnt2, Gja8, Klf17, are also listed in Table 2).The data is presented through time (passage number) relative to the average of the DMSO treated samples (D in black).The methylation pattern through time is generally reproducible across clones and laboratories.The regression lines show DNA methylation changes across the Health Canada (HC) samples in blue, including clone HC26d1.Brunel University (Pickles et al. 2016) provided samples from two laboratories, these include 12 B[a]P-exposed clones shown in red and numbered 1-12.A red line was drawn when more than one sample from the same clone were collected at different passage numbers.A second set highlighted in black include samples of naïve cells (N), DMSO treated samples (D), and an unusual DMSO-exposed MTc (illustrated with two samples connected with a black line).The HC samples at passage 0 are naïve cell samples 24-48 h after tissue dissection.across cell types, for example it is highly expressed in MCF7 cells but not in Hela or HEK293 cell lines (Freeman et al. 2022).Other DNA repeated sequences may be more sensitive to toxic insults than L1.Investigating DNA methylation in primary human fetal embryonic cells collected from amniocentesis and exposed to X-ray radiation, SAT-α DNA methylation showed a larger dynamic range than ALU or L1 at 20 population doubling (the only time-point), with no statistically significant effect on L1 (Flunkert et al. 2018).Overall, DNA methylation of retrotransposons and of other DNA repeats as cell type specific biomarkers of genetic stability deserve further investigation.

Conclusion
Starting with primary cultures of SH fetal cells exposed initially for 7 days to B[a]P, this work established the chronology of DNA methylation changes up to the time at which the cells bypassed senescence and then sustained proliferation in vitro for more than 213 days of culture.Furthermore, the study identified early DNA methylation changes (19-25 days of culture) that persisted and that were found in colonyderived cultures obtained from a different laboratory, which confirms the reproducibility of the data.A challenge is to ensure that what is an apparent SH carcinogenic step in vitro is also relevant to human in vivo.Further investigations of DMRs reported herein can provide mechanistic understanding of molecular events leading to SEN and SENbp and eventually may be used in new chemical testing strategies to support early key events (e.g., genetic instability, DNA repair deficiency, immune evasion) of cancer AOPs relevant to humans, thereby strengthening predictions of chemicals that can increase the risk of developing cancers.

Declaration of Competing Interest
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.
: Ctrl vs B[a]P, number of filtered genes = 43, Ctrl vs Nc n = 53, Ctrl vs MTc n = 108, Nc vs MTc n = 122; Fig. 5: Ctrl vs SEN n = 380, Ctrl vs SENbp n = 554, SEN vs SENbp n = 388),and the gene lists for each sections of the Venn's diagrams are provided in sheets of the supplementary Excel file.Interestingly, Fig.4indicates that most DMR-bearing genes in the contrast Ctrl vs B[a]P were transiently affected (36/43) and no longer detectable in the 19-25 day cultures necessary to generate the Nc and MTc samples.Also a list of genes (regions E plus F) showed DMRs that can differentiate Nc from MTc at 19-25 doc.DMR-bearing genes that may be involved in the mechanisms of SEN and SENbp are identified in Fig.5.The control vs SEN contrast revealed 380 genes with DMRs that may be induced either by the cell culture delay from 9 to 48 days required to reach SEN, and/or by the mechanisms inducing SEN.Interestingly, 171 genes (B:135 + E:36) still bear DMRs in the SENbp samples collected at 94 doc, which support the persistence of many DMRs from Day 48 up to 94 doc.The contrast SEN vs SENbp revealed 230/388 genes bearing DMRs that may distinguish these two events.

Fig. 2 .
Fig. 2. Global genome DNA methylation derived from RRBS analyses of cells exposed in vitro for 7 days to dimethylsulfoxide (DMSO, 0.2%) or to 20 µM benzo[a]pyrene (B[a]P).See protocol (Fig. 1).Four samples were analysed by RRBS in each group.The level of methylation in one normal colony (Nc) was 22.6% and was not considered in this graph.The Nc and morphologically transformed colonies (MTc) were allowed to proliferate until 19-25 days of culture (doc) to collect enough cells for molecular analyses.A culture in senescence (SEN, 48 doc) gave rise to the clone HC26d1 collected after senescence-bypass (SENbp, 94 doc) and later during sustained proliferation (SUSTp, 213 doc).Means with different letters are significantly different, Tukey-Kramer HSD p < 0.05.

Fig. 3 .
Fig.3.Diseases, biofunctions, and canonical pathways enriched with genes bearing differentially methylated regions in their promoter, as identified by IPA.To ensure readable text in the figure, the list of diseases and biofunctions was cut-off at B-H adjusted p < 10 − 8.5 and that of canonical pathway at B-H adjusted p < 0.01.Interestingly, of the 29 diseases and biofunctions listed, 26 were cancer related.

Fig. 4 .
Fig. 4. Venn's diagram representing the distribution of the number of differentially methylated regions-containing genes across the contrasts control (Ctrl) vs B[a]P treated cells, Ctrl vs normal colonies (Nc), Ctrl vs morphologically transformed colonies (MTc), and Nc vs MTc.The lists of genes for the region A, B, C, D, E, F, and G are provided in the supplementary Excel file.

Fig. 5 .
Fig. 5. Venn's diagram representing the distribution of the number of differentially methylated regions-containing genes across the contrasts control (Ctrl) vs senescence (SEN), Ctrl vs senescence-bypass (SENbp), and SEN vs SENbp.The lists of genes for the region A, B, C, D, E, F, and G are provided in the supplementary Excel file.
Fig. 6.Gene networks identified by IPA that include some of the 24 DMR-bearing genes detected in 4 or 5 contrasts.A) IPA pathway entitled "Cell Cycle, Cellular Development, Organ Development", including 11 DMR-bearing genes shown in grey out of 35 molecules (IPA score=25) (the score is the -log 10 p-value of the probability of finding 11 genes in a set of 35 genes from the IPA Knowledge Base Global Molecular Network.Interesting p-value are typically low e.g. 10 − 8 (https:// resources.qiagenbioinformatics.com/white-papers/IPA-netgen-algorithm-whitepaper.pdf)).B) IPA pathway entitled "Cellular Development, Cellular Growth and Proliferation, Hematological System Development and Function", including 10 DMR-bearing genes shown in grey out of 35 molecules (IPA score=22).Solid and dashed lines represent direct and indirect interactions, respectively.Molecule symbols are described with Fig. 6 B.

Table 2
Percentages of DNA methylation differences in differentially methylated regions (DMRs) across cellular events relative to the DMSO control group at 9 days of culture.Negative or positive numbers indicate hypo or hypermethylation, respectively.

Table 3
Distribution of the 253 categories of diseases and functions identified by IPA (B-H adjusted p < 0.05) using the list of 24 differentially methylated region (DMR)bearing genes.Categories that included from 5 to 22 DMR-bearing genes all had cancer as a common topics.The DMR-bearing genes are presented in two sheets of the supplementary Excel file (one sheet for all diseases, the other for cancers only).