Tissue and sex-specific programming of DNA methylation by perinatal lead exposure: implications for environmental epigenetics studies

ABSTRACT

Early developmental environment can influence long-term health through reprogramming of the epigenome. Human environmental epigenetics studies rely on surrogate tissues, such as blood, to assess the effects of environment on disease-relevant but inaccessible target tissues. However, the extent to which environment-induced epigenetic changes are conserved between these tissues is unclear. A better understanding of this conservation is imperative for effective design and interpretation of human environmental epigenetics studies. The Toxicant Exposures and Responses by Genomic and Epigenomic Regulators of Transcription (TaRGET II) consortium was established by the National Institute of Environmental Health Sciences to address the utility of surrogate tissues as proxies for toxicant-induced epigenetic changes in target tissues. We and others have recently reported that perinatal exposure to lead (Pb) is associated with adverse metabolic outcomes. Here, we investigated the sex-specific effects of perinatal exposure to a human environmentally relevant level of Pb on DNA methylation in paired liver and blood samples from adult mice using enhanced reduced-representation bisulphite sequencing. Although Pb exposure ceased at 3 weeks of age, we observed thousands of sex-specific differentially methylated cytosines in the blood and liver of Pb-exposed animals at 5 months of age, including 44 genomically imprinted loci. We observed significant tissue overlap in the genes mapping to differentially methylated cytosines. A small but significant subset of Pb-altered genes exhibit basal sex differences in gene expression in the mouse liver. Collectively, these data identify potential molecular targets for Pb-induced metabolic diseases, and inform the design of more robust human environmental epigenomics studies.

KEYWORDS:

• Toxicoepigenetics
• DNA methylation
• Developmental origins of health and disease (DOHaD)
• surrogate tissue
• sex differences
• environmental epigenomics

Introduction

It is increasingly evident that environmental exposures and nutritional perturbations during critical windows of development can influence the long-term risk of disease. One important mechanism by which early-life environment can reprogram disease risk is via modification of the epigenome. Normal human development is governed by precise spatio-temporal regulation of gene expression, and this is driven in large part by epigenetic mechanisms. Methylation of DNA at the 5th-position of cytosines (5-methylcytosine, 5mC) occurring in cytosine-phosphate-guanine (CpG) dinucleotides is critical for establishment of tissue-specific gene expression patterns, regulation of imprinted genes, maintenance of genome stability, and silencing of transposable elements[1]. 5mC undergoes distinct waves of reprogramming during early embryonic development and in primordial germ cells, and is stable and mitotically heritable [2]. Epigenetic reprogramming during these periods of rapid growth and differentiation can be altered by environmental or nutritional perturbation. 5mC can undergo further oxidation to 5-hydroxymethylcytosine (5hmC) [3]. Although 5hmC is an intermediate in the DNA demethylation pathway, recent studies suggest that the genomic distribution, downstream targets, and functional roles of 5hmC are distinct from that of 5mC [3]. Environmental exposures may also lead to reprogramming of 5hmC, and these changes are associated with altered transcription and adverse health outcomes [4–6].

Exposure to heavy metals, including lead (Pb), remains a significant public health concern, particularly in poor urban communities due to various sources of contamination (e.g., leaded paint and pipes, mining, lead-glazed ceramics). Although the neurological effects of early-life Pb exposure have been well-documented, substantially less is known about how Pb exposure impacts the liver and other metabolic organs. Notably, recent human studies suggest that blood Pb levels are associated with liver disease, metabolic syndrome, and diabetes [7–10]. Likewise, maternal blood Pb levels are associated with increased DNA methylation of the imprinted MEG3 locus in offspring blood, concomitant with rapid gains in body fat [11]. In animal models, we and others have recently reported that perinatal [12–15] and post-natal [16–18] exposure to Pb leads to sex-specific changes in body weight, metabolic function, and gut microbiota. Post-natal Pb exposure leads to hypermethylation of genes associated with glucose and lipid metabolism in the rat liver [16]. However, the effects of Pb exposure during gestation and lactation on hepatic DNA methylation, and the implications this may have for long-term disease risk, are unclear.

In human population-based environmental epigenetics studies, easily obtainable sources of DNA (e.g. hair, blood, stool, saliva) are used as surrogates for the tissue(s) targeted by the environmental exposure. It is unknown, however, to what extent exposure-induced epigenetic changes in surrogate tissues reflect those occurring in target tissues, and whether this overlap differs by sex. Likewise, whether toxicant-induced epigenomic changes persist in both target and surrogate tissues over time within individuals is unclear. To address these important questions, the National Institute of Environmental Health Sciences (NIEHS) established the TaRGET II consortium [19]. Using mouse models of human-relevant environmental exposures, the objective of TaRGET II is to assess the conservation of environment-induced epigenomic perturbations between easily obtainable surrogate and disease-relevant, but inaccessible, target tissues. As part of this effort, we used an established mouse model of perinatal Pb exposure to determine the effects of Pb on genome-wide DNA methylation in the liver and blood of 5-month old, adult animals. We hypothesized that perinatal Pb exposure would lead to sex-specific alterations in DNA methylation in liver at genes important for metabolic homoeostasis. We further hypothesized that a subset of differentially methylated sites would overlap between liver and blood.

Materials and methods

Animal exposure paradigm

Mice utilized for breeding and exposure were obtained from a colony maintained for over 230 generations with the A^vy allele passed through the male line, resulting in forced heterozygosity on a genetically invariant background with 93% identity to C57BL/6 J [20,21]. Virgin a/a females (6–8 wks old) were mated with virgin a/a males (7–9 wks old), and randomly assigned to control or Pb exposure via drinking water. For Pb exposure, Pb-acetate was mixed with water to result in a Pb concentration of 32 ppm in drinking water, to model human-relevant maternal exposure in the 16–60 μg/dL range [12]. Exposure water was made by dissolving Pb (II) acetate trihydrate (Sigma‐Aldrich) in a single batch of distilled water, and Pb concentrations were verified using inductively coupled plasma mass spectrometry with a limit of detection of 1.0 µg/L (ICPMS; NSF International). Animals were maintained on a phytoestrogen‐free modified AIN‐93 G diet (TD.95092, 7% Corn Oil Diet, Harlan Teklad). The a/a dams in the Pb group were exposed to Pb‐supplemented drinking water for 2 weeks prior to mating a/a males. Exposure was continued during gestation and lactation. After weaning, the resulting pups were weighed and switched to Pb‐free drinking water (Figure 1). A subset of pups from each exposure group, representing approximately 1–2 male and 1–2 female offspring per litter, was maintained to 5 months of age (n = 6 males and 6 females per treatment). All animals had access to food and drinking water ad libitum throughout the experiment while housed in polycarbonate‐free cages. The study protocol was approved by the University of Michigan Institutional Animal Care and Use Committee (IACUC).

