General study design

In the dose-response experiment, A549 cells were treated with three doses of DEP for four weeks, using 2–3 biological replicates per group. The DNA methylation levels of 850 K CpG sites were assessed for each treated group, and a sufficient dose of DEP was identified based on the results of differentially methylated analyses compared to the control.

In the main experiment, A549 cells were treated with the previously identified DEP dose for four weeks in four biological replicates, independent of the dose-response experiment. Differential methylation and differential expression analyses were used to estimate methylation levels of 850 K CpG sites and genome-wide transcriptional changes, respectively. The biological effects of DEP exposure on the methylome were evaluated by colocalizing differentially methylated CpG sites with other epigenetic marks and transcription factor binding sites. The impact of transcriptional changes was evaluated by assigning differentially expressed genes to the 51 gene sets representing specific, well-defined biological states or processes, as well as the ENCODE transcription factor target gene sets. To assess the coordination of DNA methylation and transcriptional changes, correlations were calculated between the DNA methylation level at individual differentially methylated cis-CpG sites and the expression level of the corresponding differentially expressed gene.

DEP collection and preparation

The same batch of DEP as used in the previous study [24] was used for all subsequent experiments. To briefly describe the DEP preparation, a 3-L common rail direct injection diesel engine with oxidation catalyst and exhaust gas recirculation system was used to generate DEP. The engine was operated at 2,040 rpm and 35 Newton-meters on a 220 kW ECDY dynamometer (Meiden-Sya, Tokyo, Japan) using a fuel (LS light oil; ENEOS, Tokyo, Japan) containing 27 ppm sulfur and lubricating oil. DEP were collected from the accumulation in the dilution tunnel of the diesel exhaust generating system. The particle size distribution of the DEP based on particle number was analyzed with a scanning mobility particle sizer (model 3938; TSI, Inc., MN, USA), and the particle size was normally distributed with a peak at 100 nm in diameter and a maximum value of less than 800 nm. To measure the concentrations of polycyclic aromatic hydrocarbons (PAHs) in the DEP used, the DEP were subjected to Soxhlet extraction with dichloromethane. The soluble fraction of the DEP was about 10%. The contents of PAHs in the soluble fraction of the DEP were analyzed by fluorescence high-performance liquid chromatography (Shimadzu, Kyoto, Japan), and the concentration of each measurable PAH such as fluoranthene in the DEP was calculated (Additional file 1: Table S1). DEP were stored at − 80 °C until the experiments were performed. The stored DEP were thawed before use and suspended in DMEM/F12 medium containing 10% fetal bovine serum (FBS) and 10 µg/ml gentamicin) for cell exposure.

Cell culture and treatment

A549 cells (ATCC, Manassas, VA, USA) were cultured in flasks containing DMEM/F12 medium supplemented with 10% FBS and 10 µg/ml gentamicin under humidified 5% CO2/95% air at 37 °C. At cell passages, cells adhering to the bottom of the flasks were detached with 0.25% trypsin/1 mM EDTA solution (Invitrogen, Carlsbad, CA, USA). A549 cells were seeded in T-25 cell culture flasks at a cell density of 0.25 × 105 cells/ml and then exposed to DEP-suspended culture medium containing DEP at concentrations of 0, 3, 10, or 30 µg/ml. After one week, cells (at a density of 3–5 × 106 cells/ml) were collected with the trypsin/EDTA solution and reseeded into a new T-25 cell flask at a density of 0.25 × 105 cells/ml. The reseeding procedure was repeated a total of three times over a four-week DEP exposure period. The DEP-suspended culture medium was replaced at the time of reseeding and four days after seeding/reseeding. After four weeks of DEP exposure, cells were harvested for use in subsequent experiments.

DNA isolation and DNA methylation arrays

