The role of DNA methylation on gene expression in the vertebrae of ancestrally benzo[a]pyrene exposed F1 and F3 male medaka

ABSTRACT

Benzo[a]pyrene (BaP) is ubiquitously present in the aquatic environment and has been identified as a bone toxicant. Previous studies have demonstrated that ancestral BaP exposure can cause transgenerational bone deformities in fish. Transgenerational effects are thought to be caused by heritable epigenetic changes, such as DNA methylation, histone modification, and non-coding RNAs. To investigate the role of DNA methylation in BaP-induced transgenerational skeletal deformities and the related transcriptomic changes in deformed vertebrae, we examined the vertebrae of male F1 and F3 medaka fish using high-throughput RNA sequencing (RNA-seq) and whole-genome bisulphite sequencing (WGBS). The histological results revealed that osteoblast numbers at the vertebral bone decreased in the BaP-derived F1 and F3 adult males in comparison with the control group. Differentially methylated genes (DMGs) associated with osteoblastogenesis (F1 and F3), chondrogenesis (F1 and F3), and osteoclastogenesis (F3) were identified. However, RNA-seq data did not support the role of DNA methylation in the regulation of genes involved in skeletogenesis since there was very little correlation between the level of differential methylation and gene expression profiles related to skeletogenesis. Although DNA methylation plays a major role in the epigenetic regulation of gene expression, the dysregulation of vertebral gene expression patterns observed in the current study is most likely to be mediated by histone modification and miRNAs. Notably, RNA-seq and WGBS data indicated that genes related to nervous system development are more sensitive to ancestral BaP exposure, indicating a more complex transgenerational phenotype in response to ancestral BaP exposure.

KEYWORDS:

• Benzo[a]pyrene
• DNA methylation
• transgenerational effects

Introduction

Benzo[a]pyrene (BaP) is a representative polycyclic aromatic hydrocarbon (PAH) ubiquitously present in the aquatic environment. Many studies have demonstrated its effects on development, reproductive function, neurobehavior, and immunity in both humans and fish (reviewed in [1] [2–4]). Notably, an increased incidence of skeletal defects (spinal and craniofacial) has been reported in feral fish from PAH-polluted waters [5]. Moreover, parental BaP exposure resulted in bone toxicity (altered body shape) in the offspring of medaka (F1) and zebrafish (F1 – F2) (10–1012 μg/g diet; 0.2–125 μg/L) [6,7]. Recent studies in our laboratory have further demonstrated that an environmentally realistic level of BaP (1 μg/L) caused transgenerational skeletal deformities (larvae: bone compression, reduced osteoblast number; adult: bone thinning) in both developing larvae and adult males (F1 – F3) following ancestral BaP exposure [8–11]. The in vivo expression of osteoblast and osteoclast marker genes (collagen 10a1 (COL10a1) and osterix (OSX)) were disrupted in developing larvae [8], while bone morphogenetic proteins (BMP2 and 3) and genes involved in the WNT/β-catenin pathway were possibly responsible for the dysregulation of bone formation and resorption processes in the vertebrae of adult males upon ancestral BaP exposure [9,10].

Transgenerational effects occurring in the F3 generation and beyond should be distinguished from multigenerational effects occurring in the F1 and F2 generations (see review by [12]). In medaka fish, which reproduce via external fertilization, the F1 offspring are directly exposed to a toxicant via F1 germ cell exposure in the gonads of F0 parents and/or possible chemical transfer into the yolk sac. Moreover, F2 offspring may also be directly exposed since the chemical deposit in an F1 yolk sac may release and affect F2 primordial germ cells (PGCs) harboured in the developing F1 embryos [11,13]. On the contrary, F3 offspring are never exposed to the chemical, and any effects observed in the F3 and subsequent generations are likely attributable to epigenetic inheritance. Considering the difference between multigenerational and transgenerational effects, the transgenerational vertebral deformities of F3 adult males can be attributed to epigenetic modifications inherited from the exposed ancestors. To explore the epigenetic modifications potentially responsible for BaP-induced transgenerational bone impairment, five functional pairs of mRNA/miRNA (i.e., OSX/miR−214, col2a1b/miR−29b, Runt-related transcription factor 2(RUNX2)/miR−204, SRY-box transcription factor 9 (SOX9b)/miR−199a−3p, and Adenomatous polyposis coli (APC)/miR−27b) were suggested to be involved in decreased osteoblast activity in the F3 generation [10]. Additionally, mRNA and miRNA transcriptomic analysis identified a few dysregulated mRNA/miRNA sets that may explain the inhibition of bone formation and promotion of bone resorption in F3 vertebrae [9].

Gene expression can also be regulated by DNA methylation and histone modifications. As one of the most-studied epigenetic mechanisms, it has been demonstrated that BaP exposure modifies global and gene-specific DNA methylation in fish in concomitance with developmental, reproductive, and neurobehavioral toxicity. Fang et al. [14] reported that the waterborne BaP exposure (24 μg/L) of zebrafish embryos significantly decreased global methylation levels and decreased global DNA methylation and aberrant changes in the promoters of many human disease-related genes were identified following parental and continued embryonic BaP exposure at 50 μg/L [6]. Furthermore, Knecht et al. ([15]) demonstrated that the developmental exposure of F0 zebrafish embryos to BaP (1262 μg/L and 2514 μg/L) could cause the transgenerational inheritance of neurobehavioral and physiological deficits in subsequent generations (hyper locomotor activity, decreased heartbeat, and mitochondrial function in F0 and F2 juvenile). Reduced global DNA methylation coupled with the downregulated expression of dnmt genes was observed in F0 embryos, which suggests that DNA methylation may play a role in the transgenerational neurophysiological effects resulting from BaP exposure. However, methylation data for F1 and F2 generations were not investigated. Deregulation of DNA methylation was highlighted in the sperm of the F1 generation upon larval neurobehavioral and-developmental impacts in response to parental BaP exposure [16].

Although BaP-induced DNA methylation changes influence the transcription levels of genes that play important roles in various physiological functions, little is known about their involvement in vertebral malformations in fish. Nevertheless, we have recently reported that BaP exposure in F0 adult fish caused transgenerational vertebral deformities in F1, F2, and F3 adult males. Therefore, it is hypothesized that DNA methylation changes are involved in the regulation of gene expressions in vertebrae and responsible for the observed vertebral deformities in male fish. More importantly, some DNA methylation marks might be conserved between generations. Although skeletal deformities were present in the vertebrae of male F1, F2, and F3 fish, the F1 and F3 generations can present multigenerational bone toxicity (directly exposed) and transgenerational bone toxicity caused by BaP, respectively.

Materials and methods

Benzo[a]pyrene

BaP (Chemical Abstracts Service no. 50-32-8; purity > 97%) was purchased from Sigma-Aldrich. A stock solution (20 000 μg/l) was prepared by dissolving BaP powder in ethanol – dimethyl sulphoxide solvent (EtOH: DMSO = 4:1), which was stored at −20°C in aliquots. The ethanol and DMSO used were analytical or high-performance liquid chromatography grade.

Medaka model

The Japanese medaka (Cab strain) used in this study originated from Prof. Christoph Winkler’s laboratory at the Department of Biological Sciences, National University of Singapore, and have been maintained in our laboratory since 2015. The breeding pairs (∼6 months old) used in this experiment were raised in charcoal-dechlorinated tap water at a constant temperature (26 ± 1°C) with a 14:10 h light:dark cycle.

BaP exposure

Briefly, 10 pairs of sexually mature medaka (six months old, F0 generation) were randomly assigned into glass aquariums containing 20-L of charcoal-dechlorinated tap water. Four tanks were assigned as solvent controls by adding a 0.0005% EtOH:DMSO (vol/vol) mixture and six tanks were treated with BaP by adding 1 μg/L BaP with a final concentration of EtOH:DMSO at 0.0005%. Thereafter, 75% of the water from both control and treatment group tanks was changed and replenished every 2 days. Fish were fed twice per day with hormone-free dry flake food (Marubeni Nisshin Feed Co., Ltd., Tokyo, Japan) and once per day with brine shrimp (Brine Shrimp Direct, CA, USA). BaP exposure was then suspended after 21 days. To eliminate the BaP from tanks, a total of five water changes were conducted on days 22 and 23. Thereafter, embryos were collected from the spawning females within two hours after the beginning of the light photoperiod. Collected embryos were assessed under a dissection stereomicroscope and only fertilized embryos were kept and raised to adults (F1 generation). For the assessment of F1 adult male bone development, 6-month-old males were collected from both the control and BaP treatment groups. The same procedure was repeated for the F2 and F3 generations. All animal handling was performed following the regulations of the Animal Ethics Committee, City University of Hong Kong.