Figure 1. Mouse treatment paradigm showing timing of Pb exposure and animal sacrifice/tissue collection. Dams were exposed to normal chow and either control or Pb‐acetate water, starting two weeks prior to mating, through gestation, until offspring were weaned on postnatal day 21. At this point, Pb exposure ceased and offspring were administered control water and chow. At 5 months of age, matched liver and blood tissues were collected for DNA methylation analysis. Pb doses and average blood lead levels are shown

Body composition and tissue collection

Body weights were measured for each mouse on a weekly basis (Mettler Toledo). Total body fat and lean mass were measured via EchoMRI (Houston, TX USA). Body fat percentage and lean mass percentage were calculated as a percent of total body weight. Upon euthanasia at 5 months of age, mesenteric fat, blood, and liver samples were collected in accordance with protocols established by the TaRGET II Consortium (Figure 1 and [19]). Briefly, mice were fasted for 6 h prior to euthanasia. Euthanasia was carried out via CO₂ asphyxiation and bilateral pneumothorax. Blood was collected by cardiac puncture, followed by whole body perfusion with cell culture grade 0.9% saline solution (Sigma Life Sciences). Next, liver and mesenteric fat were collected and weighed. Mesenteric fat was collected by dissecting white adipose tissue away from the stomach, spleen, pancreas, and intestines. Relative liver and mesenteric fat weights were expressed as a percent of total body weight. Liver was separated into individual lobes, and the left lobe was used for the analyses in this study. Liver and fat samples were immediately snap-frozen in liquid nitrogen and stored at −80° C. Blood was collected with EDTA anticoagulant, centrifuged, and plasma was removed. Red blood cells were lysed using Erythrocyte Lysis Solution (Qiagen, 79,217). The remaining white blood cell fraction was washed with PBS, pelleted by centrifugation, and resuspended in Buffer RLT (Qiagen, cat # 1,053,393). Crude cell lysate was stored at −80 C until DNA extraction.

DNA extraction and enhanced reduced representation bisulphite sequencing

The left lobe of the liver was cryo-pulverized and suspended in Buffer RLT (Qiagen, cat # 1,053,393). Blood and liver DNA extraction were performed using the AllPrep DNA/RNA/miRNA Universal Kit (Qiagen #80,224). Enhanced reduced representation bisulphite (ERRBS) was performed at the University of Michigan Epigenomics and Advanced Genomics Cores as described previously [22,23]. Briefly, 50 ng of genomic DNA was digested using MspI, a restriction enzyme that preferentially cuts CG-rich sites. The digested DNA was then purified using phenol:chloroform extraction and ethanol precipitation in the presence of glycogen, before blunt-ending and phosphorylation. A single adenine nucleotide was next added to the 3′ end of the fragments in preparation for the ligation of the adapter duplex with a thymine overhang. The ligated fragments were cleaned, then processed for size selection on agarose gel. Selected fragments were treated with sodium bisulphite to convert unmethylated cytosines to uracils, which are then replaced with thymines during PCR amplification. Both methylated and hydroxymethylated CpGs are protected from bisulphite conversion. Therefore, ERRBS detects a combination of 5mC and 5hmC. Final libraries were next cleaned up with AMPure XP beads (Product #A63880; Beckman Coulter), quantified using the Agilent TapeStation genomic DNA kit (Catalogue #G2991AA; Agilent) and Qubit High Sensitivity dsDNA (Catalogue #Q32850; Invitrogen). Single-end, 50 bp reads were obtained for each library by sequencing on the HiSeq 4000 system (Illumina). ERRBS samples were multiplexed, with three samples per sequencing lane.

Bioinformatics pipeline, quality control, and differential methylation analysis

Snakemake [24] was used to manage the bioinformatics workflow, and Anaconda (Anaconda Software Distribution) was used to manage software dependencies and create reproducible software environments. For quality control, FastQC (v0.11.3) was employed to assess the overall quality of each sequenced sample and identify specific reads and regions for trimming. TrimGalore (v0.4.5) was used to trim low-quality bases (quality score lower than 20), adapter sequences (required overlap of 6bp), and end-repair bases from the 3ʹ end of reads. Reads shorter than 20bp after trimming were discarded. For alignment and methylation calling we used Bismark [25] (v0.19.0), an integrated alignment and methylation call program that performs unbiased alignment (by converting residual cytosines to thymines prior to alignment in both reads and reference). Within Bismark, reads were aligned to the reference genome (mm10) using Bowtie2 [26] (v2.3.4) with default parameter settings (multi-seed length of 20bp with 0 mismatches). Methylation calls were reported for all nucleotides with a read depth of at least 5. We used the methylSig R package [27] (v0.5.0) to identify differentially methylated CpG positions and/or regions. Differential methylation for each comparison was tested using methylSigDSS(), which tests for differential methylation under general experimental design using a beta-binomial approach with the ‘arcsine’ link function [28]. To control for batch effects, run was included as a covariate in the model. After obtaining p-values, we adjusted for multiple testing using the FDR approach. The results of the differential methylation tests were then annotated using the annotatr R Bioconductor package [29] (v1.5.9). The annotate_regions function was used to generate genomic annotations, which include CpG annotations [CpG islands (CGI), shores, shelves, open sea (InterCGI)], genic annotations (exon, intron, promoter, 5ʹUTR, 3ʹUTR), and gene IDs. CpGs with read coverage > 1000 were removed because they were likely the result of PCR amplification, and CpGs with read coverage < 10 were removed due to decreased power to detect differential methylation. Opposite strand CpGs at the same position were combined via destranding. We performed differential methylation testing on individual CpG sites (DMCs), requiring sufficient sequencing coverage for a minimum of 4 samples from the Pb group, and 4 samples from the control group for a site to be tested. Differentially methylated regions (DMRs) were identified in 50 bp tiles using the same requirements. For both DMCs and DMRs, sites with FDR <0.05 and an absolute difference in methylation of > 10% were considered significant.

Identification of DMRs within sex-specific mouse liver enhancers

ChromHMM tracks from a 2013 study by Sugathan et al [30]. were kindly provided by Dr. David Waxman. Bedtools version 2.28.0 was used to identify the intersections between male or female DMRs and their respective sex-specific mouse liver enhancer elements. Specifically, intersections were generated between DMRs and 5 chromatin states that were identified as enhancers in this study based on a combination of DNAse hypersensitivity and the presence of enhancer-specific chromatin marks (states 5, 6, 9, 10, 11) [30].

Pathway analysis of DMRs

Poly-Enrich [31] was used to assess the biological pathways enriched among the DMRs. Separate pathway analyses were conducted for each sex, tissue, and direction of differential methylation, and were restricted to proximal promoter regions (within 1 kb of transcription start sites). The following Gene Ontology gene sets used for this analysis: Biological Process, Cellular Component, and Molecular Function. Pathways with FDR <0.05 were considered statistically significant.

Analysis of overlap in differentially methylated sites between blood and liver

Annotated lists of DMCs and DMRs for males and females were compared between blood and liver to identify chromosomal locations that directly overlapped between the two tissues. Overlapping sites with changes in methylation in the same direction in both tissues were considered biologically relevant. In a further analysis, the full lists of DMR-associated genes in blood were compared to those in liver to identify a list of genes in common between the two tissues. Separate lists were generated for males and females. These gene lists were analysed using DAVID version 6.8, Laboratory of Human Retrovirology and Immunoinformatics (LHRI),

(https://david.ncifcrf.gov) to identify pathways that were significantly enriched (FDR <0.05).

qRT-PCR

Quantitative real time reverse transcription PCR (qRT-PCR) was carried out using RNA samples from male and female livers. RNA was converted to cDNA using iScript cDNA Synthesis Kits (Bio-Rad cat# 1,708,891) following manufacturers’ instructions. Quantitative PCR reactions were set up for Gnas and 3 housekeeping genes (Actb, Hmbs, Psmc4) using the SsoAdvanced Universal SYBR Green Supermix (Bio-Rad) via the manufacturer’s instructions. Gnas primer sequences are shown in Supplemental Table S1. Primers for Actb, Hmbs, and Psmc4 were PrimePCR primers purchased from Bio-Rad. Samples were run in triplicate, and controls included no template and genomic DNA conditions. Quantitative real-time PCR was conducted using the CFX384 Real-Time PCR Detection System (Bio-Rad). The cycling protocol was as follows: 1) 2 mins at 95 C x 1 cycle, 2) 5 seconds at 95 C and 30 seconds at 60 C x 40 cycles, and 3) melt curve with 0.5 C 5 second increments from 65 C – 95 C. Gnas expression was calculated relative to the geometric mean of the Ct values for the housekeeping genes using the delta Ct method.

RNA-seq gene expression

RNA was isolated from the left lobe of PND21 or 5 month mouse liver using the Qiagen AllPrep DNA/RNA/miRNA Universal Kit. Library preparation and sequencing were conducted on postnatal day 21 and 5- month liver and blood at the University of Michigan Advanced Genomics Core. Liver and blood samples from mice that had gross liver tumours (1 female Control and 1 female Pb) were excluded. RNA Integrity Numbers (RIN) were assessed by Agilent 2200 TapeStation analysis. Electropherograms indicated that the RINs were greater than 6 for most liver RNA samples greater than 8 for most blood samples. Library preparation was conducted using the Illumina TruSeq stranded mRNA Library Prep Kit (5- month samples) or the KAPA mRNA Hyper Prep Kit with Dual Indexing Adapters (PND21 samples) according to manufacturer instructions. Quantity and quality of the prepared libraries were confirmed with the Agilent 2200 TapeStation. Sequencing (paired-end 50 base pair reads) was carried out on the Illumina HiSeq 4000 (5 month samples) and the NovaSeq S2 flow cell (PND21). Data quality control and analysis were conducted by the TaRGET II Consortium Data Coordinating Centre at Washington University in St. Louis, MO. Data were normalized using the RUVr method of the RUVseq R package [32]. Genes were first selected with more than 1 cpm in more than 10 RNA-seq samples. The resulting gene signal matrix was divided into adult liver, PND21 liver, adult blood and PND21 blood. RUVr with k = 3 was run separately for these four groups and samples were stratified base on sex and treatment. Additional methods for data processing can be found at https://targetepigenomics.org/documents/. Read count data for individual mice were plotted, and gene expression differences were analysed for statistical significance using linear mixed effects regression with the lme4 package and R version 3.5.1.

Results

Litter parameters and phenotype

Perinatal lead exposure did not significantly alter the offspring sex ratio, litter size, or mortality rates (p = 1.0, 0.9, and 0.2, respectively, Fisher’s Exact Test) (Table 1). We have previously reported that perinatal Pb exposure leads to sex-specific changes in food intake, body weight, and body fat later in life (n = 12–18 mice per sex) [13]. In this study, we observed no significant differences in body weight, total body fat, or body fat percentages in either males or females, an observation that may have been due to a smaller sample size (n = 6 mice per sex) (Supplemental Figure 1 A-F) than our previous Pb studies observing significant metabolic effects [12,13]. Nevertheless, consistent with our previous observations, we observed a trend towards increased total and relative mesenteric fat weights in males at 5 months, although this did not reach statistical significance (p = 0.06, Supplemental Figure 1 G-J). Total and relative liver weights did not differ significantly between exposed and control groups in either males or females (Supplemental Figure 1 K-N).

Table 1. Litter parameters

Exposure Group	# Litters	# Pups Born	# Pups Died	Mean Pups per Litter (± SD)^a	Pup Mortality Rate (%)^b	% Females^a
Control	6	53	3	7.57 ± 3.36	5.7	46%
Pb	7	52	7	7.43 ± 1.51	13.5	47%

Statistical analyses (Pb vs. Control) were conducted using Fisher’s Exact Test.

^ap = ≥ 0.9

^bp = 0.2

Epigenome-wide differential DNA methylation with Pb exposure

To assess the effects of perinatal Pb exposure on DNA methylation across target and surrogate tissues, ERRBS was conducted on liver and blood samples from 5-month-old Pb and control-exposed offspring. Tables 2 and 3 illustrate the total number of regions tested and the number of differentially methylated cytosines (DMCs) and DMRs identified in each exposure group. Although Pb exposure ceased at 3 weeks of age, this analysis revealed thousands of sex-specific DMCs and DMRs in the adult liver and blood (Table 2–3 and Supplemental Tables 2–9). The relative proportion of hypo- and hyper-methylated regions/cytosines was similar for each condition (Table 2–3). Although observed methylation changes were between 10 and 20% at most sites with Pb exposure, numerous sites exhibited changes in methylation that were much greater, up to 70% in females and 74% in males (Figure 2(a-d) and Supplemental Figure 2). We annotated the DMCs and DMRs to the mouse mm10 genome to determine the genomic locations of DMCs and DMRs relative to all sites tested, stratified by sex, tissue, and direction of differential DNA methylation. Compared to all regions tested, DMCs/DMRs were significantly depleted from CpG islands, promoters, 5ʹUTRs, and exons, and were significantly enriched for regions greater than 4000 bp from a CpG island (interCGI) and intergenic regions (chi squared test, p <.05) (Supplemental Figure 3 A-H). This pattern was similar across sex and tissue (Supplemental Figure 3 A-H).

Table 2. Differentially methylated cytosines (DMCs) in female and male blood and liver identified using ERRBS

Condition	Total	# Hypermethylated (% Total)	# Hypomethylated (% Total)	Total Tested
Female Blood	1,628	747 (45.9)	881 (54.1)	1,356,386
Female Liver	2,404	1,350 (56.2)	1,054 (43.8)	1,228,997
Male Blood	2,175	1,076 (49.5)	1,099 (50.5)	1,362,718
Male Liver	3,808	1,598 (42.0)	2,210 (58.0)	1,173,995

Methylation changes are expressed as lead vs. control.

Table 3. Differentially methylated regions (DMRs) in female and male blood and liver identified using ERRBS

Condition	Total	# Hypermethylated (% Total)	# Hypomethylated (% Total)	Total Tested
Female Blood	1,035	551 (53.2)	484 (46.8)	646,876
Female Liver	1,740	792 (45.5)	948 (54.5)	597,077
Male Blood	1,404	724 (51.6)	680 (48.4)	650,782
Male Liver	3,342	1,874 (56.1)	1,468 (43.9)	563,681

Methylation changes are expressed as lead vs. control.

Figure 2. Volcano plots showing differentially methylated regions (DMRs) for lead vs. control in female blood (a), female liver (b), male blood (c), and male liver (d). Blue: regions significantly hypomethylated with Pb exposure. Green: regions significantly hypermethylated with Pb exposure

The liver exhibits marked sexual dimorphism, with distinct patterns of chromatin accessibility and enhancer modifications driving sex-biased gene expression [30]. Given the sex specificity of differential methylation, we investigated whether DMRs occurred at sex-specific enhancers in the mouse liver. To this end, we stratified DMRs by sex, tissue, and direction of DNA methylation, and determined the overlap with previously published, sex-specific mouse liver enhancers [30]. Specifically, we determined whether DMRs in females fell within female-specific enhancers, and whether DMRs in males occurred within male-specific enhancers. Five chromatin states were identified as enhancers using ChromHMM, based on a combination of DNAse hypersensitivity and specific chromatin marks. These include the presence of enhancer chromatin marks H3K27 acetylation and/or H3K4 mono-methylation, and a lack of the promoter-enriched H3K4 trimethylation mark. (see Figure 3 and reference [30]). Notably, in both sexes and tissues we observed a subset of DMRs falling within their respective sex-specific enhancer regions. Chromatin states 5 and 6 are characterized by robust DNAse hypersensitivity [30], suggesting an open enhancer state [33]. In contrast, states 9, 10 and 11 exhibit low DNAse hypersensitivity [30], suggesting closed chromatin and an inactive enhancer [33]. In all conditions, the majority of DMRs were found in closed enhancer regions (Figure 3(a-d)). The relative proportions of different enhancer states overlapping DMRs were similar for both liver and blood, but differed substantially between sexes (Figure 3(a-d)). Compared to females, in male blood and liver a greater proportion of DMRs were found in open enhancer regions (Figure 3(c-d)).

Figure 3. Pie charts depicting the overlap between DMRs in each sex and their respective sex-specific enhancer states in female blood (a), female liver (b), male blood (c) and male liver (d). Totals are the number of DMR-enhancer overlaps for each condition. Published enhancers are from Sugathan and Waxman, Mol. Cell. Biol. 2013, and chromatin states are defined according to the degree of enrichment with H3K4me1, H3K27Ac, and DNAse hypersensitivity (DHS): State 5 – low K4me1, medium K27ac, high DHS; State 6 – high K4me1, high K27ac, high DHS; State 9 – high K4me1, high K27ac, low DHS; State 10 – low K4me1, medium K27ac, low DHS; State 11 – high K4me1, low K27ac, low DHS

Pathway analysis of DMRs

Having established that early life Pb exposure leads to sex-specific changes in DNA methylation in adulthood, we next determined the cellular pathways enriched among the DMRs (Supplemental Tables S10-S13). We stratified our analyses by sex, tissue, and direction of differential methylation, and focused on proximal promoter regions (within1 kb of transcription start sites). Pathways with a FDR<0.05 were considered significant. This analysis revealed that the pathways enriched among hypo-methylated sites were distinct from those enriched among hyper-methylated sites. In addition, the enriched pathways differed across sex and tissue. In spite of these differences, overall, we observed that Pb-altered methylation of CpG sites in key organ developmental pathways including mammary gland, bone, nervous system, and kidney. Interestingly, we observed enrichment of pathways involved in regulation of neurotransmitter levels and activity in female and male liver (Supplemental Tables S10-S13).

Perinatal Pb exposure and regulation of imprinted loci

Previous findings from our lab and others demonstrated that developmental exposure to BPA and Pb led to altered CpG methylation at imprinted loci [4,11,34]. Therefore, we investigated whether perinatal Pb exposure disrupted methylation at imprinted loci in liver and blood, and whether these changes overlapped between the two tissues. To do this, we compared the lists of Pb-induced DMCs/DMRs for each condition to a database of mouse-imprinted genes [35]. Of 139 imprinted genes interrogated, 112 were covered in this study using ERRBS (Supplemental Table S14). We observed Pb-induced DMRs/DMCs at 44 imprinted genes combined across tissue and sex. Specifically, we identified 3 in female blood, 18 in female liver, 22 in male blood, and 21 in male liver. The enrichment of imprinted genes among the full set of DMC-associated genes (both tissues and sexes) was statistically significant (p < 0.01, hypergeometric test). See Tables 4–5 for DMCs mapping to imprinted genes in females and males, respectively, and Supplemental Tables S15-S18 for DMRs. In liver samples, we discovered 5 imprinted genes with DMCs in both males and females, including Commd1, Gnas, Nespas, Pde10a, and Pde4d (Tables 4–5). Methylation at CpG sites in Gnas/Nespas was in the same direction in both males and females (Figure 4(a)). In blood, CpG sites in Meg3 were hypermethylated in both males and females (Figure 4(b)). We next compared DMC-associated imprinted genes between blood and liver within each sex. In females, no imprinted genes overlapped between blood and liver. In males, we observed DMCs at Begain, Cdkn1c, Commd1, Peg12, Rasgrf1, and Snrpn in both blood and liver, and found that CpGs in Begain, Peg12, Rasgrf1, and Snrpn exhibited changes in methylation in the same direction in both blood and liver (Table 5 and Figure 4(c-f)). We previously observed that BPA exposure led to changes in methylation and hydroxymethylation at the Gnas locus, concomitant with increased Gnas transcript expression [4]. Therefore, we investigated whether Pb-induced changes in methylation at the Gnas locus were accompanied by changes in Gnas transcript expression using qRT-PCR. We observed a significant reduction in Gnas expression with Pb exposure in female mice, and a trend towards reduced expression in male mice (Figure 5). Collectively, these data suggest that Pb exposure alters methylation at imprinted loci in a manner that is largely tissue and sex-dependent.

Table 4. Differentially methylated cytosines (DMCs) identified at imprinted loci in female blood and liver

DMCs in Blood
Gene	Chromosome	Position	Location	Change in methylation	FDR q-value
Meg3	12	109,547,833	Intron/1 to 5 kb	36.92	0.002
Slc22a3	17	12,454,881	Intron	−14.18	6.62E-05
Tnfrsf23	7	143,670,844	Intron/Exon/3ʹUTR	−10.12	0.002
DMCs in Liver
Gene	Chromosome	Position	Location	Change in methylation	FDR q-value
Art5	7	102,100,155	Promoter/Intron/Exon/5ʹUTR	15.81	0.021
Axl	7	25,764,245	Promoter/Exon/Intron	−33.18	0.048
Axl	7	25,781,580	Intron	24.81	0.039
Commd1	11	22,972,256	Intron	−36.26	0.007
Commd1	11	22,972,731	Intron/Exon	16.37	0.002
Fthl17e	X	9,028,755	1 to 5 kb	−14.42	0.0009
Gatm	2	122,595,535	Exon	−23.74	0.014
Gnas/Nespas	2	174,285,424	Exon	13.02	0.001
Gnas	2	174,285,427	Exon	10.32	0.024
Grb10	11	11,948,460	Intron	−19.24	0.040
H19	7	142,578,975	Promoter/1 to 5 kb	−23.47	0.033
Kcnq1	7	143,295,410	Intron	−13.35	0.001
Kcnq1/Kcnq1ot1	7	143,295,427	Intron/Exon	−16.02	0.021
Kcnq1 Kcnq1ot1	7	143,295,437	Intron/Exon	−14.64	0.0001
Kcnq1/Kcnq1ot1	7	143,295,442	Intron/Exon	−21.31	0.005
Kcnq1	7	143,395,334	Intron	22.05	0.0003
Magel2	7	62,378,976	Exon	−20.11	0.023
Pde4d	13	109,659,183	Intron	−18.56	0.039
Pde10a	17	8,754,083	Intron	−14.17	0.003
Zrsr1	11	22,972,256	Promoter/Exon/5ʹUTR	−36.26	0.007
Zrsr1	11	22,972,731	Promoter/Exon/5ʹUTR	16.37	0.002

Methylation changes are expressed as lead vs. control.

Table 5. Differentially methylated cytosines (DMCs) identified at imprinted loci in male blood and liver

DMCs in Blood
Gene	Chromosome	Position	Location	Change in methylation	FDR
Airn	17	12,822,718	Intron	−11.17	0.0007
Art5	7	102,104,325	1 to 5 kb	23.58	0.029
Axl	7	25,774,103	Intron/Exon	−10.66	0.012
Begain	12	109,056,064	Intron/1 to 5 kb	10.21	0.044
Cdkn1c	7	143,461,132	Promoter/1 to 5 kb	17.77	0.034
Commd1	11	22,973,934	Intron/Exon/1 to 5 kb	19.41	0.012
Grb10	11	12,025,803	Intron	−10.20	0.047
Grb10	11	12,038,672	1 to 5kb	13.65	0.026
H13	2	152,707,734	Exon/3ʹUTR	−10.54	0.0005
Igf2r	17	12,726,414	Intron	17.13	0.008
Meg3	12	109,540,296	Promoter/1 to 5 kb	38.50	0.0009
Mkrn3	7	62,419,838	Exon	27.44	0.030
Peg12	7	62,463,945	Exon	−24.76	0.042
Plagl1	10	13,117,338	Intron	−23.32	0.006
Rasgrf1	9	89,958,385	Intron	19.19	6.68E-05
Snrpn	7	60,012,762	Intron	21.67	0.0002
Usp29	7	6,770,496	Intron	−12.12	0.003
Zim1	7	6,678,078	Exon	−12.05	0.0004
Zrsr1	11	22,973,934	Exon	19.41	0.012
DMCs in Liver
Gene	Chromosome	Position	Location	Change in methylation	FDR
Begain	12	109,055,956	Intron/1 to 5 kb	23.07	0.012
Cdkn1c	7	143,459,490	Promoter/Intron	−49.58	0.004
Cdkn1c	7	143,460,320	Promoter/1 to 5 kb/Exon	−13.94	0.0033
Commd1	11	22,978,331	Intron/Exon	−36.13	0.001
Dlk1	12	109,449,906	1 to 5 kb	10.67	0.031
Fbxo40	16	36,990,334	Exon/5ʹUTR	−10.96	0.021
Gnas/Nespas	2	174,298,185	Intron/Exon/1 to 5 kb	15.82	0.048
Gnas/Nespas	2	174,299,492	Intron/Exon/1 to 5 kb	14.23	0.006
Htr2a	14	74,706,653	Exon/3ʹUTR	17.33	0.029
Klf14	6	30,958,628	Exon	−14.23	0.014
Pde4d	13	108,848,465	Intron	28.82	0.025
Pde4d	13	109,936,716	Exon	12.18	0.030
Pde10a	17	8,547,124	Intron	11.66	0.037
Pde10a	17	8,927,375	Intron	−12.70	0.004
Peg12	7	62,463,837	Exon	−20.88	0.007
Rasgrf1	9	89,953,923	Exon	21.35	8.36E-36
Rhox5	X	37,807,531	Intron	22.33	0.006
Sfmbt2	2	10,368,329	1 to 5 kb	−28.39	0.004
Snrpn/Snurf	7	60,005,199	Intron/Promoter	11.12	0.040
Th	7	142,895,420	Intron/Exon/3ʹUTR	14.65	0.002

Methylation changes are expressed as lead vs. control.

Figure 4. Lollipop diagrams showing differentially methylated cytosines (DMCs) with Pb exposure in: (a) Gnas in male and female liver, (b) Meg3 in male and female blood, (c) Begain in male blood and liver, (d) Peg12 in male blood and liver, (e) Rasgrf1 in male blood and liver, and (f) Snrpn in male blood and liver

Figure 5. qRT-PCR data depicting Gnas expression in livers from female (a) and male (b) mice. Gnas expression was determined relative to the geometric mean of the Ct values for housekeeping genes Actb, Hmbs, and Psmc4. Since some of the male and female mice were littermates, data were analysed using linear mixed effects regression with litter-specific random effects to account for within-litter correlation

Tissue and sex specificity of DMCs and DMRs

An important objective of the TaRGET II Consortium, and a need for the epigenetic research community at large, is to identify signatures of environmental exposures in surrogate tissues that mirror those in target tissues. To this end, we first compared the chromosomal locations for DMCs and DMRs between liver and blood to determine whether there was direct overlap in differentially methylated sites between the two tissues. In females, 4 DMCs and 3 DMRs overlapped between blood and liver; however, the Pb-induced change in methylation at each of the overlapping sites was in the opposite direction for each tissue (Table 6). In males, 4 DMCs and 5 DMRs each overlapped between blood and liver (Table 7). In contrast with females, several sites in males exhibited changes in methylation that were in the same direction in both blood and liver (Table 7). However, the overlaps of DMCs/DMRs between tissues were not statistically significant, as assessed by hypergeometric test (Tables 6–7). Next, we evaluated the overlap in genes with at least one DMC between blood and liver in males and females. We identified 188 DMC-associated genes in females and 325 genes in males that overlapped between blood and liver (Figure 6(a-b) and Supplemental Tables S19-S20). In both sexes, this overlap was statistically significant (p < 2.7x10⁻[36] and p < 1.1x10⁻[37] for females and males, respectively, hypergeometric test). Gene ontology analysis of the overlapping genes in females revealed enrichment for cell junction, membrane, and synapse pathways. In males, we observed enrichment for cell adhesion, protein binding, synapse, and neuronal process pathways (Figure 6(a-b)).

Table 6. Overlap between differentially methylated cytosines/regions (DMCs/DMRs) in female blood and liver

DMCs
Gene	Chromosome	Coordinate	Change in blood	FDR blood	Change in liver	FDR liver
Ccin	4	43,985,009	−11.08	0.00935	28.08	0.011
Lsm14b	2	180,023,832	−19.35	0.00035	16.92	0.017
Cntnap2	6	46,349,288	−17.51	0.00091	39.58	0.024
Fam122c	X	53,314,658	−18.76	0.00012	18.05	0.006
DMRs
Gene	Chromosome	Coordinate	Change in blood	FDR blood	Change in liver	FDR Liver
Glo1	17	30,597,451	11.86	0.00092	−13.38	0.045
Lsm14b*	2	180,023,801	−19.35	0.00026	16.92	0.012
Cntnap2*	6	46,349,251	−17.51	0.00070	39.58	0.018

Methylation changes are expressed as lead vs. control.

*DMRs comprise 50 bp genomic intervals, and some DMRs contain a single CpG. At these regions the change in methylation is therefore the same for both the DMC and DMR. Separate statistical tests were conducted to identify DMCs and DMRs.

Table 7. Overlap between differentially methylated cytosines/regions (DMCs/DMRs) in male blood and liver

DMCs
Gene	Start	End	Change in blood	FDR blood	Change in liver	FDR liver
Grifin	140,567,250	140,567,250	−12.83	0.025	−25.08	0.03932
Elmsan1	84,151,896	84,151,896	15.08	0.0005	11.74	0.00020
Amn	111,270,178	111,270,178	19.13	0.033	11.97	0.00060
Tbc1d30	121,294,168	121,294,168	22.54	0.047	27.76	0.01649
DMRs
Gene	Start	End	Change in blood	FDR blood	Change in liver	FDR liver
Grifin*	140,567,201	140,567,250	−12.83	0.020	−25.08	0.026
Mob2	142,016,401	142,016,450	−23.64	0.038	21.0	0.021
Dnm3	162,453,101	162,453,150	12.63	0.001	−23.95	0.049
Plekhg3	76,535,801	76,535,850	17.04	0.017	17.14	0.047
N/A	114,931,401	114,931,450	21.20	0.022	13.43	0.039

Methylation changes are expressed as lead vs. control.

DMRs comprise 50 bp genomic intervals, and some DMRs contain a single CpG. At these regions the change in methylation is therefore the same for both the DMC and DMR. Separate statistical tests were conducted to identify DMCs and DMRs.

Figure 6. Analysis of overlapping DMC-associated genes in females (a) and males (b). Panels C-D denote genes in blood (c) and liver (d) that are in common across sexes. Panel E depicts the number of genes in common across sex and tissue. Numbers in circles indicate the total number of DMC-associated genes in each tissue. All overlaps were statistically significant, as indicated by hypergeometric test (panel A: p < 2.7x10⁻⁴⁶, panel B: p < 1.1x10⁻⁷³, panel C: p < 9.459x10⁻⁵⁹, panel D: p < 2.656x10⁻⁹⁰, panel E: p < 3.227x10⁻³⁰). In the tables, # Genes indicates the number of genes in the specified GO pathway that appear in the experimental dataset. Pathway analysis was conducted using DAVID, version 6.8 (Laboratory of Human Retrovirology and Immunoinformatics (LHRI))

We next compared the sex-specific overlap between genes in both tissues. In blood and liver, respectively, we observed 186 and 354 genes in common (Figure 6(c-d) and Supplemental Tables S21-S22). Although the intersecting genes were a minority of the total number of DMC-associated genes, the overlaps were statistically significant (p < 9.46x10⁻⁵⁹ and 2.66 × 10⁻⁹⁰ for blood and liver, respectively, hypergeometric test). We then investigated the intersection of DMC-associated genes across tissues and sexes, and found a small number of genes (35 genes) were common across tissues and sexes (Figure 6(e) and Supplemental Table S23). This was significantly greater than what would be expected by chance p < 3.23x10⁻³⁰, hypergeometric test). These included several genes associated with obesity, metabolic syndrome and hepatocellular carcinoma, including Prdm16 [38–40], Fto [41,42], and Wwox [43].

Pb-induced changes in gene expression

In order to determine whether altered DNA methylation was accompanied by changes in expression of Fto, Wwox and Prdm16 in both tissues, we interrogated publicly available RNA-seq read count data from the same animals that we generated as part of the TaRGET II consortium (see Methods). Notably, Fto expression was significantly increased in female liver, but not blood at 5 months of age (Figure 7(a), top). Males showed no change in Fto expression in either tissue at this time point (Figure 7(b), top). In contrast, females showed no change in Wwox expression at 5 months of age, but expression in males was increased with Pb exposure in both tissues in adulthood (Figure 7(c-d), top). Cessation of Pb exposure occurred at postnatal day (PND) 21, and we wondered whether changes in gene expression were present at this early time point. Interestingly, Fto expression was significantly increased in male, but not female liver at PND21 (Figure 7(a-b), bottom). Wwox expression was not significantly changed at PND21 in either sex or tissue (Figure 7(c-d), bottom).

Figure 7. Effects of Pb exposure on gene expression in mouse blood and liver at two time points. Normalized RNA-seq read count data for females (a,c) and males (b,d) for Fto and Wwox. Top graphs in each panel depict gene expression in a cohort of animals sacrificed at 5 months of age (the same animals used for ERRBS), and bottom graphs show gene expression at postnatal day (PND) 21, when animals were weaned and Pb exposure stopped. A cohort of animals were sacrificed for gene expression analysis at this time point. Statistical significance was determined using linear mixed effects regression

Differential methylation at sexually dimorphic genes

Given the established sex-biased gene expression the liver [44,45], we next investigated whether genes exhibiting sex-specific methylation with Pb exposure were among those that show basal sex differences in expression in the adult mouse liver. To this end, we compared our findings with those of a previously published study that examined sex differences in mouse hepatic gene expression at multiple time points [44]. Of the 1422 genes containing DMCs with Pb exposure in the female liver, 46 were significantly differentially expressed between males and females at 8 weeks of age (Figure 8(a) and Supplemental Table S24). Likewise, 78 of the 1999 genes showing differential methylation with Pb exposure in males were also differentially expressed in the male vs. female adult mouse liver (Figure 8(b) and Supplementary Table S25). These overlaps were statistically significant (p = 0.001 and 0.02 for females and males, respectively). A subset of these genes was also differentially methylated in the blood, including 2 in common across tissue and sex (Figure 8(a-b)). Together, these data suggest that there exist sex-specific changes in DNA methylation with perinatal Pb exposure that are consistent across blood and liver in adulthood, with significant gene overlap between tissues and sexes. Furthermore, several Pb-induced changes in methylation in both males and females occur at sexually dimorphic genes.

Figure 8. Overlap between genes with differential methylation (DM genes Pb) and genes that are differentially expressed in male vs. female liver at 8 weeks of age (DE genes 8 weeks). Separate analyses were conducted in females (a) and males (b). Genes listed below each diagram are those that also exhibit differential methylation in the blood of either females (A) or males (B). Bold genes are differentially methylated with Pb exposure in both tissues and sexes, and are differentially expressed between sexes in 8 week mouse liver. P values (hypergeometric test) represent the significance of overlap between gene lists

Discussion

In this work, we have investigated the effects of perinatal Pb exposure on DNA methylation in the liver and blood of adult mice. The findings from this study are novel for several reasons. First, although increasing evidence links Pb exposure to adverse metabolic outcomes [46], the mechanistic basis for these associations is unclear. Indeed, little is known about the effects of Pb exposure on the liver. To our knowledge, this is the first report demonstrating that Pb exposure during pregnancy and lactation leads to sex-specific alterations in hepatic DNA methylation in adulthood. Second, we discovered that perinatal Pb exposure leads to changes in methylation at several imprinted loci that are also reprogrammed by perinatal exposure to another common toxicant, BPA [4]. Finally, we mapped DNA methylation genome-wide in paired blood and liver samples from the same mice after perinatal Pb exposure. Few differentially methylated sites directly overlapped between blood and liver; however, we discovered a small subset of DMC-associated genes with significant overlap between tissues in both males and females, including several that play critical roles in obesity and metabolic syndrome. This finding has important implications for the design and interpretation of human population-based epigenomics studies.

Tissue- and sex-specific changes in DNA methylation with Pb exposure

Although Pb exposure ceased at 3 weeks of age, we discovered thousands of tissue- and sex-specific DMCs/DMRs in both male and female mice at 5 months of age. The sex-specific effects of Pb on DNA methylation have been reported by our lab and others [12,13,47]. Indeed, a subset of differentially methylated regions were found in published sex-specific mouse liver enhancers [30]. As sex-specific chromatin structure drives sex differential gene expression, it is plausible that Pb exposure may reprogram the sexually dimorphic patterns of gene expression in the liver. Consistent with this, we observed altered DNA methylation at a subset of genes that are differentially expressed between male and female mice at 8 weeks of age [44]. Likewise, when comparing males and females, pathway analyses revealed enrichment of distinct biological pathways. These data suggest that sex will be a critical factor in the design, analysis, and interpretation of human- and animal-based environmental epigenetics studies.

Although we investigated the effects of perinatal Pb exposure on DNA methylation in liver and blood, pathway analysis of DMRs in promoter regions revealed enrichment for pathways involved in neurological development and function in male and female liver. Moreover, among overlapping genes between blood and liver, gene pathways associated with synapse function were significantly enriched in both males and females. The nervous system is a well-known target of Pb, and although we did not directly measure nervous system tissues in this study, it is nonetheless striking that neurological pathways were enriched in genes annotated to DMCs/DMRs in both liver and blood. However, it is also important to note that the liver receives significant innervation from the autonomic nervous system, which functions in the sensing of glucose and lipid concentrations and maintenance of metabolic homoeostasis[48]. Intriguingly, several nervous system-related genes identified in the pathway analysis are also expressed in the liver in normal or disease states. For example, Slc18a3, Rab3b, and Nos1, are expressed in normal liver, although their hepatic functions are not well-characterized [36,49,50]. Likewise, Slc6a8, Tenm4, Celsr1, and Dvl1 have all been reported to be mutated, hyper-methylated or over-expressed in hepatocellular carcinoma or colon cancer liver metastases [51–54]. Analysis of other tissues, including brain, from Pb-exposed mice is currently underway, and will provide important insight into whether there is a common subset of tissue-independent, environmentally labile genes or CpG sites.

Imprinted genes

Genomic imprinting, which results in mono-allelic, parent-of-origin expression of a subset of genes, is mediated in large part by DNA methylation, and establishment of 5mC at imprinted loci occurs during early embryonic development [55,56]. Importantly, we and others have reported that environmental exposures during this critical developmental window lead to altered 5mC and 5hmC at imprinted loci [4,11,34,57,58] and alterations in imprinted gene expression [4,58]. Imprinted genes play a critical role in growth and metabolism, and abnormal DNA methylation and expression of imprinted genes are associated with several diseases, including growth, metabolic and developmental disorders, and cancer [59]. Notably, we discovered that Pb exposure during gestation and lactation leads to changes in DNA methylation at 44 imprinted genes, and this enrichment is statistically significant. The Gnas locus encodes for three distinct protein products, and disruptions in the regulation of this locus have been implicated in a number of disorders, including early-onset obesity [60]. We recently reported that developmental BPA exposure leads to changes in DNA methylation and hydroxymethylation at the Gnas locus, concomitant with increased Gnas expression [4]. In this study, hypermethylation at the Gnas locus in both males and females was accompanied by reduced Gnas transcript expression. Thus, Gnas may be an important target of multiple environmental exposures, in both males and females. When comparing male liver and blood, we observed DMCs at several imprinted genes. For example, Snrpn exhibited changes in methylation in the same direction in both blood and liver. Snrpn encodes for an RNA binding protein that is critical for neurological function, and disruption of Snrpn expression is implicated in Prader-Willi syndrome, a neuroendocrine disorder associated with intellectual disability and obesity [38]. Work from us and others suggests that environmental exposures are associated with changes in methylation of the Snrpn locus [4,58,61,62]. Given that the Srnpn locus is hypermethylated in both blood and liver with gestational Pb exposure, it is plausible that it may represent an important biomarker of early life Pb exposure. We recently reported that perinatal BPA exposure leads to reprogramming of DNA methylation and hydroxymethylation at 12 imprinted genes [4]. Interestingly, 9 of the imprinted loci altered by Pb exposure were among those reprogrammed by BPA exposure. These data suggest that a core set of imprinted genes may be vulnerable to epigenetic reprogramming by multiple environmental exposures. Recent work suggests that the late stages of oocyte development and early embryonic development are particularly vulnerable to altered genomic imprinting by environmental exposures [58]. As both Pb and BPA exposure have been implicated in adverse metabolic health outcomes [46,57,63], it is plausible that reprogramming of genomic imprinting may represent a common mechanism linking both toxicants to these outcomes. The mechanistic consequences of environment-mediated dysregulation of imprinted genes, and the implications this may have for human health, should be investigated in future studies.

Potential mechanisms for altered DNA methylation

The mechanism(s) by which developmental Pb exposure alter DNA methylation are currently unknown. DNA methyltransferases (DNMTs) methylate the 5-position of cytosine bases in DNA using S-adenosylmethionine (SAM), which is generated via one-carbon metabolism [64]. Pb may disrupt DNA methylation via alteration in DNMT expression/activity or through perturbation of SAM levels. Consistent with the first hypothesis, Pb exposure in vitro and in vivo leads to inhibition and altered expression of DNMTs [65]. In support of the second possibility, administration of exogenous SAM mitigates the deleterious effects of Pb exposure, suggesting that Pb may deplete endogenous levels of SAM [66,67]. DNA methylation can be erased through an active process that includes 5hmC as an intermediate, and altered 5hmC has also been reported with Pb exposure [68]. 5hmC is generated by the activity of TET methylcytosine dioxygenases (TETs), using the cofactor alpha-ketoglutarate (α-KG), which is generated by the TCA cycle [69]. Interestingly, Pb exposure has been shown to reduce activity of TCA cycle enzymes, including isocitrate dehydrogenase, which catalyzes the formation of α-KG [70,71]. In spite of the substantial differences between Pb and BPA, similar mechanisms have been reported for BPA. Indeed, the effects of developmental BPA exposure on DNA methylation can be countered by methyl donor supplementation [72]. Likewise, BPA has been reported to alter expression of DNMT1 in other tissues [37,73]. These mechanistic similarities, combined with the similar developmental window of exposure, may explain why BPA and Pb exhibit effects on DNA methylation at some common gene loci.

It is additionally notable that Pb-induced changes in DNA methylation are sex-specific. Although DNA methylation undergoes reprogramming during fertilization and primordial germ cell development, the kinetics of, and mechanisms underlying, this reprogramming differ between male and female genomes [2,74]. In the preimplantation embryo, demethyation of the paternal genome occurs via an active, TET-mediated process [74]. In contrast, demethylation of the maternal genome occurs through a slower, passive loss of DNA methylation over the course of cell divisions [74]. In the primordial germ cells, similar sex differences occur during the restoration of sex-specific DNA methylation. In sperm, this process is completed prior to birth; however, in oocytes, sex-specific remethylation continues until puberty [2,74]. Thus, sex differences in Pb-mediated alterations in DNA methylation may be due to these mechanistic and temporal differences. The precise mechanisms underlying the observed Pb-induced changes in DNA methylation are currently under investigation.

Utility of blood as a surrogate tissue for Pb effects on liver

When comparing the chromosomal locations of DMCs/DMRs between blood and liver, only a nonsignificant number of sites directly overlapped between the two tissues. However, we observed a significant number of DMC-associated genes in common between blood and liver, and 35 genes showed differential methylation across tissues and sexes. These overlaps were statistically significant. Importantly, altered methylation at Fto and Wwox was accompanied by changes in gene expression which were time, tissue and sex-dependent, highlighting the complex nature of sex- and tissue-specific gene regulation. Fto expression was elevated in the female liver at 5 months, and in the male liver at PND21. Moreover, Wwox expression was elevated in blood and liver of males at 5 months of age, suggesting that it may serve as a biomarker for Pb exposure. SNPs in the FTO gene are associated with increased energy intake and obesity in human and rodent studies [42]. WWOX is a tumour suppressor that is lost in hepatocellular carcinoma and is also de-regulated in metabolic syndrome [43,75]. Whether altered regulation of these genes plays a mechanistic role in Pb-induced adverse metabolic outcomes requires further investigation. The genes identified herein, particularly those that are commonly differentially methylated and expressed in blood and liver, may represent candidates for analysis in human environmental epigenetics studies. Future studies are needed to validate the plausibility of these loci as potential biomarkers of Pb exposure, as well as to investigate the functional and health-related consequences of altered DNA methylation at these sites.

Limitations

This study has several important limitations. First, we utilized ERRBS, a method that, while providing resolution at the level of individual CpG sites, does not cover the entire genome. The use of whole-genome bisulphite sequencing techniques may yield additional insight into the effects of Pb exposure on the liver, and identify additional Pb-induced signatures that overlap between blood and liver. Second, ERRBS does not differentiate between 5mC and 5hmC. Thus, the differentially methylated loci identified in this study may be a combination of both 5mC and 5hmC. Third, changes of interest in DNA methylation may only occur in specific subpopulations of cells [76]. Thus, our approach using unfractionated tissues may have missed other potential markers of interest for Pb exposure. Additionally, blood cells undergo rapid turnover, as well as exposure-, age-, and disease-dependent changes in cellular composition [77–79]. Future studies using single-cell transcriptomics and epigenomics approaches will help to address issues of tissue-specific cellular heterogeneity and composition [76,80]. It is also important to note that DAVID gene ontology analysis does not consider the differing number of CpG sites per gene, and therefore may exhibit bias. Although we found that altered methylation was accompanied by changes in expression of a few genes, we observed few significant changes in gene expression with RNA-seq (reference https://targetepigenomics.org/). Thus, the functional significance of altered methylation is unclear based on this study. However, recent work suggests that developmental exposures may induce epigenetic changes that are silent until a further physiologic stimulus or environmental insult is encountered [81,82]. Finally, as we did not measure Pb levels in the blood and liver of offspring mice, we cannot exclude the possibility that some Pb remains in these tissues at 5 months of age. Thus, some changes in methylation and gene expression observed at this time point may be due to ongoing exposure to residual Pb and not to a programming event that occurred during early development. However, the majority of Pb is absorbed in bone, with a small fraction found in blood and soft tissues [83,84]. While Pb is stably stored in bone, children are less likely to retain Pb in bone compared to adults, and turnover is much more dynamic in soft tissues [83,84]. Ongoing release of Pb from bone is of concern; however, this occurs primarily during pregnancy, stress, and age-related bone loss [85,86]. Thus, given the early exposure window, the relatively young age of the adult mice and the dynamic turnover of Pb in soft tissues, we would expect little Pb to be present in blood and liver at 5 months of age. Nevertheless, whether Pb remains in the body after early developmental exposure or not, our data suggest that Pb exposure during early life causes epigenetic and gene expression changes that may have important implications for health across the life course.

Conclusions

In summary, we have investigated the effects of perinatal Pb exposure on DNA methylation in paired samples of mouse blood and liver. There are several key strengths of this study, including the side-by-side comparison of blood and liver to identify signatures of Pb exposure present in the blood that reflect changes in the liver. Future studies should investigate the reproducibility of these findings, the mechanistic and functional consequences of persistent, Pb-induced changes in hepatic methylation, and the implications these changes have for human health.

Author Contributions

All authors provided critical feedback and helped shape the research, analysis and manuscript. KN, CAR, TRJ, CL performed experiments. LKS, RC, CAR, ZT, SL, JMG, JAC, MAS, DCD analyzed data. LKS wrote manuscript in consultation with DCD, JAC, MAS, and JMG.

Disclosure statement

The authors declare no conflict of interest.

Data Availability

ERRBS sequencing data have been uploaded to GEO (accession number GSE150670). RNA-seq data are available from the TaRGET II Consortium: https://targetepigenomics.org/.

Supplementary material

Supplemental data for this article can be accessed here.

Additional information

Funding

This work was supported by the National Institute of Environmental Health Sciences [P01 ES022844/RD83543601]; National Institute of Environmental Health Sciences [ES026697]; National Institute of Environmental Health Sciences [ES028802]; National Institute of Environmental Health Sciences [T32 ES007062]; National Institute of Environmental Health Sciences [P30 ES017885].