DNA isolation was performed from each of the A549 cells exposed to different concentrations of DEP for four weeks (total of 18 biological replicates). The experimental design was as follows: dose-response experiment: 0 µg/ml (n = 3; control), 3 µg/ml (n = 3), 10 µg/ml (n = 2), 30 µg/ml (n = 2); main experiment: 0 µg/ml (n = 4; control), 30 µg/ml (n = 4). The genomic DNA of the cells exposed to DEP for 4 weeks was extracted using the DNeasy Blood & Tissue Kit (Qiagen, Hilden, Germany). The concentration of the extracted genomic DNA was measured using a Qubit fluorometer (Thermo Fisher Scientific, Waltham, MA, USA). Epigenome-wide DNA methylation profiling was performed using the Infinium MethylationEPIC BeadChips (Illumina, San Diego, CA, USA), according to the manufacturer’s protocol (for the dose-response experiment) or by outsourcing to Macrogen Japan (Kyoto, Japan) (for the main experiment). Briefly, after performing bisulfite treatment using the EZ DNA Methylation Kit (Zymo Research, Irvine, CA, USA), the processed DNA samples were hybridized, and fluorescently labeled on the Infinium MethylationEPIC BeadChips. The labeled beads were visualized using the iScan system (Illumina, San Diego, CA, USA), and the DNA methylation status of each sample was quantified for each target-specific probe based on the fluorescence intensity.

Differential methylation analysis

Raw microarray data in idat format were imported into the R programming environment (v.4.1.0) and processed using the minfi package (v.1.38.0) [25]. Annotations for each CpG site were added using the getAnnotation() function and the IlluminaHumanMethylationEPICanno.ilm10b4.hg19 package. The quantile method was used to normalize the data. As a quality control (QC), probes with detection p-value > 0.01 in one or more samples, probes on sex chromosomes, and probes with a single-nucleotide polymorphism at the CpG site were excluded from subsequent analyses. The β-value and M-value for each CpG site were obtained using the getBeta() and getM() functions, respectively. Principal component analysis and its visualization were performed using the plotMDS() function with the “gene.selection = common” parameter. Probe-wise differential methylation analysis for each pair of treated groups was performed using the R package limma (v.3.48.0) [26]. The M-value and batch-adjusted design matrices were subjected to the lmFit() function to fit the linear model. The resulting model was fit to a given set of contrasts with the contrasts.fit() function, followed by the empirical Bayes adjustment eBayes() function to obtain p-values for the respective CpG sites. P-values were adjusted for multiple testing using the Benjamini–Hochberg procedure. An adjusted p-value < 0.05 was considered statistically significant. In addition, Bonferroni correction was applied in the dose–response experiment with three pairwise comparisons (control vs. groups exposed to 3, 10, 30 µg/ml DEP). The final threshold for the adjusted p-value in the dose–response experiment was set at 0.05/3 = 0.016667.

RNA isolation and RNA sequencing (RNA-seq)

RNA was isolated from the cells exposed to DEP for four weeks in the main experiment using the RNeasy Mini Kit (Qiagen, Hilden, Germany) following the manufacturer’s protocol (control (n = 4); DEP 30 µg/ml (n = 4)). The extracted total RNA was quantified and qualified with the NanoDrop (Thermo Fisher Scientific), the Qubit RNA Assay (Thermo Fisher Scientific) and the TapeStation RNA ScreenTape (Agilent Technologies, Inc., Santa Clara, CA, USA). RNA-seq libraries were generated from 500 ng of total RNA for samples with RNA integrity numbers > 8. Briefly, the NEBNext Poly(A) mRNA Magnetic Isolation Module and the NEBNext Ultra II Directional RNA Library Prep Kit for Illumina (New England BioLabs, MA, USA) were used for polyA RNA selection and subsequent RNA-seq library preparation, respectively, according to the manufacturer’s instructions. The RNA-seq libraries from each sample were pooled in equimolar amounts and sequenced on the Illumina NovaSeq 6000 (Illumina), yielding approximately 20 M paired-end reads (2 × 150 bp) per sample. Library preparation and sequencing were conducted by Azenta Life Sciences (Tokyo, Japan).