DNA library preparation for whole-genome bisulphite sequencing

The genomic DNA (gDNA) of vertebrae from male F1 and F3 fish was extracted for WGBS using the phenol-chloroform-isoamyl alcohol method. For each gDNA sample (n = 1; gDNA of 15 fish), DNA quantity, integrity, and purity were determined using 1% agarose gel electrophoresis and a Qubit DNA Assay on a Qubit 3.0 fluorometer (Life Technologies, CA, USA). The gDNA was fragmented into 200–300 bp fragments by sonication. Terminal repair, A-ligation, and methylation sequencing adapter ligation were then performed on the DNA fragments. DNA library bisulphite treatment was conducted using the EZ DNA Methylation Gold Kit (Zymo Research Corp, Irvine, CA, USA). The resulting single-stranded DNA was amplified using a HiFi HotStart Uracil + ReadyMix polymerase chain reaction (PCR) Kit (KAPA Biosystems, Boston, MA, USA). Thereafter, the library concentration and quality were determined using a Qubit 3.0 Fluorometer and an Agilent 2100 Bioanalyser (Agilent Technologies, Santa Clara, CA, USA), respectively. Libraries were sequenced using the HiSeq X Ten System at Novogene Co., Ltd (Beijing, China) to generate 150-bp paired-end reads delivered in FASTQ format, which provided sufficient genome coverage for differential methylation analysis.

WGBS analysis

Sequencing files that were generated by the sequencer first underwent initial quality control [17] and were then subjected to quality and adaptor trimming by using Trim Galore [17] with the following parameters: (a) removal of base calls with a Phred score of 20 or lower; (b) removal of reads shorter than 20 bp; (c) removal of reads with an adaptor sequence. After retrieving qualified reads from raw data, Bismark (v0.22.3) [18] was used for WGBS alignments (http://www.bioinformatics.babraham.ac.uk/projects/bismark/). Before alignments could be performed, the genome of Japanese medaka (Oryzias latipes, genome assembly: ASM223467v1) was bisulphite-converted in-silico and indexed using the bismark_genome_preparation command. Filtered reads were aligned to the bisulphite-converted genome using the bismark command, which takes advantage of the Bowtie2 aligner. Generated Bismark alignment BAM files underwent a deduplication step using the deduplicate_bismark command, which removed all reads but one, which was aligned to the very same position and orientation. Additionally, alignments with a mapping quality (MAPQ) value of less than 10 were discarded. To extract the methylation call for each cytosine (C), the bismark_methylation_extractor command was used on the result files. The position of every C was then written out to a new output file depending on its content (i.e., cytosines (C) followed by guanine (G) residues (CpG), C followed by any base except G and G (CHG), or C followed by any base except G (CHH)), and only CpG content was kept for downstream analysis.

The downstream analysis and generation of figures for methylation data were conducted in R [19]. The Bioconductor package DSS (Dispersion Shrinkage for Sequencing data) was used for differential methylation assessment [20]. The core of DSS is a procedure based on the Bayesian hierarchical model to estimate and shrink CpG site-specific dispersions. The minimum read depth was 30 × . Read depth variation has been taken into consideration through the DSS pipeline, which includes an algorithm to account for technical and biological variability [21] [22] Wald tests were then performed for differential methylation detection. Since only one biological sample in each treatment was sequenced, the method used to analyse the differentially methylated locis (DMLs) or regions was limited. However, one of the advantages that DSS provides is that the statistical test can be conducted even when there are no biological replicates because the neighbouring CpG sites can be considered pseudo-replicates due to smoothing, while dispersion can still be calculated with reasonable precision [23]. Data from the bismark output were prepared in the format required by DSS, which contains the following information: chromosome number, genomic coordinate, the total number of reads, and the number of reads with methylation. The callDML and callDMR functions from the DSS package were used to extract DML and regions with methylation differences greater than 10% (DMLs: p < 0.001 and DMR: p < 0.01). Differently methylated regions (DMRs) were filtered based on their length (≥50 bps) and number of CpG sites (≥4). The neighbouring DMRs were merged if the distance between them was less than 50 bps.

Gene ontology enrichment analysis (DMR-Associated genes)

Genes that had DMRs on their promoter regions (4 kb upstream and 200 bp downstream of the transcription start site (TSS)) and gene body (from TSS to transcription end site (TES)) were considered DMR-associated genes. These DMR associated genes were used in the Gene Ontology (GO) (biological process and molecular function) and Kyoto Encyclopedia of Genes and Genomes (KEGG) enrichment analysis in Metascape [24] and Ingenuity Pathway Analysis (IPA, version: 60467501, Qiagen, US).

RNA-Seq

Total RNA extracted from the vertebrae of F1 and F3 males (n = 3 and each sample was pooled from 6 fish) were sent to BGI (Wuhan, China) for RNA sequencing. A library was prepared by BGI and sequencing was conducted on an Illumina Nova-seq platform (Illumina, US). At least 70 million sequences with 150-base paired-end single-reads were generated per sample. Reads were mapped to the Japanese medaka HdrR genome (assembly version: ASM223467v1) using Hisat2 [25]. Successfully aligned reads were counted using the featureCounts programme from the Subread package [26]. Read counts were processed and analysed using DESeq2 to calculate the expression changes of all successfully aligned genes [27]. Differentially expressed genes (DEGs) were filtered by adjusted p-value (padj <0.05). GO enrichment analysis and IPA were performed on these DEGs using Metascape [24] and IPA (Qiagen, US), respectively. Transcriptome raw data for the F3 generation were retrieved from PR- JNA757563 [9].

Bone histology

Whole medaka fish (6 months old) were fixed after the removal of hard tissues (i.e., fins, operculum, skull roof, and otoliths). Fixed samples were processed according to the protocol described in Kong et al. [28]. Briefly, adult males were fixed in GPHS [0.05% glutaraldehyde, 2% paraformaldehyde, 80% Histochoice MB fixative (Cat# 64115–04, Electron Microscopy Sciences, USA), 1% sucrose and 1% CaCl₂ solution] overnight at 4°C. Fixed fish were dehydrated in a 70/80/95/100% graded series of methanol followed by cleaning in chloroform. Cleaned samples were infiltrated and embedded in melted paraffin (Fisher Scientific, US). Serial sections of individual fish were cut (7 μm) on a rotary microtome (Leica RM2125, Germany), mounted onto SuperFrost® Plus slides (Menzel-Gläser, Germany), and dried on a slide warmer at 33°C overnight. Tissue sections were observed under a compound microscope, and only sections showing vertebral columns were selected and number coded for specific staining.

Identification of osteoblasts via Goldner’s trichrome staining

Osteoblasts (OBs) were identified by measuring alkaline phosphatase (ALPase) activity using a modified version of Goldner’s trichrome staining protocol [29]. Briefly, medaka tissue sections were deparaffinized and immersed in Weigert’s haematoxylin for nuclei staining (20 minutes), followed by immersion in Ponceau/acid fuchsin/azophloxine for 5 minutes. The staining colour was differentiated with 0.5% acid alcohol for 20 seconds. The stained sections were treated with phosphomolybdic acid and light green stain for 5 and 20 minutes, respectively. Stained sections were dehydrated in a series of ethanol solutions (70, 80, 100%) and xylene (two times) and then mounted for microscopic observation and analysis.

Quantitative histological analysis of osteoblast (OB)

Six sections (n = 5 in F1) containing vertebral segments 15–29 were selected to assess OB numbers at the intervertebral segments from the control and BaP-derived offspring in both F1 and F3 generations [10].

Statistical analysis

For the osteoblast population, Student’s unpaired t-test were used to test the null hypothesis that the means of the control and BaP-derived group are equal. Significant differences are indicated by * < 0.05, ** < 0.01, and *** < 0.001. Statistical analyses were performed using R and figures were plotted using Microsoft Excel and R.

Results

DNA methylation profiling of vertebrae

Approximately 31 Gb (range 30–33) of raw data were generated from each pooled vertebrate gDNA sample from the F1 and F3 adult males (Table 1). The average number of sequenced reads for each sample was 200 million with an average read depth above 30× and the bisulphite conversion efficiency was over 99.9% (Table 1). The unique mapping rates of these reads to the bisulphite-converted medaka genome were between 62% and 65% in vertebrae. The global methylation levels for the control and BaP groups of the F1 and F3 generations were 7.47%, 7.26%, 7.32%, and 7.38%, respectively. CpG content, which presents 3% of the medaka genome [30], had the highest methylation rate with 74.6%. Only 0.27% and 0.3% of cytosines were methylated in the CHH and CHG context.

Table 1. Overview of whole-genome bisulphite sequencing status in vertebrae.

Male		Raw bases (G)	Clean bases (G)	Clean reads	Effective rate^a (%)	Q20^b (%)	GC percentage (%)	Bisulfite conversion rate (%)
F1V	Ctrl	30.2	29.8	198841826	98.7	96.9	22.6	99.9
	BaP	33.2	32.8	219043308	98.9	96.8	22.2	99.9
F3V	Ctrl	31.3	30.9	205555016	98.6	96.6	22.2	99.9
	BaP	32.3	31.6	211225764	98.0	96.0	22.3	99.9

Note: ^aEffective rates (%) indicate the ratio of clean reads/bases to the total number of raw bases. ^bQ20 indicates the percentage of bases whose Qphred > 20 (inferred base call accuracy is higher than 99%).

BaP-induced multigenerational changes in DNA methylation in the F1 vertebrae

A total of 2093 DMRs were identified in the F1 vertebrae in response to BaP exposure of the F0 adults. 826 DMRs were annotated to a promoter region or a gene body (Figure 1a). 217 genes had a differently methylated promoter region (hypermethylated: 118 and hypomethylated: 99) and 608 genes showed a changed methylation pattern in the gene body (hypermethylated: 342 and hypomethylated: 266; Figure 1a).

Figure 1. Number of methylated genes in vertebrae. Number of genes with differential methylation region on promoter and gene body in the vertebrae of F1 (a) and F3 (b).

GO analysis of the differentially methylated genes (DMGs) revealed enrichment of biological processes in different categories. The enriched processes in the F1 male vertebrae associated with skeletogenesis were chondrocyte differentiation (14 genes, p = 2.07E–06), skeletal system development (33 genes, p = 7.81E–06), cartilage development (18 genes, p = 8.98E–06), regulation of chondrocyte differentiation (9 genes, p = 1.07E–05), regulation of cartilage development (10 genes, p = 2.07E–05), chondrocyte development (5 genes, p = 1.40E–03), bone development (14 genes, p = 1.70E–03), negative regulation of chondrocyte differentiation (4 genes, p = 4.12E–03), skeletal system morphogenesis (14 genes, p = 4.16E–03), negative regulation of cartilage development (4 genes, p = 8.48E–03), ossification (27 genes, p = 6.15E–05), regulation of ossification (12 genes, 1.09E–04), osteoblast differentiation (16 genes, p = 9.69E–04), regulation of osteoblast differentiation (11 genes, p = 1.36E–03), negative regulation of ossification (5 genes, p = 4.75E–03), positive regulation of chondrocyte differentiation (4 genes, p = 1.27E–03), positive regulation of cartilage development (5 genes, p = 1.40E–03) and cartilage condensation (4 genes, p = 2.41E–03; Figure 2a).

Figure 2. GO terms and ingenuity pathway analysis of DMGs in the F1 (a, b) and F3 (c, d) vertebrae upon ancestral BaP-exposure. Left panel represents GO biological process terms and the right panel shows canonical pathways. Go term analysis was performed on DMGs using Metascape and significant terms were selected with p-value <0.01 and q-value <0.25. Canonical pathway analysis was performed on DMGs by Qiagen IPA tool and statistically significant pathways were filtered at p-value <0.05 level.

In addition, based on the IPA, DMGs were found to be significantly enriched in many canonical pathways, including AMP-activated protein kinase (AMPK) signalling (p = 1.51E–04, z-score 0.905), calcium signalling (p = 2.63E–04, z-score 0), transforming growth factor-β (TGF-β) signalling (p = 2.14E–03, z-score 1.134), xenobiotic metabolism CAR signalling pathway (p = 3.09E–03, z-score 2.324), sonic hedgehog signalling (p = 3.24E–03, z-score −2.236), xenobiotic metabolism PXR signalling pathway (p = 3.55E–03, z-score 2.84), osteoarthritis pathway (p = 4.17E–03, z-score 0.277), Wnt/β-catenin signalling (p = 8.71E–03, z-score 0), HIPPO signalling (p = 1.05E–02, z-score 0), BMP signalling pathway (p = 1.20E–02, z-score 1.89), ErbB signalling (p = 1.82E–02, z-score −0.707) and aryl hydrocarbon receptor signalling (p = 2.63E–02, z-score 0.707; Figure 2b).

BaP-induced transgenerational changes in DNA methylation in the F3 vertebrae

In the F3 vertebrae, 2396 DMRs were found after F0 BaP exposure. 213 DMRs were within promoter regions (hypermethylated: 102 and hypomethylated: 111) and 691 were inside gene bodies (hypermethylated: 336 and hypomethylated: 355) (Figure 1b).

The enriched GO terms (biological processes, Figure 2c) associated with skeletogenesis were ossification (38 genes, p = 3.65E–09), osteoblast differentiation (25 genes, p = 7.27E–08), regulation of osteoblast differentiation (13 genes, p = 2.52E–04), negative regulation of osteoblast differentiation (7 genes, p = 1.04E–03), skeletal system development (39 genes, p = 1.91E–07), chondrocyte differentiation (15 genes, p = 1.19E–06), bone development (19 genes, p = 1.63E–05), cartilage development (18 genes, p = 3.16E–05), bone morphogenesis (11 genes, p = 1.46E–04), skeletal system morphogenesis (18 genes, p = 1.94E–04), chondrocyte development (5 genes, p = 2.12E–03), regulation of chondrocyte differentiation (6 genes, p = 4.67E–03), regulation of cartilage development (7 genes, p = 5.19E–03), cartilage development involved in endochondral bone morphogenesis (4 genes, p = 9.01E–03), bone remodelling (10 genes, p = 5.74E–04), regulation of bone resorption (6 genes, p = 2.70E–03), bone resorption (7 genes, p = 3.66E–03), regulation of bone remodelling (6 genes, p = 5.69E–03), osteoclast differentiation (10 genes, p = 7.47E–04), notochord development (4 genes, p = 1.80E–03), osteoblast development (4 genes, p = 1.80E–03), biomineralization (13 genes, p = 2.29E–03), bone mineralization (10 genes, p = 3.53E–03), biomineral tissue development (12 genes, p = 5.69E–03), osteoblast proliferation (5 genes, p = 3.81E–03).

The enriched canonical pathways included AMPK signalling (p = 9.55E–07, z-score −0.243), calcium signalling (p = 3.16E–04, z-score −0.535), ErbB signalling (p = 1.07E–03, z-score 0.302), osteoarthritis pathway (p = 4.79E–03, z-score 0.632), hepatocyte growth factor (HGF) signalling (p = 3.399E–02, z-score 0.378), Wnt/Ca⁺ pathway (p = 3.80E–02, z-score 0.447) and vascular endothelial growth factor (VEGF) signalling (p = 3.89E–02, z-score 0.707) (Figure 2d), where positive z-score and negative z-score indicates the pathway is activated and inactivated, respectively.

Methylation changes across generations

To identify the transmission of DNA methylation status across generations, the DMGs between the two generations were compared. The total number of DMGs in the F3 was higher (904) than that in the F1 (825) (Figure 3). There were 21 common DMGs with changes in the promoter region methylation and 96 common DMGs with a modified methylation pattern of the gene body in the F1 and F3 generation (Table 2 & 3). Of the 21 common DMGs with methylation changes in the promoter region, 12 were hypermethylated and 9 were hypomethylated (Figure 3). Of the 96 common DMGs with methylation changes in the gene body, 37 were hypermethylated and 59 were hypomethylated (Figure 3). Functional analysis of conserved DMGs showed the enrichment of skeletogenesis-associated biological processes, including chondrocyte differentiation, cartilage development, and skeletal system development (Appendix I). The DMGs consistent in the F1 and F3 were provided in Appendix II.

Figure 3. Venn diagrams for hypermethylated, hypomethylated, and total DMGs in promoter and gene body, showing numbers of common DMGs between the F1 and F3 generations in the BaP treatment.

Table 2. Common DMGs and DEGs showing expected methylation difference and expression directions in the F1 vertebrae.

Gene	OL_Gene_ID	Log2FC	Meth.diff_Promoter^a	Meth.diff_Genebody^b
SLC43A2	ENSORLG00000005019	−1.384	.297
HIF3A	ENSORLG00000002004	−0.963		−0.332
LRRN4CL	ENSORLG00000002917	−1.065		−0.365
MBP	ENSORLG00000020727	−1.077		−0.309

Note: ^a methylation difference in promoter region between control and BaP and ^b methylation difference in gene body region between control and BaP.

Table 3. Common DMGs and DEGs showing expected methylation difference and expression directions in the F3 vertebrae.

Gene	OL_Gene_ID	Log2FC	Meth.diff_Promoter^a	Meth.diff_Genebody^b
CHSY1	ENSORLG00000026489	0.809	−.317
COTL1	ENSORLG00000015454	0.828	−.667
CST3	ENSORLG00000011963	0.548	−.505
FAM171B	ENSORLG00000025190	−0.960	.334
GNA13	ENSORLG00000017378	1.064	−.383
NDFIP1	ENSORLG00000006004	1.372	−.326
NFASC	ENSORLG00000029596	−0.719	.358
NPDC1	ENSORLG00000028127	0.616	−.465
RNF220	ENSORLG00000012587	−0.733	.406
RPS4X	ENSORLG00000009584	0.645	−.371
SPARC	ENSORLG00000011369	0.711	−.313
ADAM33	ENSORLG00000012240	−0.709		−0.286
ANK2	ENSORLG00000024491	−0.626		−0.377
CD81	ENSORLG00000015561	0.401		0.297
CPNE5	ENSORLG00000006371	−1.2569		−0.645
FREM2	ENSORLG00000003520	−0.925		−0.364
GFAP	ENSORLG00000002346	−0.756		−0.522
GPHN	ENSORLG00000011754	−0.694		−0.337
HDAC7	ENSORLG00000001083	−0.600		−0.295
LDHB	ENSORLG00000009988	0.503		0.343
LLGL1	ENSORLG00000013225	−0.475		−0.335
NDUFS2	ENSORLG00000011324	0.438		0.360
NXPE3	ENSORLG00000028212	−0.786		−0.397
SCFD2	ENSORLG00000000475	0.749		0.329
SHISA6	ENSORLG00000022792	−0.925		−0.465
SLC38A4	ENSORLG00000028149	0.502		0.335
SRCIN1	ENSORLG00000014080	−0.971		−0.349
TJP1	ENSORLG00000005217	−0.706		−0.355
TTC28	ENSORLG00000011900	−0.494		−0.348
UQCC1	ENSORLG00000002962	0.532		0.353
WNK1	ENSORLG00000013832	−0.495		−0.292

Note: ^a methylation difference in promoter region between control and BaP and ^b methylation difference in gene body region between control and BaP.

BaP-Induced transcriptional changes in the F1 and F3 vertebrae

The transcriptome profiles of F1 and F3 vertebrae were studied using RNA-seq. 23.5 M and 21.8 M quality trimmed reads were obtained from the samples of the F1 generation, which generated 2.2 Gb and 2.1 Gb clean reads for control and BaP groups, respectively. The analysis of the samples from the F3 generation resulted in 21.8 M and 22.6 M filtered reads, yielding 2.2 Gb clean reads for the F3 control and F3 BaP groups, respectively.

A total of 23,209 transcripts (out of 34,461 mRNAs in the medaka reference genome (ASM31367v1) were identified in the vertebrae of the F1 generation, of which 230 (including up-regulated: 39 and down-regulated: 191, Figure 4a) were differentially expressed (adjusted p < 0.05). The functional analysis and pathway analysis revealed diverse categories of biological processes and pathways. The selected biological processes and pathways that are involved in skeletogenesis and stimuli responses were plotted (Figure 5a & b). Only one biological process involved in skeletogenesis was enriched: osteoblast differentiation (7 genes, p = 1.15E–03) (Figure 5a). The significantly enriched IPA canonical pathways included HIPPO signalling (p = 3.63E–03, z-score NA), xenobiotic metabolism PXR signalling pathway (p = 1.41E–02, z-score 2.236) and calcium signalling (p = 1.86E–02, z-score 1.342) (Figure 5b).

Figure 4. Vertebral gene expression changes. Volcano plots showing DEGs in F1 (a) and F3 vertebrae (b). Vertical dashed lines represent the -log10(adjusted p-value) of 0.05. Red and green spots represent statistically up- and down-regulated genes, respectively.

Figure 5. GO terms and ingenuity pathway analysis of DEGs in the F1 (a, b) and F3 (c, d) vertebrae. Left panel includes GO biological process terms and the right panel reveals canonical pathways. Go term analysis was performed on DEGs using Metascape and significant terms were selected with p-value <0.01 and q-value <0.25. Canonical pathway analysis was performed on DEGs by Qiagen IPA tool and statistically significant pathways were filtered at p-value <0.05 level.

In the F3 generation, 23409 genes were identified in the transcriptome and 1094 of them were differentially expressed. Among the 1094 DEGs, 731 were upregulated and 363 downregulated (Figure 4b). The DEGs showed functional enrichment of biological processes and pathways in various categories. Among them, GO enrichment analysis revealed two skeletogenesis biological processes, cartilage condensation (5 genes, p = 5.42E–04) and skeletal system development (34 genes, p = 6.36E–04) (Figure 5c). The significantly enriched canonical pathways by DEGs included calcium signalling (p = 1.20E–08, z-score 0.688), CXCR4 signalling (p = 4.79E–05, z-score 2.5), role of osteoblasts, osteoclasts and chondrocytes in rheumatoid arthritis (p = 2.04E–04, z-score NA), role of macrophages, fibroblasts and endothelial cells in rheumatoid arthritis (p = 1.38E–03, z-score NA), ERK/MAPK signalling (p = 2.24E–03, z-score 0.277), ErbB signalling (p = 2.45E–03, z-score 1.265), AMPK signalling (p = 2.57E–03, z-score 0.707), mTOR signalling (p = 5.13E–03 z-score 1.134), RANK signalling in osteoclasts (p = 5.25E–03, z-score 3), HIPPO signalling (p = 1.32E–02, z-score −1.342), BMP signalling pathway (p = 1.32E–02, z-score 2.828), aryl hydrocarbon receptor signalling (p = 1.74E–02, z-score 1.342), xenobiotic metabolism general signalling pathway (p = 1.74E–02, z-score 0.905), sonic hedgehog signalling (p = 2.09E–02, z-score NA), TGF-beta signalling (p = 2.57E–02, z-score 2.828) and xenobiotic metabolism AHR signalling pathway (p = 3.80E–02, z-score 1.89) (Figure 5d).

Relationship between DMGs and DEGs

In the F1 generation, seven common genes were both significantly differentially methylated and expressed: the amino acid transporter slc43a2 (SLC43A2) showed a negative correlation between promotor DNA methylation and gene expression, and three genes (hypoxia inducible factor (HIF3A), LRRN4 C-terminal-like protein (LRRN4) and myelin basic protein (MBP) showed a positive correlation between gene body methylation and gene expression (Table 2).

In the F3 generation, 79 genes had significantly altered DNA methylation and gene expression. Of these, 11 genes showed a negative correlation between promotor DNA methylation and gene expression, and in 20 genes gene body methylation and gene expression were positively correlated (Table 3).

Reduced osteoblast population in F1 and F3 vertebrae

The histomorphometric assessments revealed a significantly reduced OB population after parental BaP exposure in comparison with the control group (p = 0.0006, Figure 6).

Figure 6. Representative histology images of Goldner’s trichrome staining for the detection of osteoblasts (black arrowheads), osteoid (asterisks) and mineralized bones (number sign) in adult male Japanese medaka O. latipes. (a) the locations of corresponding cells on the serial tissue sections are indicated by arrows. (b) Representative segment in F1 control. (c) Representative segment in F1 BaP. (d) number of osteoblasts at the intervertebral segments of individuals from the control and ancestral 1 μg/L BaP exposed group in F1; (n = 5; mean ± SD; * and ***: p < 0.05 and 0.001, respectively; Student’s t-test).

Discussion

Whole-genome methylation analysis revealed both CpG- and DMR-specific changes in the DNA methylation profiles of adult fish (F1 and F3) after ancestral exposure to an environmentally relevant BaP concentration. BaP exposure modified DNA methylation patterns in both the promoter regions and gene body regions of the vertebrae genome in the F1 and F3 generation. Over 100 DMGs were consistent in the F1 and F3 generation. Nine DMGs (MAPK14, GLI3, HOXB3, NFIB, SKI, CHST11, KIAA1217, CREB3I2, CST) were enriched in diverse physiological pathways including chondrocyte differentiation, osteoblast differentiation and ossification, which are associated with skeletogenesis in both generations, but did not appear differently regulated in the transcriptome.

Transcriptomic profiling uncovered enrichment of (I) skeletogenesis related pathways (F1, F3); and (II) classical BaP exposure response pathways: activation of xenobiotic metabolism PXR signalling pathway (F1), xenobiotic metabolism general signalling pathway (F3) and aryl hydrocarbon receptor signalling (F3). Although both DMGs and DEGs were enriched in similar or common pathways involved in skeletogenesis, no correlated overlap between DMGs and DEGs was conserved between the F1 and F3 generation. In the F1 generation, SLC43A2 and HIF3a are possibly associated with decreased bone mass. In the F3 generation, chondroitin sulphate synthase 1 (CHSY1), histone deacetylase 7 (HDAC7) and secreted protein acidic and rich in cysteine (SPARC) have a confirmed role during skeletogenesis and a set of 19 were somehow associated with osteoblast differentiation, function and bone mineral density. The results suggest ancestral BaP exposure could alter both DNA methylation and gene expression in the directly exposed offspring (F1) and the unexposed offspring (F3). Moreover, similar methylation patterns were observed in nine loci associated with genes involved in skeletogenesis in the F1 and F3. This corroborates the initial assumption that DNA methylation may play a role in regulating specific gene expressions associated with the transgenerational skeletal deformities upon ancestral BaP exposure. Persistence of the bone phenotype was confirmed with significantly reduced OB numbers in the ancestrally exposed F3 generation, which were previously reported [8,11].

Differential methylation in the F1 vertebrae

In the vertebrae of the F1, genome-wide methylation analysis revealed a variety of DMGs. The enriched biological processes involved in skeletogenesis (Figure 2a) are major regulators of osteoblastogenesis and chondrogenesis (Figure 7).

Figure 7. Relationship between enriched biological process (solid line box) involved in skeletogenesis and potential dysregulated pathways (dash line box; orange colour: activated pathway; white colour: neither activation nor inactivation; grey colour: no predication on status) associated with the biological process in the F1. DMGs involved in specific biological process are presented and only additional genes are indicated in following tier.

TGF-β signalling and BMP signalling are crucial for osteoblast differentiation, chondrogenesis and bone formation [31,32]. Among the DMGs involved in TGF-β and BMP signalling pathways, type I BMP receptor (BMPR1b) was found to be hypermethylated in the gene body in the F1 vertebrae. BMPR1b is indispensable for chondrogenesis, including chondrocyte condensation, proliferation, differentiation, survival, and function [33]. In addition, the transcription factor Suppressor of Mothers against Decapentaplegic 7 (SMAD7) was hypermethylated in the promoter, which interferes with both TGF-β and BMP signalling [34]. The regulation of skeletogenesis is also controlled by multiple cytokines, which requires dynamic crosstalk between TGF-β and BMP signalling and other signalling pathways, such as hedgehog, Wnt and Notch signalling. For instance, TGF-β can promote Wnt/β-catenin signalling through upregulating the expressions of Wnts (Wnt−2, −4, −5a, −7a, and −10a) [35]. The cooperation of TGF-β and Wnt stimulate the differentiation of both osteoblast and chondrocyte [36]. Among the DMGs associated with Wnt/β-catenin signalling, Wnt5a was hypermethylated in the gene body. In contrast to positive regulatory role of TGF-β, BMPs can inhibit Wnt signalling in osteoblast via BMPR1a, which negatively regulate bone mass [37]. BMP-induced osteoblast differentiation can be stimulated by notch activation [38], while TGF-β and BMP2 can coordinate expression of notch signalling proteins [39]. Jagged Canonical Notch Ligand 1 (Jag1) and Neurogenic locus notch homolog protein 1 (Notch1) were found to be hypermethylated in the gene body and they are required for triggering Notch signalling [40]. Evidence of jag1 participating in bone development in fish is limited, but the gene is required for normal trabecular bone formation in mice [41]. Both TGF-β and BMP signalling were found to be activated, which suggested a positive regulation towards osteoblastogenesis and chondrogenesis. However, histomorphometric assessments revealed a significant reduction in osteoblast numbers (Figure 6d). Compensatory changes have also been observed in rockfish with increasing BaP exposure concentrations upon direct exposure [42]. Thus, the observed elevated TGF-β and BMP signalling might present a compensatory response to re-establish bone homoeostasis, which is impacted by ancestral BaP exposure.

PAHs, including BaP, have been reported to induce expression changes of genes involved in detoxification pathways [43,44]. BaP-responsive pathways (xenobiotic metabolism CAR signalling pathway, xenobiotic metabolism PPAR signalling pathway and aryl hydrocarbon receptor signalling) were persistently activated in both the F1 and F3 generation. The constitutive androstane receptor (CAR), pregnane X receptor (PXR) and aryl hydrocarbon receptor (AHR) are evolutionary conserved cellular sensors for both endogenous and exogenous stimuli and modulate numerous genes involved in a wide range of cellular processes, such as cell growth, differentiation, metabolism and stress response [45]. The AhR plays a confirmed role in bone homoeostasis and direct exposure experiments revealed a concomitant deregulation of AhR and bone genes in fish [46, 47]. In the three BaP-responsive pathways, Glutathione S-Transferase Omega 2 (GSTO2), UDP Glucuronosyltransferase Family 1 Member A3 (UGT1a3), Superoxide dismutase 3 (SOD3) and Sulfotransferase Family 1A Member 1 (SULT1a1) were differentially methylated. All of them are part of the chemical defensome defined for model fish and their transcript changes have been reported in response to direct BaP exposure in fish [48,49].

Differential methylation in the F3 vertebrae

The enriched GO terms related to skeletogenesis in the F3 generation (Figure 2c) were classified into three categories: (I) osteoblastogenesis, (II) chondrocytogenesis and (III) osteoclastogenesis (Figure 8). During osteoblastogenesis, the inactivation of AMPK signalling may point towards negatively regulated bone formation in the F3. Previous studies have shown that osteoblast differentiation and bone mineralization could be positively regulated by AMPK activation [50,51]. The calcium signalling which is associated with osteoblast differentiation was significantly inhibited in the F3 methylome. Ca2⁺ is a prominent and ubiquitous intracellular messenger that is responsible for guiding numerous cellular processes [52]. Many of its cellular regulation, which includes osteoblast proliferation and differentiation, are orchestrated by the Ca2⁺ binding protein, calmodulin (CaM). The CaM-dependent kinase II (CAMKII) was found to interact with the osteoblast-specific transcription factor OSX, and enhanced the transcriptional activity and protein levels of OSX during osteoblast differentiation [53]. Deregulation of calcium signalling is corroborated by the hypomethylation of calcium/calmodulin-dependent protein kinase (camk1d and camk2b), which might lead to silencing of these genes. The third inactivated pathway enriched by DMGs in the F3 was VEGF signalling. VEGF can regulate osteoblast differentiation and bone formation from periosteal progenitor cells [54]. The observed DMGs providing a possible molecular explanation for the reduced osteoblast differentiation and activity in response to ancestral BaP exposure.

Figure 8. Relationship between enriched biological process (solid line box) involved in skeletogenesis and potential dysregulated pathways (dash line box; orange colour: activated pathway; blue colour: inactivated pathway) associated with the biological process in the F3. DMGs involved in specific biological process are presented and only additional genes are indicated in following tier.

Regarding chondrocytogenesis, activation of ErbB signalling was revealed. Erb-B2 Receptor Tyrosine Kinase 4 (ErbB4, hypomethylated in the gene body), one of the receptors inducing ErbB signalling, is expressed in chondrocytes and has been reported to play a role in chondrocyte growth and differentiation [55]. Ancestral BaP exposure-induced perturbed chondrocyte differentiation was suggested to lead to defective notochord sheath and further disrupted bone formation in developing larvae [11].

In terms of osteoclastogenesis, the activation of HFG signalling in the F3 generation implied that osteoclastogenesis was promoted and subsequently increased in bone resorption. HGF influences various responses in numerous cells including epithelial cells, osteoblasts and osteoclasts, with a confirmed role in osteoclastogenesis [56]. Jiezhang [8] found significant reduction of bone thickness in the male fish of F3 generation following ancestral BaP exposure. The balance between bone formation and bone resorption determines whether the skeleton undergoes normal development. The BaP-induced transgenerational bone phenotype suggest that bone homoeostasis was skewed towards bone resorption by either inhibited bone formation or promoted bone resorption or both. The number of osteoblasts was significantly reduced in the BaP-derived F3 individuals, which could partially explain the bone thinning. The here described methylation changes in the F3 vertebrae indicate that reduced bone thickness could be attributed to the downregulation of bone formation and upregulation of bone resorption, which is in-line with the transcriptome [9].

Conservation of DMGs during skeletogenesis in the F1 and F3 vertebrae

To unravel epigenetic regulation and regulated genes responsible for the bone phenotype inheritance, the data were mined for conserved methylation pattern changes. The comparison of DMGs between the F1 and F3 generation revealed 117 conserved DMGs (Figure 3). Nine of them were associated with the biological process skeletal system development (Appendix I): mitogen-activated protein kinase 14 (MAPK14), GLI family zinc finger 3 (GLI3), homeobox B3 (HOXB3), nuclear factor I B (NFIB), SKI Proto-Oncogene (SKI), carbohydrate sulfotransferase 11 (CHST11), sickle tail protein homolog (KIAA1217), cAMP responsive element-binding protein 3-like 2 (CREB3L2) and CTS telomere maintenance complex component 1 (CTC1). MAPK14 is a key factor for OSX and RUNX2 phosphorylation to regulate osteoblast activity and subsequently affects bone mineralization [57,58]. GLI3 is associated with bone development and serves as a key effector for hedgehog signalling in the developing skeleton [59]. HOXB3 has been identified as a candidate gene for bone growth and development in chicken [60]. NFIB is an important regulator of the Insulin-like growth factor-binding protein 5 expression in osteoblasts, and thus, plays a critical role in osteoblast proliferation, differentiation, and survival as well as cartilage development [61]. SKI is one of the SMAD-interacting proteins that negatively regulate the TGF-β signalling pathway during tissue homoeostasis [62]. CHSTT11 expression is required for proper chondroitin sulphate location, modulation of TGF- β and BMP signalling and cartilage development [63]. KIAA1217 deregulation has resulted in vertebral malformations in mammals [64]. Mutations in CREB3L2 induced skeletal deformities during zebrafish development through impacted notochord sheath cells [65], while CTC1 has been associated with the clinical feature of recurrent bone fractures [66].

In general, conserved methylation pattern between F1 and F3 generation have been reported in other ecotoxicological transgenerational studies [67,68]. These consistently differently methylated genes are potential markers for bone phenotype inheritance and bone health over generations in response to ancestral BaP exposure.

Differential transcription in the F1 and F3 vertebrae

The transcriptome data further confirmed the BaP-induced transgenerational dysregulation of biological processes involved in skeletogenesis, including osteoblast differentiation in the F1 (Figure 5a) and cartilage condensation and skeletal system development in the F3 (Figure 5C). Of the seven differently regulated genes involved in osteoblast differentiation (ACHE, LRP3, NPPC, PENK, SEMA4d, GDPD2 and TP53INP2) in the F1 vertebrae, six of them play specific roles in promoting osteoblast differentiation. The expression pattern in the bone tissue of the ancestrally exposed individuals suggested an inhibited osteoblastogenesis event during skeletogenesis possibly responsible for the significantly reduced osteoblast number. Both ACHE (down-regulated) and TP53INP2 (down-regulated) play a role in facilitating osteoblast differentiation and mineralization in vitro by activating Wnt/β-catenin signalling [69–71]. LRP3, an LRP family member, was found to be inhibited and possibly indicative for reduced osteogenesis [72]. C-type natriuretic peptide (NPPC, down-regulated) was found to stimulate osteoblast proliferation and collagen-X expression [73] and preproenkephalin (PENK, down-regulated) was reported to be involved in bone mineralization in mice [74]. Glycerophosphodiester phosphodiesterase domain containing 2 (GDPD2) was found to be inhibited after parental BaP exposure, which serve as an osteoblast differentiation-promoting factor in the mouse osteoblast-like cells [75]. Semaphorin 4D (SEMA4d, down-regulated) promotes bone resorption by inhibiting the function of osteoblasts [76].

In the F3 generation, a larger number of DEGs were observed compared to the F1 generation suggesting that ancestral BaP exposure had a more comprehensive impact on the unexposed offspring. The F3 transcriptome data are discussed at an extent in Jiezhang [9]. DEG-enriched GO biological processes and IPA canonical pathways were critical to skeletogenesis (e.g., cartilage condensation, calcium signalling, AMPK signalling, mTOR signalling, RANK signalling in osteoclasts, BMP signalling, AhR signalling, TGF-beta signalling; Figure 5c & d).

To determine if transcriptional changes were persistent between the F1 and F3 generation, a comparison was made between the DEGs of the F1 and F3 with over four times more DEGs identified in the F3 than in the F1 generation vertebrae but no conserved gene expression pattern, which is in line with other transgenerational studies. Nguyen et al. [77] reported that maternal intake of high-fat food during pregnancy of F0 induced more changes in the transcriptome of the F3 than the F1 generation, suggesting different mechanisms of inheritance in the F1 and F3 generation. Jing et al. [68], found that ancestral TCDD exposure of F0 pregnant rats caused a greater decrease in mRNA expression of the hepatic insulin-like growth factor 2 in the F3 than in the F1 offspring. Notwithstanding, the transcriptome of the F1 and F3 generation reveal the deregulation of the same biological processes, indicating the importance of functional pathways versus individual marker gene expression for tissue and organism-level phenotypes.

Integrated methylome and gene expression analysis reveals potential methylation driven transcriptome changes in the F3 vertebral deformities

DNA methylation has been proposed to be one of the pivotal epigenetic mechanisms through which environmental stressor exposure of one generation can lead to persistent phenotypic changes in the subsequent generations. The presented data revealed only a small number of DMGs, which were also differentially expressed per generation, which is consistent with recent studies showing a small degree of overlap between DEGs and DMGs due to environmental stressors including increased temperature in fish [78,79] and chemical exposure in fish [80].

It is suspected that BaP exposure-induced double-strand breaks are guided by the DNA methylation pattern and thus occur more often in methylated DNA sequences [81]. Direct BaP exposure induces dynamic, sequence-specific hypo- and hypermethylation events as demonstrated in human cell lines, affects DNA methyltransferase activity and the BaP metabolite benzo[a]pyrene-trans−7,8-dihydriol−9,10-epoxide has a binding preference for methylated DNA [82–84]. Recent research in mice has highlighted that BaP-induced modification may affect the methylation profile of imprinting genes in sperm, which was partly persistent upon the F2 generation [85]. There is mounting evidence that enhancement of the AHR mediated pathway is further promoting DNA hypomethylation, indicating that the AHR activation is modulating DNA methylation pattern in response to stress, which are proposed to be ‘fixed’ upon long-term exposure [86]. The here described differential methylation pattern and transcriptional regulation of the Bap-responsive pathways in the F1 and F3 corroborate the hypothesis that the AHR may regulate the DNA methylation over generations in a tissue-specific manner [87,88].

No conserved DMGs/DEGs were identified between the F1 and the F3 generation. The poor correlation between DNA methylation and gene expression indicates a potential role of DNA methylation on trans-acting elements that mediate gene transcription. DMRs could be located on regulatory elements that are outside of the gene regions used in this study (4 kb upstream to TES), or the inheritance of bone deformities upon ancestral BaP exposure on gene expression may be not due to DNA methylation but might be attributed to changes in histone modification or miRNA abundance. It has to be noted that the present methylome data originate from a preliminary screen with reduced replication.

In the F1 generation, the amino acid transporter 43A2 (SLC43A2), hypoxia-inducible factor (HIF3a), LRRN4 C-terminal-like protein (LRRNC4CL) and myelin basic protein (MBP) (Table 2) showed a correlation between de/methylation and expression. Transcription of the methionine transporter SLC43A2 has not been reported in osteoblasts, but methionine restriction affects osteoblast differentiation and bone mineralization in a male-biased manner [89]. Hypoxia inducible factors are a potential target for low bone mass treatment regulating osteoblastogenesis, although HIF3a expression has not been measured in bone tissue, HIF3 is thought to regulate HIF1/2 activity [90–92].

No function for LRRN4CL and MBP were specifically reported for bone tissue and osteoblast function and differentiation. The results indicate that DNA methylation is not the sole epigenetic factor responsible for the observed bone phenotype in the ancestrally BaP-exposed F1 generation.

In the F3 generation, 33 genes showed the expected correlation (Table 3), out of which secreted protein acidic and rich in cysteine (SPARC), chondroitin sulphate synthase 1 (CHSY1), and histone deacetylase 7 (HDAC7) were considered the most interesting for the presented inherited bone phenotype. SPARC is among the most abundant non-collagenous proteins expressed in mineralized tissues and is secreted by osteoblasts during bone formation [93]. SPARC is one of the marker genes for the osteoblast lineage [94] and is involved in osteoblast differentiation and mineral incorporation in the bone matrix [95]. In medaka, SPARC is present during the onset of bone development and regulates sclerotomal cell migration to the notochord [96]. The downregulation of SPARC may be responsible for the inhibited osteoblastogenesis and lead to the observed reduced osteoblast number in the tissue sections of BaP-derived F3 offspring, which is corroborated by reports that SPARC polymorphism with idiopathic osteoporosis in men [97]. CHSY1 is required for bone development and CHSY1 null mice exhibited a defect in chondrogenesis and delayed endochondral ossification [98]. Up-regulation of CHSY1 was found in the F3, suggesting a positive effect on chondrogenesis and hampering osteoblast differentiation. During bone resorption, HDAC7 plays a repressive role in osteoclast formation and is needed for osteoblast maturation [99]. In mice, conditional deletion of HDAC7 promoted osteoclast formation in vivo, associated with increased bone resorption and lower bone mass [100]. As a result, inhibited expression of HDAC7 might lead to an enhancement in osteoclast formation possibly entailing increased bone resorption. In addition, HDAC7 represses RUNX2 mediated histone acetylation, thus, impacting osteoblastic lineage commitment [101,102].

A literature review further revealed potential functional involvement of the identified differentially methylated and expressed genes in the ancestrally exposed F3 male fish. CD81, FAM171B, TJP1 and FREM2 have been demonstrated to be expressed in osteoblast and play a role in osteoblast differentiation and proliferation [103–107]. LDHB is a major regulator of osteoblast lactate metabolism, whose inhibition affects osteoblast differentiation in a sex-specific manner, while its expression is downregulated in osteoclasts during bone resorption [108–110]. CST3 is an osteoclast coupling factor regulating osteoblast differentiation [94,111]. ANK2and SRCIN are regulators of osteoblast and osteoclast differentiation, and deficiency of the former has been associated with reduced bone mass and calcification [112–114]. GNA13 deficiency results in osteoclast hyperactivation, while upregulation may protect from pathological bone loss [115]. NDUFS2 is part of the mitochondrial complex I assembly, which has been demonstrated to have an osteoclastogenic effect and is controlled by oestrogen receptors [116].

RNF220 suppresses GLI3 processing and enhances its degradation, thus affecting the SHH signalling while promoting the WNT/β-catenin signalling [117–119]. SHISA6 is a WNT/β-catenin signalling inhibitor, although its role in bone tissue or bone cells has not been described. GFAP, SLC38A4 and UQCC1 are regulation of bone development and growth [120] [107,121–123]. LLGL1 and NPDC1, are crucial factors for bone repair and potential candidate gene for fracture risk assessment in humans [124,125] and NDFIP1 is a regulator of bone development and repair through the BMP2/TGF-β pathway [126]. TTC28 and WNK1 expression is correlated with bone mineral density in mammals [127,128]).

Overall, the observed methylation and transcriptome changes, although not on the same genes, but in the same skeletogenesis pathways, indicate an ancestrally BaP exposure-induced transition from osteoblastogenesis towards chondrocyte differentiation, as these cells originate from the same sclerotome progenitors. The reduction of osteoblast numbers in both generations (reported here and in [8,9]) may result in a lack of osteoblastic mediator release to orchestrate osteoclastogenesis and bone resorption. Moreover, direct Bap exposure resulted in impaired stem cell self-renewal, cell replication and differentiation in multiple in vitro studies affecting cell growth and function [129–132]. Changes on the cell subpopulation levels may be responsible for modified gene expression, which are not driven by DNA methylation, and aberrant DNA methylation, which is not affecting gene expression in adult male fish [133], resulting in the here observed reduced alignment between gene expression and DNA methylation. As demonstrated by our group, the persistent impact of ancestral BaP exposure has developmental origins, thus, it is hypothesized that differential DNA methylation and changed gene expression are downstream compensation mechanisms in the adult tissue, which may be variable between generations. Moreover, the modification of DNA methylation and gene expression may have different origins, furthering the idea of the multidimensional nature of transcriptomes and methylomes in response to cellular and micro-environment [134]. Blake et al. [135] stated that conserved tissue-specific gene expression levels are more likely to be regulated through DNA methylation, which may account for same pathways found enriched in this study.

The data highlight the need for cell population- and differentiation stage-specific analyses to possible decipher cause–effect relationships between osteoblasts, osteoclast and chondrocytes. SPARC, CHSY1 and HDAC7 are highlighted for further investigation of this transgenerational bone phenotype.

Conclusion

The present study demonstrated the vertebral DNA methylation and concomitant gene expression changes in medaka fish with skeletal deformities after an ancestral BaP exposure. It is known that transgenerational effects can be regulated by different epigenetic mechanisms, including DNA methylation. However, genome-wide DNA methylation profiling and their relationship with corresponding gene expressions in the deformed bone were still lacking. The results suggest that the effects of ancestral BaP exposure on DNA methylation are distinct from gene expression changes, especially in the directly exposed generation (F1). The weak correlation between the degree of methylation changes and gene expression of overlapping DMR-associated genes and DEGs suggests that the regulation of gene expression is complex and DNA methylation alone is not enough to regulate the expressions of observed DEGs. However, DNA methylation in bone was still found to play an important role in BaP-induced skeletal deformities by mediating the bone-related physiological pathways and/or biological processes in a multigenerational and transgenerational manner. Additionally, transcriptional expressions of SPARC and HDAC7 were probably regulated by DNA methylation, which can serve as candidate markers for bone formation and resorption, respectively, during BaP-induced transgenerational skeletal deformities. Confirmation is needed if these DNA methylation changes originate from the germ lines.

Acknowledgements

The authors would like to acknowledge the generosity of Prof. Christoph Winkler of the National University of Singapore in providing us with the transgenic medaka fish line used in this study; Dr K. M. Cheung for his guidance on bioinformatics; Ms. K. L. Wong for help with fish dissection.

Disclosure statement

No potential conflict of interest was reported by the authors.

Data availability statement

Raw reads of the methylome have been deposited in NCBI’s BioProject accession number PRJNA845945. The RNA sequencing raw data of this study have been submitted to the NCBI Sequence Read Archive (SRA) (http://www.ncbi.nlm.nih.gov/sra) under the accession number PRJNA757563

Additional information

Funding

This work was supported by funding from the Southern Marine Science and Engineering Guangdong Laboratory (Guangzhou) (SMSEGL20SC02) and the National Natural Science Foundation of China (No. 41977371). Dr F.S. was supported by the National Institute of Environmental Health Sciences of the National Institutes of Health under award number 1R15ES032936-01.

Appendix I.

Selected top GO (BP) terms of conserved differentially methylated genes in the F1 and F3 vertebrate

Table

Biological process	LogP	p.adj	Gene
cell morphogenesis involved in differentiation	−8.928	.000	ACTB|GLI3|NFIB|NRCAM|PLXNA2|PRKCA|PTPRM|TESK1|LRP8|SYNGAP1|CYFIP1|TMEFF2|DISC1|CARMIL1|PARD3|SPTBN4|LINGO1|CCDC88C
regulation of nervous system development	−7.617	.000	GLI3|HOXB3|PLXNA2|SKI|SMARCD3|LRP8|SYNGAP1|CYFIP1|TENM4|DISC1|PARD3|HOOK3|LINGO1
cell junction organization	−6.644	.000	ACTB|MAPK14|GRIN2B|MYO9A|NRCAM|PRKCA|LRP8|SYNGAP1|CYFIP1|ADGRL3|DISC1|PARD3|CDH24|MDGA1|DOK7
modulation of chemical synaptic transmission	−5.114	.006	ABR|GRID1|GRIN2B|NPY2R|LRP8|SYNGAP1|RIMS2|CYFIP1|DISC1|MCTP1
gliogenesis	−4.603	.016	GLI3|NFIB|SKI|LRP8|TENM4|DISC1|PARD3|ADGRG6
cell surface receptor signalling pathway involved in cell-cell signalling	−4.283	.029	MAPK14|GLI3|GRIN2B|NPY2R|SKI|RIMS2|PTPRU|DISC1|RBMS3|RNF220|CCDC88C
chondrocyte differentiation	−4.124	.037	MAPK14|GLI3|NFIB|CHST11|CREB3L2
cartilage development	−3.897	.050	MAPK14|GLI3|HOXB3|NFIB|CHST11|CREB3L2
skeletal system development	−3.697	.066	MAPK14|GLI3|HOXB3|NFIB|SKI|CHST11|KIAA1217|CREB3L2|CTC1
connective tissue development	−3.244	.136	MAPK14|GLI3|HOXB3|NFIB|CHST11|CREB3L2
positive regulation of neuroblast proliferation	−4.055	.039	GLI3|SMARCD3|DISC1
cell junction assembly	−3.380	.113	ACTB|MYO9A|NRCAM|PRKCA|ADGRL3|PARD3|CDH24|MDGA1
behavior	−3.346	.113	GLI3|GRID1|GRIN2B|NPY2R|SYNGAP1|HIPK2|SPTBN4|VWA1|NEGR1

Appendix

Appendix II Conserved DMGs in the F1 and F3 vertebrae

Gene	F1V	F3V	Gene	F1V	F3V	Gene	F1V	F3V
	Promoter			Gene body			Gene body
A2ML1	-0.614	-0.383	ABCC4	-0.38	-0.381	CTC1	-0.319	-0.331
ABR	-0.498	-0.513	ABCC9	-0.323	-0.391	CTDSPL	0.326	0.43
ASF1B	0.501	0.352	ABHD4	-0.573	-0.501	CYFIP1	-0.351	-0.461
ATF3	0.345	0.309	ACTB	0.44	0.491	DISC1	0.312	0.328
C1orf109	0.371	0.521	ADARB2	-0.356	-0.528	DOK7	-0.618	-0.286
CLK3	-0.341	-0.477	ADCK1	-0.349	-0.35	DPP10	-0.36	0.349
CLSPN	0.602	0.341	ADGRG6	0.319	0.345	DYNC1H1	0.417	0.474
COL6A6	0.41	0.337	ADGRL3	-0.36	-0.29	EPPK1	-0.423	-0.3
ECHS1	0.32	0.336	AKR1D1	-0.418	-0.356	ERAP1	0.377	0.472
GOT2	0.483	0.489	ARID2	-0.312	-0.563	GLI3	0.357	0.403
IQGAP3	-0.41	-0.37	ARRDC1	-0.478	-0.492	GNG4	0.353	0.325
MPV17	-0.424	-0.447	ASAP1	-0.537	-0.337	GRID1	-0.438	-0.289
MRC1	-0.486	-0.371	B4GALT1	-0.39	-0.51	GRIN2B	0.343	0.596
NPY2R	0.328	0.323	CARMIL1	-0.428	-0.33	H1-6	0.516	0.474
NRCAM	0.345	0.422	CCDC88C	-0.344	-0.598	HIPK2	0.308	0.347
PTGFRN	0.332	0.425	CDH24	0.468	0.389	HOOK3	-0.399	-0.367
PZP	-0.614	-0.383	CHD1L	-0.538	-0.44	HOXB3	-0.361	-0.336
RNF220	0.301	0.406	CHST11	-0.735	-0.34	HSPA14	-0.448	-0.335
SMYD1	-0.55	-0.505	CPNE5	0.292	0.645	IPMK	0.305	0.276
TESK1	-0.311	-0.343	CREB3L2	-0.275	-0.438	KCNIP3	0.413	0.354
TULP4	0.342	0.345	CSMD3	-0.522	-0.344	KIAA1217	-0.408	-0.449
LINGO1	0.491	0.823	PDE5A	0.62	0.544	SLC2A13	0.341	0.404
LRP8	-0.321	-0.396	PDK4	0.385	0.426	SLC41A1	-0.378	-0.289
LSS	0.319	0.363	PLEKHA5	-0.409	-0.417	SLC4A1	0.47	0.378
MAPK14	-0.378	-0.334	PLEKHM2	-0.534	-0.428	SLC9A7	-0.42	-0.443
MARCHF9	-0.312	-0.32	PLXNA2	-0.362	-0.388	SMARCD3	-0.332	-0.348
MCTP1	0.362	0.374	PPT2	0.581	0.505	SPTBN4	-0.336	-0.281
MDGA1	0.73	0.679	PRKCA	-0.377	-0.287	SVEP1	0.345	0.423
MLPH	0.45	0.347	PTPRM	-0.403	-0.363	SYNGAP1	0.315	0.364
MTG1	-0.32	-0.336	PTPRU	-0.498	-0.342	TAX1BP1	-0.658	-0.341
MYO9A	-0.44	-0.396	RASSF5	0.325	0.378	TBC1D10B	0.691	0.39
MYT1	-0.507	-0.303	RBFOX3	-0.305	-0.317	TENM4	-0.313	-0.3
NCAM2	-0.426	-0.326	RBMS3	-0.437	-0.353	THSD7B	-0.448	-0.355
NCOA3	-0.476	-0.348	RFTN1	0.088	0.364	TMEFF2	0.466	0.355
NEGR1	-0.411	-0.337	RIMS2	-0.384	-0.339	VPS13B	-0.409	-0.34
NELL1	0.299	0.411	RNF31	-0.448	-0.301	VWA1	0.35	0.39
NFIB	-0.352	-0.335	SATB1	0.304	0.323	WDR55	-0.38	-0.526
OSBPL10	0.303	0.316	SH3D19	-0.642	-0.402	ZDHHC14	0.383	0.37
PARD3	0.338	0.631	SKI	-0.291	-0.294	ZNF609	-0.345	-0.38