Differential expression analysis

The raw fastq files obtained from RNA-seq were trimmed to remove adapters and bases with < 20 quality. STAR aligner [27] was used with the Ensembl_Homo_sapiens.GRCh37 reference genome to align the filtered reads and generate alignment files (bam files). Alignment metrics for the bam files were evaluated using the Samtools software (v.1.3.1) (http://samtools.sourceforge.net), and the results were included into Additional file 4: List 1. HTseq-union was used for gene counting [28]. Genes with counts per million less than 0.5 in more than 4 samples were filtered out from the subsequent differential expression analysis. Afterwards, trimmed mean of M-values (TMM) normalization [29, 30] and log2 transformation were applied to obtain a normalized gene count (expression) matrix. Differential analysis was performed using the limma workflow with trend adjustment. The genes with Benjamini–Hochberg adjusted p-value < 0.05 and |log2 fold-change (logFC)| > 0.6 were considered to be differentially expressed.

Verification of RNA-seq results using microarray analysis

To ensure reproducibility of the transcriptional changes in the main experiment, total RNAs for microarray analysis were extracted from five samples of A549 cells exposed to 0 (n = 2) or 30 µg/ml (n = 3) concentration of DEP for four weeks using the RNeasy Mini Kit (Qiagen, Hilden, Germany), independently of the main experiment. For each sample, total RNA was converted to cyanin-3-labeled cRNA using the Low Input Quick Amp Labeling Kit (Agilent Technologies) according to the manufacturer’s instructions. The labeled cRNA was used for hybridization to the Agilent SurePrint G3 Human GE Microarray 8 × 60K (v.3.0). After hybridization, the microarray slide was washed, and then scanned with an Agilent DNA microarray scanner according to the manufacturer’s protocol. Fluorescence signal intensity was quantified from the scanned image using Feature Extraction software (Agilent Technologies). The resulting signal intensity files were imported into the R programming environment (v.4.1.0) and processed using the limma package (v.3.48.0) [26]. Briefly, signal intensities were background corrected (method = ‘normexp’), normalized (method = ‘quantile’), filtered for failed probes, and used for further differential expression analysis. The differential expression analysis was performed in the same way as the probe-wise differential methylation analysis above, using linear model fitting and empirical Bayes adjustment to obtain p-values for the respective genes. The genes with Benjamini–Hochberg adjusted p-value < 0.05 and |logFC| > 0.6 were considered to be differentially expressed.

Enrichment analysis

Region-based enrichment analysis was performed using the runLOLA() function of the LOLA R package (v.1.24.0) [31]. This function used one-sided Fisher’s exact test and required three files as input: query regions of interest, universe region set, and region set database. For regulatory region enrichment analysis and TFBS enrichment analysis, the list of genomic positions of the examined CpG sites (all differentially methylated CpGs (DMCs); hypomethylated DMCs (∆β < 0); hypermethylated DMCs (∆β > 0); positively correlated DMCs; negatively correlated DMCs) was provided as query regions of interest. As the universe region set, we used either all QC-passed CpG sites or all DMCs identified in the main experiment. Genomic positions of the regions tested for enrichment were extracted from various databases. Specifically, the genomic positions of CpG islands were downloaded from UCSC ftp using the hg19 genome reference. The genomic positions of DNase hypersensitive sites (DHSs) broad peaks for A549 cells were obtained from the ENCODE database (GEO accession numbers: GSM736580, GSM736506). To call the genomic positions of the peaks associated with H3K9me3, H3K9ac, H3K4me3, H3K4me1, H3K36me3, H3K27me3, H3K27ac marks in A549 cells, raw ChIP-seq fastq files were downloaded from European Nucleotide Archive (study accession number: PRJDB2453; run accession numbers: DRR016932, DRR016933, DRR016934, DRR016935, DRR016936, DRR016937, DRR016938, DRR016939 (input)). Subsequently, the ChIP-seq data were pre-processed, and the genomic positions of the broad or narrow peaks below the q-value threshold of 0.01 were called using the MACS3 [32]. The genomic positions of the TFBSs in A549 cells were downloaded from the ENCODE database (wgEncodeRegTfbsClusteredWithCellsV3_hg19.bed). After applying the runLOLA() function to the entire set of files, regions of interest with a q-value less than 0.05 were considered as DMC enriched.

To compare the transcriptional changes obtained in the main and verification experiments, Gene Set Enrichment Analysis (GSEA) [33] was performed using the GSEA software (v.4.1.0) with default settings except for the following parameters: “Collapse/Remap to gene symbols” was set to “No collapse” and “Permutation type” was set to “Gene set”. The input data matrix consisted of normalized log2(TMM) gene expression values obtained from the main RNA-seq experiment. Significantly up- or downregulated genes obtained in the verification experiment were used as predefined gene sets. The degree to which genes in the up- or downregulated gene set were overrepresented at the top or bottom of the logFC-based ranked list of genes from the RNA-seq experiment was assessed using the enrichment score (ES) and its probability at a 5% significant level.

The WebGestalt tool [34] was used to perform gene-based overrepresentation analysis. In brief, the GMT file containing the Hallmark gene sets [35] was downloaded from the Human Molecular Signatures Database (https://www.gsea-msigdb.org/gsea/msigdb/collections.jsp), and then the Senescence gene set [36] was added to the downloaded GMT file. The GMT file containing the ENCODE ChIP-seq transcription factor target gene sets was downloaded from the ChEA3 website (https://maayanlab.cloud/chea3/). Each GMT file, along with lists of differentially expressed genes (DEGs) and reference genes (14,604 genes that passed quality control), was uploaded to the WebGestalt for the overrepresentation analysis. The minimum and maximum number of genes for the gene set were set to 5 and 2000, respectively. Gene sets with Benjamini–Hochberg adjusted p-value of less than 0.05 were considered to be significantly enriched.

Correlation analysis

Correlations between DNA methylation changes in the main and dose-response experiments, as well as transcriptional changes in the main and verification experiments were examined using the cor.test() R function with the “method = spearman” parameter.

To assess the correlation between DNA methylation and gene expression in the main experiment, cis-CpG–gene pairs were generated based on the genomic locations of the CpG site and the transcript start site (TSS) of the gene. The location of the CpG site was determined based on the Illumina Human Methylation EPIC annotation; the TSS location of the gene was determined based on the Ensemble_Homo_sapiens.GRCh37.87.gtf annotation, taking into account the direction of the transcription. A cis-CpG–gene pair was formed only if the TSS of the gene was located within a 1-Mbp window to the left or right of the CpG site location. A CpG site paired with a gene in this manner was termed a cis-CpG. To calculate the correlation between the DNA methylation level and gene expression of each cis-CpG–gene pair, the cor.test() R function with “method = kendall” parameter was used for a set of DNA methylation (M-values)–expression (normalized gene counts) pairs obtained from eight samples (four controls and four DEP-treated). To find significantly correlated pairs, the significance threshold was set at a nominal p-value < 0.001 and the correlation coefficient |rk| > 0.9. For significantly correlated cis-CpG–gene pairs, cis-CpG sites with a correlation coefficient of rk > 0.9 or < -0.9 for all paired genes were considered as positively or negatively correlated, respectively.

Additional statistical tests

One-way analysis of variance (ANOVA) was performed using the GraphPad Prism version 9.0.1 for Windows (GraphPad Software, Boston, Massachusetts, USA). Linear regression analysis was performed using the lm() R function with default settings.

Heatmaps for DNA methylation and gene expression data were generated using the Heatmap() R function with default settings.