Effects of cross-fostering and developmental exposure to mixtures of environmental contaminants on hepatic gene expression in prepubertal 21 days old and adult male Sprague-Dawley rats

ABSTRACT

The notion that adverse health effects produced by exposure to environmental contaminants (EC) may be modulated by the presence of non-chemical stressors is gaining attention. Previously, our lab demonstrated that cross-fostering (adoption of a litter at birth) acted as a non-chemical stressor that amplified the influence of developmental exposure to EC on the glucocorticoid stress-response in adult rats. Using liver from the same rats, the aim of the current study was to investigate whether cross-fostering might also modulate EC-induced alterations in hepatic gene expression profiles. During pregnancy and nursing, Sprague-Dawley dams were fed cookies laced with corn oil (control, C) or a chemical mixture (M) composed of polychlorinated biphenyls (PCB), organochlorine pesticides (OCP), and methylmercury (MeHg), at 1 mg/kg/day. This mixture simulated the contaminant profile reported in maternal human blood. At birth, some control and M treated litters were cross-fostered to form two additional groups with different biological/nursing mothers (CC and MM). The hepatic transcriptome was analyzed by DNA microarray in male offspring at postnatal days 21 and 78–86. Mixture exposure altered the expression of detoxification and energy metabolism genes in both age groups, but with different sets of genes affected at day 21 and 78–86. Cross-fostering modulated the effects of M on gene expression pattern (MM vs M), as well as expression of energy metabolism genes between control groups (CC vs C). In conclusion, while describing short and long-term effects of developmental exposure to EC on hepatic transcriptomes, these cross-fostering results further support the consideration of non-chemical stressors in EC risk assessments.

KEYWORDS:

• Rat
• liver
• developmental exposure
• organochlorines
• persistent organic pollutants
• mixtures
• cross-fostering
• gene expression microarray

Introduction

Humans are continuously exposed to mixtures of environmental contaminants (EC) through various sources including food, drinking water, and ambient air. In the United States and Canada pregnant women were found to contain detectable levels of EC including polychlorinated biphenyls (PCB), organochlorine pesticides (OCP), perfluorinated compounds, phenols, polybrominated diphenylethers (PBDE), phthalates, polycyclic aromatic hydrocarbons (PAH), and perchlorate (Canadian environmental contaminants biomonitoring programs 2018; Woodruff, Zota, and Schwartz 2011). Several epidemiological and clinical studies demonstrated a link between EC exposure and increased risks of health impairments and diseases including neurological [Parkinson’s (Chhillar et al. 2013; Ismail et al. 2017; Weisskopf et al. 2010) and Alzheimer’s (Richardson et al. 2014)], reproductive (Han et al. 2016; Karwacka et al. 2017; Lenters et al. 2015; Veras et al. 2010), developmental (Li et al. 2018; Torres-Sánchez et al. 2013; Williams, Watson, and Desesso 2012), respiratory (0Gascon et al. 2014), metabolic diseases (Vafeiadi et al. 2016), as well as cancer (Bencko et al. 2008; VoPham et al. 2015). Recent epidemiological investigations reported positive correlations between exposure to persistent organic pollutants (PCB, OCP), metabolic diseases, and dyslipidemia in populations of Northern Canada (Inuits, First Nations from Northern Manitoba and Ontario) (Marushka et al. 2018; Singh and Chan 2017, 2018). Although single contaminants or simple mixtures (two to three chemicals) are often studied (Cory-Slechta 2005; Jung et al. 2013), environmental exposures more typically involve low levels of complex mixtures of chemicals.

Cross-fostering (adoption of a litter at birth) is an animal study procedure used in neuroscience (Fish et al. 2004) and toxicology (Crofton et al. 2000) to isolate the consequences of in utero events, such as chemical exposure, from those occurring after birth. Previously, testing a mixture reconstituted based upon reported chemicals detected in maternal blood of Northern Canadians (Butler Walker et al. 2003; Muckle et al. 2001), Desaulniers, Xiao, and Cummings-Lorbetskie (2013) demonstrated that cross-fostering enhanced the adverse effects of developmental exposure to the mixture of EC on the glucocorticoid stress-response (GSR) in F1 (first generation) adult rats. In mixture-exposed groups, corticosterone concentrations remained elevated 30 min after stress created by tail bleed in the cross-fostered group, but not in those that were not cross-fostered. The selected doses were (1) low to intermediate relative to those of previous studies (Bowers et al. 2004; Chu et al. 2008); (2) did not markedly affect dam body weights or reproductive success; and (3) induced no overt signs of toxicity in the offspring. In addition, chemical residue levels in F1 adult males were comparable to those detected in humans (Desaulniers et al. 2017; Desaulniers, Xiao, and Cummings-Lorbetskie 2013). Using liver samples from the same rats, the objective of the current follow-up study was to determine the effect of cross-fostering in combination with developmental exposure to the EC mixture on hepatic gene expression profile at postnatal day 21 (PND21) and at PND78-86 in F1 males. The detection of a cross-fostering confounding effect on hepatic transcriptome might support our previous GSR observations (Desaulniers, Xiao, and Cummings-Lorbetskie 2013), and provide information for risk assessors to interpret other cross-fostering studies. Notably, it might improve evidence that non-chemical stressors (cross-fostering as a model) might modulate chemical-mediated effects.

Materials and methods

Animal exposure, cross-fostering procedures, and sample collection

Sprague-Dawley rats were acquired from Charles River Laboratories (St. Constant, Quebec, Canada). Treatment of rats and their housing conditions was approved by the Animal Care Committee of Health Canada and conformed to the Guidelines of the Canadian Council on Animal Care. Details on all animal manipulations, housing, acclimatization, breeding, assignment of females to treatment groups, mixture composition, dosing procedures with cookies, cross-fostering at birth, and culling procedures, were previously reported (Desaulniers, Xiao, and Cummings-Lorbetskie 2013).

The current manuscript focused on effects of cross-fostering on hepatic transcriptomes analyzed from rats of undisturbed (C, M) and cross-fostered litters (CC, MM), treated or not with the mixture (M), and reports data from F1 male pups at postnatal day (PND) 21, and in older siblings during the interval of PND78-86. The last dose was administered to the dams at PND20, and euthanasia at PND21 prior to weaning ensured maximal exposure without the potential confounding effect of weaning stress. Euthanasia of F1 adult males was performed during the interval of PND78-86 to ensure that the glucocorticoid stress-response data (Desaulniers, Xiao, and Cummings-Lorbetskie 2013) could be compared to the existing literature from young adult males (McCormick et al. 2000; Weaver et al. 2004).

Table 1 describes the treatment groups, while the supplementary Table S1 presents the detailed chemical content of the mixture M, which was prepared based upon the maternal environmental contaminant blood profile measured in Canadian Inuits (Butler Walker et al. 2003; Muckle et al. 2001). The dose of 1 mg/kg/day was selected for the cross-fostering experiment based upon a previous study indicating that 5 mg/kg/day was toxic but not 0.5 mg/kg/day (Chu et al. 2008). A discussion of the dose relative to the lowest observable effect levels and margin of exposure (MOE), and of the achieved residue levels, is provided in Supplementary Table S6. To test the effects of cross-fostering, some litters from C and M treatment groups were transferred to different mothers at birth to create the following in utero/postnatal exposure combinations, CC, and MM. The first letter represents in utero exposure from the biological mother, whereas the second letter reflects exposure during nursing by the adopting mother.

Table 1. Description of the group names, treatments, and cross-fostered groups (adapted from Desaulniers et al. 2013).

Groups undisturbed with no cross-fostering (dams treated from GD0-PND20)			Numbersof litters
1) Control (C)	Corn oil		13
2) M	Complex mixture at 1.0 mg/kg/day		10
Groups with cross-fostering^4:
	Biological mother In utero exposure GD0-PND0	Fostering mother Exposure during nursing PND0-PND21
3) CC	Corn oil	Corn oil	13
4) MM	M	M	9

Briefly, the solution of chemical mixture dissolved in corn oil was pipetted at a volume of 1 µl/g body weight and placed on cookies (Teddy Graham, Nabisco Ltd., Toronto, ON, Canada; approximately 2 g) that absorbed the entire solution. This procedure was developed as a stress-free dosing method for neurodevelopmental studies (Bowers et al. 2004). The volumes were adjusted daily for the dam’s body weight, and rats were dosed one cookie per day starting on gestation day 0 (identified by the presence of a vaginal plug or sperm in vaginal smear) until postnatal day (PND) 20. Dams were observed to ensure that the cookies were consumed entirely, and wastage was not an issue as mothers would eagerly eat the entire cookie immediately when given as treats apart from their normal diet. Livers were harvested and wrapped in aluminum foil, immediately flash frozen in liquid nitrogen, and stored at −80°C until RNA isolation.

RNA isolation and purification

RNA was purified from liver samples using the QIAGEN RNeasy (Qiagen, Mississauga, ON, Canada) extraction kit according to the manufacturer’s protocol. The RNase-free DNase (QIAGEN) kit was employed for on-column digestion of DNA during the RNA purification process. Purified RNA was eluted in 50 µl RNase-free water and stored at −80°C until use. The concentration of RNA was measured utilizing a NanoDrop ND-1000 spectrophotometer (NanoDropTechnologies, Wilmington, DE). The quality and integrity of RNA was assessed using an Agilent 2100 Bioanalyzer (Agilent Technologies Inc., Mississauga, ON, Canada). RNA samples with RNA integrity number (RIN) above 8.5 were used for subsequent microarray and RT-qPCR analyses.

Microarray analyses

Microarray experimental design

Samples from 8 male pups from 8 different litters per treatment group were randomly selected at PND21 and PND78-86. RNA from pairs of rats in the same treatment group were randomly pooled together and thus each treatment group was represented by 4 microarrays. The analysis was performed using Agilent 8X60K high-density slides (SurePrint G3 Rat GE 8x60k Kit # G4853A), which contained 8 microarrays per slide.

Cross-fostering effects are reported here based upon comparison C vs CC and M vs MM. However, the original study tested multiple mixtures with 9 treatment groups, among which 5 were undisturbed and tested dose-response effects and increasing abundance of aryl hydrocarbon-receptor (AhR)-agonists (C, 0.5M, M, AhR, 0.5MAhR), whereas 4 groups tested effects of cross-fostering (CC, MC, CM, MM) (Desaulniers, Xiao, and Cummings-Lorbetskie 2013). The groups MC, CM, at PND21 and at PND78-86, as well as the groups 0.5M and 0.5MAhR at PND78-86 were not examined using microarrays. Some microarray data from the groups 0.5M, AhR, and 0.5MAhR are presented here in the main document only to support the performance of the microarray assays (i.e. sensitivity to detect effects of AhR-agonists and dose-response effects across the mixtures AhR, 0.5M, 0.5MAhR, and M). The chemical content and discussion of the effects of these additional mixtures are described in the supplementary document. The microarray designs were as follows.

At PND21, 7 treatment groups were investigated (C, AhR, 0.5M, M, 0.5MAhR, CC, and MM). Samples (4 samples/group) from groups C, AhR, 0.5M, and M were hybridized using two slides with two microarrays from each treatment group on each slide. For the cross-fostering comparison, the 4 microarrays for group CC were compared to those of group C utilizing one slide, and similarly for the group MM and M comparison (4 microarrays per treatment group).

At PND78-86, effects of treatment on gene expression were expected to be subtle given that the rats were exposed only developmentally, therefore comparisons were established within slides. Three aliquots were prepared from group C, and treatment groups compared as follows: group M and C were hybridized to one slide; group AhR and C on a second slide; group CC and C on a third slide; and finally, group MM and M on a last slide. Agilent Feature Extraction Software version 11.0.1.1 (Agilent Technologies) was employed for data extraction and complete dataset may be found in the Gene Expression Omnibus (http://www.ncbi.nlm.nih.gov/geo/, accession number GSE109019).

cDNA synthesis and array hybridization

Total RNA (200 ng) from exposed and control groups, as well as universal reference total RNA (Agilent Technologies) were employed to synthesize complementary DNA (cDNA) and cyanine labeled complementary RNA (cRNA), using Low Input Quick Amp Labeling Kit (Agilent Technologies) according to manufacturer’s instruction. Cyanine-labeled cRNAs were synthesized using T7 RNA polymerase and in vitro transcription kits (Agilent Technologies). cRNA was then purified using RNeasy Mini Kits (Qiagen). Cyanine 5-CTP (Cy-5) was used to label the sample, and universal reference RNAs were labeled with cyanine 3-CTP (Cy-3).

Cy-3 labeled reference cRNA (300 ng) was then mixed with equal amounts of Cy-5 labeled cRNA samples and hybridized to microarrays (Agilent Technologies) for 17 hr in a hybridization chamber at 65°C with a rotation speed of 10 rpm. Microarray slides were washed and scanned on an Agilent G2505C scanner at 3 μm resolution, according to manufacturer’s instructions.

Data analysis

A reference design (Kerr 2003; Kerr and Churchill 2001) with the sample labeled in Cy5 and universal reference RNA labeled in Cy3 was used to analyze data. All data analyses were conducted using the R (R Core Team 2015) statistical environment. The median signal intensities were normalized using the global LOcally WEighted Scatterplot Smoothing (LOWESS) method (Yang et al. 2002) using the MAANOVA library (Wu et al. 2002). Differentially expressed genes (DEG) were identified using the MAANOVA library. The ANOVA model included the main effect of treatment and array as a block effect. The Fs statistic (Cui et al. 2005), a shrinkage estimator, was employed for gene-specific variance components and associated p-values for all the statistical tests were estimated using the permutation method with residual shuffling (30,000 permutations). The p-values were then adjusted (FDR p-values) for multiple comparisons using the false discovery rate approach (Benjamini and Hochberg 1995) to reduce false positives. The least-squares means (Goodnight and Harvey 1978; Searle, Speed, and Milliken 1980) were utilized to estimate the fold change (FC) for each pairwise comparison. Pairwise comparisons for individual genes were considered if all samples (n = 4) for one of the groups had signal intensities above the background distribution estimated by the trimmed mean plus 3 standard deviations of the negative control probes (3xSLv1) (Dong et al. 2007). To elucidate the influence of chemical mixtures the following comparisons were conducted: M vs C, 0.5M vs C, AhR vs C, 0.5MAhR vs C at PND21, and M vs C, AhR vs C at PND78-86. To determine the effect of cross-fostering, CC vs C and MM vs M for both age groups were tested.

To test the consistency of microarray results, hierarchical cluster analysis (HCA) and principal component analysis (PCA) were conducted using the DEG (FDR p-value <0.05; FC > |1.5|). A scatterplot of the first two principle components was generated. Each point on the scatterplot represents one sample (i.e. one array). This enabled visualization of the consistency between replicates and differences between treatment groups. With the list of DEG, an HCA was conducted using average linkage and the one minus the Spearman correlation dissimilarity metric. To analyze the effects of cross-fostering by PCA, data for the control (C) and the cross-fostered control samples (CC) were centered using the median of the control samples (C). Similarly, the mixture (M) and cross-fostered mixture samples (MM) were centered utilizing the median of mixture samples (M). This approach, as well as having the C and CC microarrays on one slide and the M and MM groups on a different slide, was adopted to optimize sensitivity given that only small effects of cross-fostering were expected. All analyses were visualized employing heat maps and dendrograms.

Functional analysis

In order to interpret the results in the context of biological processes and identify regulatory networks potentially affected in the liver of chemical-treated and cross-fostered rats, bioinformatics analysis was performed using Qiagen’s Ingenuity® Pathway Analysis (IPA®, Qiagen Redwood City, www.qiagen.com/ingenuity). The DEG with their log₂ fold changes and FDR p-values (Table S2) for each comparison were analyzed in IPA using upstream regulator analysis, canonical pathway, and network analysis. In these analyses, the SurePrint G3 Rat GE 8x60k microarray (V2) was employed as a reference set for transcriptomics data and included all data sources from Ingenuity Knowledge Base (IKB), while restricted to specific species (Rat and Mouse) and tissues (liver and cell lines).

Validation of microarray data by quantitative real-time PCR (RT-qPCR)

DEG were selected for RT-qPCR validation (Suppl. Table S3) using custom-made PCR arrays containing proprietary primers (RT² Profiler PCR Array: CAPR10984, Qiagen and PrimePCR SYBR green assays, Bio-Rad). The RNA samples used in the microarray experiment were also employed for RT-qPCR. Total RNA was reverse transcribed utilizing the RT² First Strand cDNA kit according to the company’s protocol. RT-qPCR was performed on a C-1000 thermal cycler (Bio-Rad) according to the protocol provided by the RT² SYBR green kit (Qiagen/Bio-Rad). Expression values were normalized to the mean of the three endogenous controls: β-actin (Actb), Ribosomal protein large P1 (Rplp1), and Ribosomal protein L13A (Rpl13a). Based upon microarray and PCR data, treatments did not markedly affect these endogenous control genes. The quality of the Ct values was assessed using online software provided by SABiosciences at http://pcrdataanalysis.sabiosciences.com/pcr/arrayanalysis.php. Fold changes were calculated according to the 2^−ΔΔCt method (Livak and Schmittgen 2001).

Results

This manuscript focused on the effects of cross-fostering based upon comparisons M vs C, CC vs C, and MM vs M in the PND21 and PND78-86 F1 males. Partial results related to the effects of the other mixtures (AhR, 0.5M, 0.5MAhR) are presented here only to support the assay performance; however, the complete dataset and discussion on the effects of these mixtures are presented in the supplementary document.

Effects of mixtures as indication of assay performance

The reliability of the microarray data is supported by the dose-response effects observed across treatment groups based upon (1) magnitude in gene expression, (2) the numbers of DEG, and (3) PCA and HCA. The Cyp gene family, particularly Cyp1a1, is inducible by many of the organochlorines present in the mixture, and provide an example of a dose-response pattern detected by the microarray analyses at PND21 (Figure 1). The fold changes relative to control for Cyp1a1 were 5.5, 13.9, 20.9, and 25.7 in the AhR, 0.5M, 0.5MAhR, and M groups, reflecting the 2,3,7,8-tetrachlorodibenzo-p-dioxin-toxic equivalent (TCDD-TEQ) doses of 1.8, 4.1, 6.3, and 8.2 ng/kg/day in the AhR, 0.5M, 0.5MAhR, and M group, respectively (Table S1). The fold changes for Cyp2b2 were 3.7, 4.2, and 8.7 in the 0.5M, 0.5MAhR, and M groups, respectively. Other Cyp genes (1a2, 2a3, 2c6, and 2c7) were expressed to a greater extent in the M than in 0.5M or 0.5MAhR groups, also reflecting a dose-response pattern. These fold changes and other alterations are provided in the supplementary Table S2, which also indicates accession numbers, gene symbols, names, and FDR p-values for the list of all DEG.

Figure 1. Fold changes in Cyp genes found to be differentially expressed (FDR p ≤ 0.05; FC ≥ |1.5|) by microarray analyses in AhR, 0.5M, 0.5MAhR, and M treated groups relative to the control group at PND21. Bars are missing in absence of a statistically significant effect (data in Table S2 Loc293989 = Cyp2c6).

Dose-response changes in the number of DEG (FDR p ≤ 0.05; FC ≥ |1.5|) were observed at PND21 (Figure 2), with one, 4, 37, and 56 DEG in the AhR, 0.5M, 0.5MAhR, and M groups, respectively, with the chemical dose increasing from 39.8 ng/kg/day for the AhR mixture, to 1 mg/kg per day for the mixture M (Table S1). Secondly, the AhR treatment was associated with a larger number of DEG at PND78-86 (n = 31) than at PND21 (n = 1), suggesting a delay in the induction of DEG.

Figure 2. Numbers of differentially expressed genes (DEG) (FDR p ≤ 0.05; FC ≥ |1.5|) between the control and the AhR, 0.5M, 0.5MAhR, M, and the CC groups, as well as between the M and the MM group, from the microarray analyses in the male rat livers at PND21 and PND78-86. The numbers of DEG showed a dose-response effect at PND21, reflecting the amount of chemicals and TCDD-TEQ present in the mixtures (1.8, 4.1, 6.3, and 8.2 ng/kg/day TCDD-TEQ for the AhR, 0.5M, 0.5MAhR, and M groups, respectively). The groups 0.5M and 0.5MAhR at PND78-86 were not analysed by microarray.

The PCA and HCA applied to the lists of DEG demonstrate the reproducibility of the gene expression profiles between rats within treatment groups, and ability of microarrays to detect treatment-related effects (Figure 3). The PCA of the AhR, 0.5M, and M vs C comparisons at PND21 (Figure 3(a)) demonstrated that 71.3% of the variability in the dataset was explained by the first two principle components (PC). The first PC accounted for 64.5% of the variability, mostly as a contrast between the M and C treatment groups, and the second PC 6.8% of the variability. The heatmap (Figure 3(b)) illustrates that samples within treatment groups cluster together indicating high correlations for those samples within the same treatment groups. Similar to PND21, the PCA at PND78-86 differentiated the treatment groups (Figure 3(c)) mostly due to a contrast between the C and the M group. Here, PC1 and PC2 accounted for 64.8% (PC1: 52.1% + PC2: 12.7%) of the variability in the dataset. The HCA shows that samples of the M treatment clustered generally together with the highest degree of similarity. Effects measured at PND78-86 in the AhR and M treatment groups are due to the remaining chemical residues at that time (Desaulniers et al. 2017), and to the chronic effects of developmental exposure that may generate data variability.

Figure 3. Principal Component Analysis (PCA) represented by two-dimensional scatter plot (Figure 3(a,c)), and hierarchical cluster analysis (HCA) (Figure 3(b,d) based on the list of DEG (FDR p ≤ 0.05; FC ≥ |1.5|) from the comparisons of the treatment groups at PND21 (Figure 3(a,b), and at PND78-86 (Figure 3(c,d). Treatment groups were manually circled to facilitate identification and clustering in the PCA plots. In the hierarchical cluster analyses the rows represent the probes for genes with increased (red) and decreased (green) expressions (black unchanged relative to controls).

Overall, the dose-response effects of the mixtures on Cyp expressions, on (1) number of DEG, (2) PCA and HCA, as well as (3) list of DEG validated by RT-qPCR (Table S3), demonstrate the reliability of the methodologies.

Influence of cross-fostering by Venn diagram analyses

Figure 4(a) uses Venn diagram to illustrate that among the 45 genes affected by cross-fostering at PND21 (CC vs C), three genes (Cyp2c7, and Ubd in region B, and Sult1b1 in region D) were also differentially expressed in the M vs C comparison. However, the mechanisms affecting these three genes by treatment M or by cross-fostering may not be equivalent; while the fold change for Sult1b1 was in the same direction (up-regulated) for the CC and M groups, it was in the opposite direction for Cyp2c7 and Ubd (increased by M but reduced by cross-fostering in the CC group, Table S2). Figure 4(a) also shows that cross-fostering modulated expression of 16 genes in rats treated with M (evident in the MM vs M comparison). Cross-fostering either diminished the magnitude of the M-induced expression of genes, as observed with Cyp2c7 and Ubd (region B), or amplified the effect of M (in one direction or the other) generating a list of 14 DEG (region A) that were not significantly altered in the M vs C comparison (region E). The expression of the aryl hydrocarbon-receptor nuclear translocator-like gene (Arntl) was not markedly affected in the comparison M vs C at PND21 (FC = −1.06), but significantly lower (1.8-fold) in the MM group relative to M (Table S2).

Figure 4. Venn diagram illustrating DEG (FDR p ≤ 0.05; FC ≥ |1.5|) in the comparisons of group (a) CC vs C, M vs C, and MM vs M at PND21; and (b) CC vs C, M vs C, and MM vs M at PND78-86. In Fig. 4A, the 14 genes in region “a” were: Akr1c1, Akr1c12, Arntl, Cyp2d3, Cyp2j3, Dbp, Hsd17b6, Loc689316, Lcn2, Ppapdc1a, Retsat, Rup2, Tff3, Tpt1. The two genes in region “b” were Cyp2c7 and Ubd, while the 42 genes in region “c” were: Abcb11, Actn3, Ccl21, Cmya5, Cox8b, Cpa1, Cpa2, Cpb1, Cpt1b, Cryab, Ctrb1, Des, Eno3, Fabp3, Fdps, Fgb, Fhl1, Gck, Hamp, Hrc, Igtp, Jsrp1, Klkb1, Loc688173, Lmod2, Mgc108823, Mb, Myh8, Myl4, Mylpf, Pla2g1b, Ppp1r1a, Rnase1, Slc25a30, Slc35c2, Spink3, Tnnc1, Tnnc2, Tnni2, Tnnt1, Tpm1, Upk1b. Region “d” included Sult1b1 and region “e” included 53 genes: A2m, Alas1, Aldh1a1, Aldh1a7, Aldh1b1, Apoa2, Ces1c, Ces2a, Ces2i, Ces2j, Cesl1, Cyb5a, Cyp17a1, Cyp1a1, Cyp1a2, Cyp2a1, Cyp2a3, Cyp2b2, Cyp2c37, Cyp3a23/3a1, Defa10, Entpd5, Ephx1, Gsta3, Gsta4, Gsta5, Gstm2, Loc100359554, Loc293989, Loc310902, Loc494499, Loc680955, Maged2, Mmd2, Nat8, Nfil3, Nqo1, Nr0b2, Onecut1, Pbld, Pir, Prkcdbp, Rgd1564906, Ratnp-3b, Selenbp1, Slc16a11, Slc34a2, Tsku, Ugt1a6, Ugt2a3, Ugt2b1, Ugt2b37. In Figure 4(b), the 16 genes in region “a” were: Aldh1a7, Cyp2b2, Cyp2c12, Fgl1, Hbb-b1, Inhbc, Loc305806, Loc689316, Oas1a, Oas1k, Pcolce, Rt1-n1, Rab17, Srd5a1, Tfr2, Wfdc3. The ten genes in region “b” were: Acot5, Apoa5, B3galt1, Cdh17, Dbp, Rgd1560527, Rt1Aa, Slc13a5, Sult2al1, Ubd. The two genes in region “c” were Akr1c12 and Akr1c2 while the one in region “d” was Slc35c2. Finally, the 59 genes in region “e” were: Abcg5, Acsm3, Akr1c1, Bhlha15, C6, Casp12, Cav1, Cidea, Ctnnal1, Ctnnd1, Cyp4a2, Dcn, Dusp26, Efemp1, Efna1, Eif4ebp3, Esm1, F9, Fabp4, Fam111a, Fgf21, Frzb, Gfra1, Gfra2, Gpm6a, Gpsm1, Hist1h1d, Hspa1b, Hspb1, Igfbp6, Inmt, Isg20, Kcnk2, Loc287167, Loc310926, Loc501437, Loc678704, Loc690226, Loc690435, Lox, Lrp5, Ltc4s, Mlc1, Myo7b, Nfatc4, Nrarp, Parp14, Psd, Rgd1561849, Rgd1564865, Rt1-bb, Rgs4, Sdr16c6, Serpinb7, Tap2, Thumpd1, Tpd52l1, Usp9x, Vom2r34.

The list of DEG obtained by the CC vs C comparison at PND21 (Figure 4(a)) varies from that obtained at PND78-86 (Figure 4(b)), indicating that cross-fostering yielded different hepatic gene expression profiles in both age groups. This observation is analogous to the long-term effects of the mixture M, which induced different DEG at PND21 than at PND78-86. Figure 4(b) reveals that three DEG were detected in the CC vs C and M vs C comparisons with similar fold changes in magnitude and direction (Akr1c12 and Akr1c2 in region C with decreased expression, and Slc35c2 in region D with enhanced expression). The MM vs M comparison demonstrates that cross-fostering modulated the influence of M. The Akr1c12 and Akr1c2 (region C) expression levels were reduced by M (in M vs C comparison, 1.7- and 2.1-fold down-regulated, respectively), while expression of these genes rose in group MM relative to M (1.9- and 2.3-fold up-regulated, respectively) (Table S2), supporting that cross-fostering antagonized the effects of M for these two genes.

Influence of cross-fostering by principal component and hierarchical cluster analyses

The PCA at PND21 (Figure 5(a)) illustrated that there is less of a clear separation between treatment groups, although 62.2% of the variability in the datasets was attributed to the first two PC (PC1 48.0% + PC2 14.2%). Here, PC1 might account for the variability of the two CC arrays compared to the rest of the samples. PC2, on the other hand, appears to be a contrast between two MM samples and the rest of the data. The HCA (Figure 5(b)) shows two CC and two MM samples clustering early, but variation among samples is higher, and lower correlations were found within treatment groups. The C and M groups appear uncorrelated, but as indicated in the method, in order to compare the cross-fostering effect from CC to C with MM to M comparison, data were further normalized using their respective controls. In this analysis, the C and M conditions represent the baseline expression centered at zero. Thus, the distribution of the C and M samples in the cross-fostering analyses (as well as at PND78-86) represents uncorrelated random noise. Overall, at PND21 both PCA and HCA suggest that cross-fostering effect was subtle and perhaps only some rats may have been affected by the cross-fostering procedures.

Figure 5. Principal Component Analysis (PCA) represented by two-dimensional scatter plot (Figure 5(a,c), and hierarchical cluster analysis (Figure 5(b,d) using the DEG (FDR p ≤ 0.05; FC ≥ |1.5|). The DEG lists were derived from the comparisons of the treatment groups CC vs C and MM vs M at PND21 (Figure 5A and B); and CC vs C and MM vs M at PND78-86 (Figure 5(c,d). Treatment groups were manually circled to facilitate identification and clustering in the PCA plots. In the hierarchical cluster analyses the rows represent the probes for genes with increases (red) and decreases (green) in expression (black unchanged relative to controls).

At PND78-86 (Figure 5(c,d), the PCA and HCA noted some effects of cross-fostering. The PCA accounted for 44.8% (PC1 32.8% + PC2 12.0%) of the variation in the dataset. Here, PC1 contrasts MM group with the rest of the treatment groups, whereas PC2 appears to be a contrast of CC and MM groups to C and M treatment groups, respectively, with the exception of one of the MM samples. Both the PCA (Figure 5(c)) and HCA (Figure 5(d)) show the CC and MM groups forming two clusters. These observations, based upon a small number of DEG, support subtle cross-fostering effects, and perhaps a delay in induction of a more homogeneous cross-fostering influence on the liver transcriptome when results at PND78-86 as compared to those at PND21.

Functional analyses of differentially expressed genes by ingenuity pathway analysis (IPA)

The following sections compare effects of the mixture M at PND21 and at PND78-86, as well as cross-fostering effects (CC vs C and MM vs M), through identification of enriched upstream regulators (UR), canonical pathways (CP) and network/top functions, as predicted by IPA.

Comparison of M vs C at PND21 and at PND78-86

DEG used to predict upstream regulators (URs with Z score ≥ |2| and p < 0.01) and canonical pathways (CP with [–log (p < 0.01)] >2) potentially altered in group M at PND21 are identified by “x” in Table 2. Altered activities from 8 UR are suggested in group M; 6 of these (AHR, NR1I2, NR1I3, NFE2L2, HNF4, and PPARGC-1α) are predicted to be activated, and two to be inhibited (IL-1α, and IL-1β). Analysis of the DEG revealed 22 altered CP including AHR signaling, PXR/RXR activation, LPS/IL-1 mediated inhibition of RXR function, NRF2-mediated oxidative stress response, and PPARα/RXRα activation pathways. Globally, these enriched UR and CP may be broadly categorized into 5 functions: xenobiotic metabolism, inflammation, oxidative stress, lipid/energy metabolism, and endocrine system.

Table 2. Significant canonical pathways (CP) and upstream regulators (UR in bold) identified from DEG (FDR p ≤ 0.05 and FC ≥ |1.5|) in PND21 M rat livers. Arrows indicate if the DEG were up (↑), or down (↓) regulated.

	CP	UR	CP	CP	UR	UR	UR	UR	CP	UR	CP	CP	UR	UR	CP
PND21 M vs C	Xenobiotic metabolism signaling	AHR	Aryl hydrocarbon receptor (AHR) signaling	PXR/RXR activation	NR1I2 (PXR)	NR1I3 (CAR)	IL-1α	IL-1β	LPS/IL-1 mediated inhibition of RXR function	NFE2L2 (NRF2)	NRF2-mediated oxidative stress response	Glutathione-mediated detoxification	HNF-4α	PPARGC-1α	PPARα/RXRα activation
A2m ↓		x						x
Alas1 ↑				x	x	x			x	x				x
Aldh1a1 ↑	x		x	x	x	x	x		x				x
Aldh1b1 ↓	x	x	x						x
Apoa2 ↑													x		x
Ces1c ↑					x
Ces2a ↑					x	x
Cyp17a1 ↓		x					x
Cyp1a1 ↑	x	x	x		x		x	x					x	x
Cyp1a2 ↑	x	x	x	x	x	x	x	x					x	x
Cyp2a1 ↑									x
Cyp2a3 ↑		x		x	x	x			x				x	x
Cyp2b2 ↑	x	x		x	x	x	x	x	x				x	x
Cyp2c7 ↑	x	x		x	x	x			x				x		x
Cyp3a23 ↑	x			x	x	x	x	x	x				x	x
Entpd5 ↑					x					x
Ephx1 ↑						x				x	x		x
Gsta3 ↑	x	x	x						x	x	x	x
Gsta4 ↑					x	x				x		x
Gsta5 ↑	x	x	x		x	x			x	x	x	x	x
Gstm2 ↑	x		x	x	x	x			x	x	x	x
Nat8 ↑					x							x	x
Nfil3 ↓								x
Nqo1 ↑	x	x	x					x		x	x
Nr0b2 ↑			x	x	x			x	x	x			x		x
Onecut1 ↓													x
Sult1b1 ↑	x								x
Ubd ↑								x
Ugt1a6 ↑	x	x			x					x			x
Ugt2a3 ↑													x
Ugt2b1 ↑	x					x	x
Ugt2b37 ↑	x												x
-log[p-value]*	16.40	9.19	9.35	12.70	21.92	17.13	5.34	3.50	14.50	6.52	4.11	8.59	3.65	3.84	2.28
Ratio (z-score)**	15/165	(2.78)	9/115	9/54	(2.92)	(3.39)	(-2.06)	(-2.42)	13/162 (-1.0)	(2.80)	5/131	5/24	(2.21)	(2.42)	3/127 (-1.0)

*The –log(p-value) from the Fisher’s exact test (cutoff <0.01) assessing the overlap of observed and predicted genes.

**Ratio indicates the number of DEG in a CP divided by the total number of genes in that CP. Z-scores are presented within brackets and evaluates the relationship of observed and predicted patterns (up/down-regulation), Z ≥2 are predicted to be activated whereas Z ≤-2 are potentially inhibited pathways (2 is the default cutoff value for IPA CP and UR).

Table 2. (Continued.)

	CP	CP	CP	CP	CP	CP	CP	CP	CP	CP	CP	CP	CP	CP	CP
PND21 M vs C	Fatty acid α-oxidation	Melatonin degradation I	Nicotine degradation III/II	Serotonin degradation	Bupropion degradation	Acetone degradation I	Estrogen biosynthesis	Thyroid hormone metabolism II	Dopamine degradation	Oxidative ethanol degradation III	Histamine degradation	Putrescine degradation III	Ethanol degradation IV/II	Tryptophan degradation X	Noradrenaline and adrenaline degradation
A2m ↓
Alas1 ↑
Aldh1a1 ↑	x			x					x	x	x	x	x	x	x
Aldh1b1 ↓	x			x					x	x	x	x	x	x	x
Apoa2 ↑
Ces1c ↑
Ces2a ↑
Cyp17a1 ↓
Cyp1a1 ↑		x	x		x	x	x
Cyp1a2 ↑		x	x		x	x	x
Cyp2a1 ↑
Cyp2a3 ↑		x	x		x	x	x
Cyp2b2 ↑		x	x		x	x	x
Cyp2c7 ↑		x	x		x	x	x
Cyp3a23 ↑		x	x		x	x	x
Entpd5 ↑
Ephx1 ↑
Gsta3 ↑
Gsta4 ↑
Gsta5 ↑
Gstm2 ↑
Nat8 ↑
Nfil3 ↓
Nqo1 ↑
Nr0b2 ↑
Onecut1 ↓
Sult1b1 ↑		x		x				x	x
Ubd ↑
Ugt1a6 ↑		x	x	x				x
Ugt2a3 ↑		x	x	x
Ugt2b1 ↑		x	x	x				x
Ugt2b37 ↑		x	x	x				x
-log[p-value]*	4.59	18.4	16.5	10.9	10.8	10.8	9.63	4.35	5.84	4.88	4.73	4.47	4.47	4.35	4.47
Ratio (z-score)**	2/11	11/35	10/34	7/32	6/19	6/19	6/27	2/13	3/16	2/9	2/10	2/12	2/12	2/13	2/12

*The –log(p-value) from the Fisher’s exact test (cutoff <0.01) assessing the overlap of observed and predicted genes.

**Ratio indicates the number of DEG in a CP divided by the total number of genes in that CP. Z-scores are presented within brackets and evaluates the relationship of observed and predicted patterns (up/down-regulation), Z ≥2 are predicted to be activated whereas Z ≤-2 are potentially inhibited pathways (2 is the default cutoff value for IPA CP and UR).

Table 3 summarizes top functions identified by gene networks generated with IPA based on the relationships among DEG derived from the comparisons of M vs C groups at PND21 and at PND78-86. In M treated rats at PND21, drug metabolism and glutathione depletion were identified as top hepatic functions potentially affected. These are functions similar to xenobiotic metabolism and oxidative stress previously identified in Table 2. In addition, lipid metabolism and small molecule biochemistry were predicted to be affected in both age groups. Other top functions not detected at PND21, but identified in the M group at PND78-86, included molecular transport, cell cycle, cancer, and gastrointestinal disease. Moreover, bile acid biosynthesis was a CP identified in the M group at PND78-86 (Table 3) but not at PND21 (Table 2). There was no significant UR identified for the comparison M vs C at PND78-86. Overall, these observations suggest transient induction of xenobiotic/drug metabolism and oxidative stress, but persistent metabolic disturbances evident through different DEG at both time points, all of which may be conducive to disease development.

Table 3. Top functions identified by gene networks, and canonical pathway (CP), associated with DEG (FDR p ≤ 0.05 and FC ≥ |1.5|) in the liver of M treated rat at PND21 and PND78-86. Arrow represents up (↑) or down (↓) regulation

PND 21 M vs C								PND 78-86 M vs C
Top functions								Top functions								Canonical Pathway
DEG in Network 1	Drug metabolism	Glutathione depletion in liver	Protein synthesis	DEG in Network 2	Drug metabolism	Lipid metabolism	Small molecule biochemistry	DEGs in Network 1	Lipid metabolism	Molecular transport	Small molecule biochemistry	DEG in Network 2	Cell cycle	Cancer	Gastrointestinal diseases	DEG in Pathway	Bile acid biosynthesis
A2m ↓				Apoa2 ↑		x	x	Acsm3 ↓	x		x	C6 ↑		x	x	Akr1c1 ↓	x
Alas1 ↑				Cyb5a ↑	x	x	x	Akr1c1 ↓	x	x	x	Cav1 ↑	x	x	x	Akr1c2 ↓	x
Aldh1a1 ↑	x			Cyp17a1 ↓	x	x	x	Akr1c2 ↓	x	x	x	Eif4ebp3 ↓		x
Aldh1b1 ↓				Cyp1a1 ↑	x	x	x	Cyp4a2 ↑	x	x	x	Fabp4 ↑		x
Ces1c ↑				Cyp1a2 ↑	x	x	x	Dcn ↑				Fami11A ↓		x	x
Ces2a ↑				Cyp2a3 ↑	x	x	x	Esm1 ↑				Frzb ↑		x	x
Cyp2b2 ↑				Cyp2b2 ↑	x	x	x	F9 ↓				Gfra2 ↓		x	x
Cyp2c7 ↑	x			Cyp3a23 ↑	x	x	x	Fgf21 ↑	x	x	x	Gpm6a ↑		x	x
Entpd5 ↑			x	Ephx1 ↑	x	x	x	Igfbp6 ↑				Hspb1 ↓	x	x	x
Gsta3 ↑	x	x		Maged2 ↑				Inmt ↓				Kcnk2 ↑		x	x
Gsta4 ↑	x	x	x	Nqo1 ↑	x	x	x	Rgs4 ↑	x	x	x	Tpd52l1 ↑	x	x	x
Gsta5 ↑	x	x	x	Nr0b2 ↑		x	x	Rt1-Bb ↑				Usp9x ↑		x	x
Gstm2 ↑		x	x	Onecut1 ↓			x	Serpinb7 ↓
Nat8 ↑	x		x	Selenbp1 ↑
Nfil3 ↓			x	Slc34a2 ↓		x
Sult1b1↑	x			Ugt1a6 ↑	x		x
Ubd ↑				Ugt2a3 ↑
				Ugt2b1 ↑	x	x	x
				Ugt2b37 ↑	x	x	x
-log[p-value]*		40				37				24				22		-log[p-value]*	2.96
Focus genes		17				19				13				12		Ratio**	2/19

* The –log(p-value) from Fisher’s exact test that assesses the overlap of observed and predicted genes (default cutoff <0.01 for CP and <0.001 for gene network)

** The ratio indicates the number of DEG in a CP divided by the total number of genes in that CP.

Influence of cross-fostering at PND21 and PND78-86

Table 4 summarizes IPA predictions for effects of cross-fostering in untreated rats (CC vs C) at PND21 and at PND78-86. The analysis at PND21 revealed 11 enriched CP and one UR (LEP - predicted to be inhibited) potentially affected by cross-fostering. Because of the few DEG at PND78-86, there was no significant UR and only one CP (TR/RXR activation) associated with cross-fostering. Energy production was enriched only at PND21, and molecular transport was enhanced only at PND78-86. Interestingly, lipid metabolism and small molecule biochemistry are top functions predicted to be perturbed in both age groups but by different DEG (except for Ubd down-regulation in both age groups).

Table 4. Top functions identified by gene networks, canonical pathway (CP), and upstream regulators (UR, bold), respectively, associated with DEG (FDR p ≤ 0.05 and FC ≥ |1.5|) in the livers of CC rats at PND21 and PND78-86. Arrows represent up (↑) or down (↓) regulation

PND21 CC vs C																PND78-86 CC vs C
	Top functions			CP	UR	CP	CP	CP	CP	CP	CP	CP	CP	CP	CP		Top functions			CP
DEGs in Network/Pathways	Lipid metabolism	Small molecule biochemistry	Energy production	Calcium signaling	LEP	LPS/IL-1 mediated Inhibition of RXR Function	PXR/RXR activation	Trehalose degradation II	GDP-glucose biosynthesis	Glucose and glucose-1-phosphate degradation	Intrinsic prothrombin activation pathway	Coagulation system	Melatonin degradation I	Geranylgeranyldiphosphate biosynthesis	Trans, trans-farnesyl diphosphate biosynthesis	DEGs in Network/ Pathways	Lipid metabolism	Small molecule biochemistry	Molecular transport	TR/RXR activation
Abcb11↓	x	x			x	x	x									Acot5↓	x	x
Ccl21↑																Akr1c12↓
Cpt1b↑	x	x	x		x	x										Akr1c2↓	x	x	x	x
Cyp2c7↓	x	x	x		x	x	x						x			Apoa5↓	x	x		x
Des↑		x														B3galt1↑	x	x	x
Fabp3↑	x	x	x		x	x										Dbp↑	x	x
Fdps↑	x	x												x	x	Sult2al1↑	x	x
Fgb↓											x	x				Ubd↓
Fhl1↑
Gck↓	x	x	x		x			x	x	x
Igtp↓		x
Klkb1↓	x	x									x	x
Myl4↑				x
Pla2g1b↑	x	x	x
Ppp1r1a↑
Sult1b1↑						x							x
Tnni2↑				x
Tpm1↑				x
Ubd↓
-log[p-value]*	42			3.01	3.03	4.18	2.01	2.56	2.3	2.08	2.82	2.51	2.31	2.25	2.07		19			2.61
Ratio** (Focus genes) [Z-score]***	(18)			3/72	[-2.18]	5/162	2/54	1/1	1/2	1/3	2/21	2/30	2/38	1/2	1/3		(8)			2/64

* The –log(p-value) from Fisher’s exact test that assesses the overlap of observed and predicted genes (default cutoff <0.01 for CP and <0.001 for gene network)

** The ratio indicates the number of DEG in a CP divided by the total number of genes in that CP.

*** Z-scores are presented within brackets and evaluate the relationships of observed and predicted patterns (up/down-regulation), Z ≥ 2 are activated, Z ≤ -2 are inhibited pathways (2 is the default cutoff value for IPA CP and UR).

Table 5 describes the results of IPA for MM vs M comparisons in both age groups. Lipid metabolism, small molecule biochemistry, and vitamin and mineral metabolism were identified as enriched top functions potentially modulated by cross-fostering in both age groups. These top functions were identified based upon DEG that differed with age. Five enriched CP were affected by cross-fostering at PND21 based upon three DEG (Akr1c1, Cyp2c7, Cyp2d3; Table 5). The comparison of the expressions of these three genes and of Ubd in the MM vs M and M vs C groups suggests a pattern for cross-fostering effects at PND21. Notably, two of these DEG were not affected in the M vs C comparison (Akr1c1 and Cyp2d3), while the third one (Cyp2c7) was increased by M (Table 2) but down-regulated at PND21 in the MM compared to M group. Moreover, Ubd was upregulated in the M group (M vs C, Table 2), but downregulated in MM vs M comparison (Table 5). Therefore, these 4 examples at PND21 suggest that cross-fostering reduced the magnitude of changes in gene expression induced by M. No significant UR was identified at PND21, and no significant CP or UR were identified at PND78-86.

Table 5. Top functions identified by gene networks and canonical pathway (CP), associated with DEG (FDR p ≤ 0.05 and FC ≥ |1.5|) in the livers of cross-fostered M treated rats at PND21 and PND78-86. Arrows represent up (↑) or down (↓) regulation.

PND21 MM vs M									PND78-86 MM vs M
	Top functions			CP	CP	CP	CP	CP		Top functions
DEGs in Networks/ Pathways	Lipid metabolism	Small molecule biochemistry	Vitamin and mineral metabolism	Estrogen biosynthesis	Bupropion degradation	Acetone degradation I	Nicotine degradation III/II	Melatonin degradation I	DEGs in Networks/ Pathways	Lipid metabolism	Small molecule biochemistry	Vitamin and mineral metabolism
Akr1c1↑	x	x	x	x					Akr1c12↑
Akr1c12↑									Akr1c2↑	x	x	x
Arntl↓	x	x							Aldh1a7↓		x
Cyp2c7↓	x	x	x	x	x	x	x	x	Cyp2b2↓	x	x	x
Cyp2d3↓	x	x	x	x	x	x	x	x	Cyp2c12↑
Dbp↑	x	x							Fgl1↓	x	x	x
Lcn2↓	x	x							Inhbc↓
Retsat↓	x	x	x						Pcolce↓
Tff3↓	x	x							Srd5a1↓	x	x	x
Tpt1↓		x							Tfr2↓		x
Ubd↓
-log[p-value]*		27		4.87	3.28	3.28	2.8	2.8			22
Ratio** (Focus genes)		(11)		3/31	2/22	2/22	2/37	2/38			(9)

* The –log(p-value) from Fisher’s exact test that assesses the overlap of observed and predicted genes (default cutoff <0.01 for CP and <0.001 for gene network)

** The ratio indicates the number of DEG in a CP, divided by the total number of genes in that CP.

Overall, cross-fostering effects in control (Table 4, CC vs C) and M treated rats (Table 5, MM vs M) were detected, but considering the small number of DEG and fold changes (Table S2), these effects appear to be subtle.

Discussion

The current study demonstrated that cross-fostering induces differences in gene expression between control groups (CC vs C) and modulates effects of developmental exposure to the EC mixture on hepatic transcriptomes at PND21 and at PND78-86. To our knowledge, this is the first investigation of influence of cross-fostering on the hepatic transcriptome. These data are in addition to previously reported observations from the same experiment, which indicated that cross-fostering amplified the adverse effects of exposure to the mixture on the glucocorticoid-stress response in adult rats (Desaulniers, Xiao, and Cummings-Lorbetskie 2013). Collectively, these findings suggest an interaction between treatment and cross-fostering, and warrants careful interpretation of other cross-fostering studies in the risk assessment process. Importantly, this investigation shows that non-chemical stressors may modulate chemical-mediated effects.

The current PND21 hepatic transcriptome data highlight two other significant concepts. First, the dose-response pattern in number of DEG and in Cyp1a1 mRNA abundance reflected the TCDD-TEQ doses across groups. Consequently, these results support the use of the TCDD-TEQ system (Van Den Berg et al. 2006) in risk assessment procedures to assess biological potencies of AhR-agonists in complex mixtures. Second, in all comparisons the DEG differed between age groups. In general, the early gene sets at PND21 were associated with xenobiotic detoxification, oxidative stress, energy and retinoid metabolism, and immune responses, whereas at PND78-86 different gene sets were associated with retinoids and energy metabolism as potential long-term consequences. The chronic effects were subtle, perhaps given that the chemical residue levels in these rats declined to concentrations similar to those encountered in some human populations (Desaulniers et al. 2017). Although species differences obscure interspecies extrapolation (Fader and Zacharewski 2017), the DEG at PND78-86 might assist in predicting mechanisms conducive to disorders previously described in EC exposed human populations (Andersen et al. 2014; Kuwatsuka et al. 2014; Mitoma et al. 2015; Warner et al. 2013), including dyslipidemia in Northern Canadian populations (Marushka et al. 2018; Singh and Chan 2018).

The following sections discuss the transcriptomic effects among groups C, CC, M, and MM, with supportive biochemical measurements previously published from the same rats (Desaulniers, Xiao, and Cummings-Lorbetskie 2013; 2017), or from studies testing other preparations of the same or similar mixtures (Chu et al. 2005, 2008; Elabbas et al. 2014).

Influence of developmental exposure to mixture M at PND21

Functional analysis of the list of DEG in IPA identified nuclear receptors (AHR, PXR, CAR, RXR, PPARGC-1α, NRF2, HNF4α) that are upstream regulators of pathways predicting induction of oxidative stress, xenobiotic/drug, hormone (thyroid, steroids), and lipid metabolisms. The IPA prediction of induction of oxidative stress is in agreement with studies observing elevated hepatic glutathione S-transferase (GST) activities in Sprague-Dawley rats exposed to TCDD (Baars, Jansen, and Breimer 1978) or to non-ortho or ortho-chlorinated PCB (Lamartiniere, Dieringer, and Lucier 1979; Twaroski, O’Brien, and Robertson 2001), which are present in the mixture. The current rats exhibited increased activities of the major liver enzymes involved in xenobiotic metabolism, including CYP1A, 2B, and 3A (Desaulniers et al. 2017). These enzymes alter steroid metabolism. It is of interest that developmental exposure to a different preparation of the same mixture favored the production of catechol estrogens (Desaulniers et al. 2012). Finally, the occurrence of altered lipid and retinoid metabolism is well-established. Developmental exposure studies with a different preparation of the current mixture reported that rats at PND35 exhibited elevated serum cholesterol levels at a dose of 5 mg/kg bw/day, but not at 0.5 mg/kg bw/day (Chu et al. 2008), and a dose of 0.05 mg/kg bw/day decreased abundance of hepatic retinol (Elabbas et al. 2014). The latter effect on retinoids occurs at the lowest observable effect level (LOEL) induced by this or similar mixtures. The up-regulation of Cyp2c7 [alias Cyp2c39 (Andreola et al. 2004; Fan et al. 2004)], and Ces (Mello et al. 2008) support the presence of altered retinoid (vitamin-A and its metabolites) metabolism.

IPA suggests inhibition of the upstream regulators interleukin-1α (IL-1α) and IL-1β, supporting alterations of the immune response pathways. Other investigators found that IL-6, IL-1α, and IL-1β down-regulation was induced by PXR activation in mouse liver (Moreau et al. 2007; Zhou et al. 2006), or by AhR-agonist exposure (Li et al. 2017; Stockinger et al. 2014), and immunosuppression by cis- and trans-nonachlors (two compounds present in the mixture M) exposure in rodents (Tryphonas et al. 2003). Chu et al. (2005, 2008), using a mixture and exposure protocol similar to the current experiment but at a higher dose (1.3 or 5 mg/kg bw/day), observed mild inflammation in histological sections from livers of rat offspring. In the current study, some DEG at PND21 (down-regulated Nfil3, α2M, Mmd2, Defa10, Ratnp-3b, and up-regulated Ceacam10) indicate impairment of innate immunity and of hepatic acute-phase protein system (Suppl. section S3.4). Collectively, these findings suggest that the immune system is less efficient in EC-exposed animals. It is interesting to note that a high rate of infectious disease was observed in Northern populations exposed to persistent pollutants (Kazbariene, Krikstaponiene, and Monceviciute-Eringiene 2006).

Finally, there are numerous markers of preneoplastic foci measured at PND21 including up-regulation of Gsts, α2M, Ubd (alias FAT10) (French 2010; Kimura et al. 2016; Ren et al. 2006), Mage-d2 (Hashimoto et al. 2015), and Entpd5 (Vogiatzi et al. 2016), and down-regulation of Slc34a2 (Li, Chen, and Lu 2016; Yang et al. 2014). The fact that chronic exposures to many of the chemicals present in the mixture are known to induce hepatocarcinomas in rodents (Brown et al. 2007), raises interest in the possibility that such potential preneoplastic markers might be detectable as early as PND21. However, the detection of neoplastic markers is not necessarily an indication of cancer development given that neoplastic sites might be transient.

Effects of developmental exposure to mixture M at PND78-86

In contrast to rats at PND21, mRNA levels of detoxification enzymes at PND78-86 were not significantly different from control, which is in agreement with the low chemical residue levels and absence of induced p450 activities in these rats (Desaulniers et al. 2017). Nevertheless, DEG at PND78-86 suggest effects on energy (lipid/glucose) and retinoid metabolism, molecular transport, small molecule biochemistry, bile acid biosynthesis, pro-inflammatory mechanisms, cell cycle, and contribution to diseases (Tables 2 and 3).

DEG at PND78-86 indicated altered retinoid metabolism and bile acid biosynthesis. Akr1c1 (AKR1C4 in human), Akr1c2 (AKR1C1 in human) and murine-specific Akr1c12, were down-regulated (indicative of cholesterol and bile acid metabolism deregulation), while Sdr16c6 was up-regulated in the M group, potentially contributing to altered retinoid metabolism. These aldo-keto reductase transcripts (AKR families) are phase I metabolizing enzymes that catalyze conversion of aldehydes and ketones to their corresponding alcohols, and thus reduce a variety of lipophilic substrates including xenobiotics, steroids, prostaglandins, bile acids, and retinoids (Rizner and Penning 2014; Ruiz et al. 2012). In the opposite enzymatic direction, the short-chain dehydrogenase/reductase SDR16C family members (Sdr16c6 up-regulated here) are thought to oxidize retinol to retinaldehyde (Adams et al. 2017). Depletion of vitamin-A (fat soluble retinoids) is one of the most sensitive markers of exposure to TCDD (Fletcher, Hanberg, and Håkansson 2001). Morse and Brouwer (1995) noted effects on the retinoid system until at least PND90 after gestational exposure to the highly toxic PCB mixture Aroclor 1254 at doses of 5 or 25 mg/kg body weight. Using a mixture similar to ours, Elabbas et al. (2014) observed decreased abundance of hepatic retinol in the livers of male offspring at PND35, but increased abundance at intermediate doses (0.05 and 0.5 mg/kg/day) relative to control in both F1 males and females at PND77. The persistent deregulation of retinoid metabolism and its health impacts warrants further investigation.

DEG at PND78-86 suggested altered energy metabolism, steatosis, and inflammation, which may be conducive to diseases. The relationships among the up-regulated DEG Fabp4, Cav1, Cyp4a2, and Fgf21, are indicative of activation of the PPAR pathway (Su et al. 2005), and long-term altered lipid and glucose metabolism. In addition, down-regulation of Acsm3, which activates fatty acids by forming a thioester with coenzyme A (Watkins et al. 2007), and down-regulation of Cidea coding for a lipid droplet component (Xu et al. 2016), are evidence of altered energy metabolism. Lrp5, Frzb, Ctnnal1, and Slc35c2 are DEG that may further support effects on lipid metabolism and/or early carcinogenic events (Suppl. section S4.1). LRP5 plays a role in lipolysis as an endocytotic receptor mediating the clearance of apoB- and apoE-containing lipoproteins (Julve et al. 2016). The up-regulation of Frzb and Lrp5, and down-regulation of Ctnnal1, support changes in the Wnt/β-catenin canonical signaling pathway which regulate the expression of genes involved in cholesterol and glucose metabolism (Fujino et al. 2003), liver injury, and tumorigenesis (Fan et al. 2011; Merdek, Nguyen, and Toksoz 2004; Monga 2015; Zhang et al. 2016). Slc35c2 (solute carrier family 35 member C2) when up-regulated (as noted here) was found to compete with Slc35c1, and such alterations were linked to cancer (Lu et al. 2010; Mehta et al. 2016). Mixture-induced inflammatory events may be attributed to interactions among the DEG Fabp4 [a multiple role hepatic gene (Hu et al. 2017; Leroux et al. 2012)], with C6, and Casp12 (Jimenez and Lamkanfi 2015; Sarma and Ward 2011) (Suppl. section S4.1). Collectively, DEG noted at PND78-86 might provide an indication of intermediary steps linked to steatohepatitis and/or early carcinogenic events. This is supported by observations of enhanced hepatic inflammation and steatosis associated with cholesterol accumulation and alterations in fatty acid profile in rats treated with similar mixtures of contaminants (Chu et al. 2005, 2008; Mailloux et al. 2014), and by increases in liver cancer incidence in rats chronically exposed to industrial PCB mixtures (Brown et al. 2007).

Influence of cross-fostering (CC vs C and MM vs M)

Data suggested subtle cross-fostering effects in both age groups on energy metabolism (lipids, glucose) in control groups (CC vs C), and a modulation of energy, retinoid, and steroid (glucocorticoid) metabolisms among mixture treated rats (MM vs M).

DEG related to energy (glucose, lipid) and retinoid metabolism in CC vs C at PND21

DEG that are down-regulated (Gck, Abcb11, and Cyp2c7) and up-regulated (Cpt1b, Fabp3, Fdps) in CC compared to C group suggest altered energy (glucose/lipid) metabolism. Gck and Cpt1b are major energy metabolism genes. Gck codes for glucokinase, which catalyzes a rate-limiting step for the entry of glucose into the glycolytic pathway producing pyruvate to fuel the mitochondrial TCA cycle, while carnitine palmitoyltransferase 1B (Cpt1b, up-regulated here) catalyzes transport of long-chain fatty acyl-CoAs from the cytoplasm to carnitine for translocation across the mitochondrial inner membrane (McGarry and Brown 1997). The other DEG involve altered lipid metabolism (up-regulation of Fabp3, down-regulation of Cyp2C7), cholesterol (the precursor of bile salts) and bile salt metabolism (down-regulation of Abcb11 and up-regulation of Fdps) (Suppl. section S5.1). These energy metabolism transcriptomic anomalies might be linked to the tendency to produce heavier males in cross-fostered groups at PND21 (mean body weights were: C 55.1 g, CC 57.4 g, M 53.1 g, MM 55.8 g, p = 0.06 (Desaulniers, Xiao, and Cummings-Lorbetskie 2013)).

DEG supporting altered lipid metabolism in CC vs C at PND78-86

At PND78-86, the transcriptomic data support the hypothesis that cross-fostering alters lipid metabolism, through up-regulation of Dbp and Slc13a5 and down-regulation of Akr1c2, Akr1c12, Acot5, and Apoa5 in the CC group. First, Dbp (a gene influenced by circadian rhythm also detected in the MM vs M comparison, see next section) up-regulates expression of PPARγ which is an important regulator of adipocyte differentiation and lipid metabolism (Takahashi et al. 2010). In addition, up-regulated Slc13a5 (solute carrier family 13 member 5) is known to be involved in fat and glucose metabolism (Pesta et al. 2015; Willmes, Helfand, and Birkenfeld 2016). As previously discussed, aldo-keto reductases (AKR1C2, AKR1C12) reduce a variety of lipophilic substrates (steroids, prostaglandins, bile acids, and retinoids (Rizner and Penning 2014; Ruiz et al. 2012)), whereas acyl-CoA thioesterase 5 (ACOT5) catalyzes the hydrolysis of acyl-CoA esters to their corresponding non-esterified fatty acid and coenzyme A (Brocker et al. 2010). Finally, apolipoprotein A5 (APOA5) is a component of several lipoproteins (VLDL, HDL, chylomicrons) and acts as a potent regulator of triglyceride metabolism (Julve et al. 2016). Overall, the DEG observed in the CC group at PND78-86 are indicators of changes in lipid and energy metabolisms. (Suppl. section S5.2).

Hydroxysteroid sulfotransferase, Sult2a1, up-regulated in CC, plays an important role in xenobiotic and hormone metabolism. Maiti et al. (2004) in an in vivo rodent study suggested that physical and/or chemical stress imposed on 5 to 7 week old SD male rats resulted in increased hepatic Sult2a1 expression. Further, Duanmu et al. (2002) reported induction of human and murine SULT2A1 enzyme by glucocorticoids. These findings of up-regulation of Sult2a1 associated with stress and glucocorticoid hormone are consistent with those of this cross-fostering experiment.

DEG related to retinoid and energy (lipid, glucose) metabolism in MM vs M at PND21

The down-regulation of Cyp2c7, Retsat, Lcn2, and Tff3, and up-regulation of Akr1c1 in the MM group relative to M suggests that cross-fostering modulated the influence of M on mechanisms involving retinoid and lipid metabolism. As discussed previously, CYP2C7 is a retinoic acid and steroid hydroxylase enzyme (Fan et al. 2004), whereas AKR1C1 (AKR1C4 in human) is a member of the aldo-keto reductase family that catalyzes the reduction of retinaldehyde to retinol (Ruiz et al. 2011). Retinol saturase (Retsat) is an oxidoreductase that is also involved in retinoid metabolism, catalyzing conversion of retinol to 13,14-dihydroretinol (Moise et al. 2004). Lipocalin-2 (LCN2) is a transport protein for lipids, steroid hormones, and retinoids. Recently Guo et al. (2016) reported that Lcn2 deficient mice displayed diminished levels of retinol, retinyl ester, and all-trans retinoic acid in adipose tissue. As previously discussed, a different preparation of the same mixture produced a reduction in hepatic retinol (Elabbas et al. 2014). Taken together, these results suggest that cross-fostering modulated the effect of M on retinoid metabolism; however, retinoid biochemistry analyses would be required to determine if cross-fostering amplified or lowered the M-induced retinoid deficiency. Implication of Retsat and Tff3 on lipid and glucose metabolism is presented in Suppl. S5.3.

Interestingly, cross-fostering might deregulate hepatic circadian rhythm and thus lipid metabolism, due to down-regulation of Arntl (alias Bmal1 or Mop3) and up-regulation of Dbp in the MM vs M group. Indeed, Dbp governs the circadian transcription of hepatic metabolic and detoxification enzymes in preparation for food intake (Schmutz, Albrecht, and Ripperger 2012; Stratmann et al. 2010; Takahashi et al. 2010), whereas ARNTL, as a heterodimer ARNTL/CLOCK, regulates the circadian expression of genes (Jaeger and Tischkau 2016; Leliavski et al. 2014; Zvonic et al. 2006). Thus, potential effects of cross-fostering on the CLOCK gene network may exert adverse metabolic consequences.

DEG in MM vs M at PND78-86

The down-regulation of Srd5a1, Cyp2b2, Inhbc, Tfr2, and Fgl1, and up-regulation of Akr1c2 and Akr1c12, in the MM relative to the M group at PND78-86 support that cross-fostering modulated the effects of M on steroid and bile acid metabolisms. The steroid-5α-reductases type 1 gene (Srd5a1) catalyzes A-ring reduction of pregnene steroids including androgens and glucocorticoids. Rodents deficient in SRD5A1 (knockout mice, or treated with enzyme inhibitors) accumulate hepatic excesses of glucocorticoids and lipids, develop hepatic steatosis, insulin resistance, and features of metabolic syndrome (Livingstone et al. 2015). Perhaps the down-regulation of Srd5a1 contributed to the deregulated glucocorticoid stress-responses in rats of the MM group (Desaulniers, Xiao, and Cummings-Lorbetskie 2013). Cyp2b2 (CYP2B6 in human) is involved in drug metabolism and in the synthesis of steroids, cholesterol, and other lipids. Chen and co-workers (2014) reported that expression of CYP2B6 was decreased in diseased liver, including progressive cirrhosis and hepatocellular carcinoma. The down-regulation of Cyp2b2 in MM relative to M suggests that liver functions may be more severely affected in MM than M group. This finding is strengthened by observations that Inhbc, Tfr2, and Fgl1 were not markedly altered by any other treatment, and their down-regulation in MM relative to M is consistent with impaired liver regenerative response (Esquela et al. 1997), iron metabolism (Delima et al. 2012), and tumor suppressor functions (Nayeb-Hashemi et al. 2015), respectively. Effects on these genes and on Akr1c2 and Akr1c12 are further discussed in the Suppl. section S5.4.

Cross-fostering as a form of non-chemical stressor

Overall, these results support an influence of cross-fostering in control rats and a modulation of the effects of M on hepatic metabolism genes (glucose, lipids, bile acid, retinoids, glucocorticoids). The mechanisms by which cross-fostering affected the hepatic transcriptome are unknown but may include altered transient nursing quality (post-natal maternal behavior, milk quality) with interactions with mixture effects. Other experiments exploring stress effects also noted altered energy metabolism. It is known that psychological stress induces dyslipidemia and hepatic steatosis in mice (Depke et al. 2008), and restraint-stress modifies glucose metabolism (Cohen et al. 2009; Nirupama, Devaki, and Yajurvedi 2010), including hepatic Gck down-regulation (Tsoi et al. 2011). Interestingly, in non-human primates, stress imposed on nursing mothers during early development led to peripubertal obesity and insulin resistance in juveniles (Kaufman et al. 2007). Globally, the current transcriptomic data support the hypothesis that cross-fostering is a form of non-chemical stressor that induces hepatic metabolic disturbance in the developing liver. This is associated with subtle transcriptomic differences in control and mixture treated rats of young and adult ages, and might promote intermediary changes contributing to disease development.

Conclusions

This study in SD male rats showed that developmental exposure to EC mixture with a composition relevant to human exposure, induced changes in hepatic gene expression at PND21 that are indicative of effects on xenobiotic metabolism, oxidative stress, and energy (lipid, glucose) metabolism. Remarkably, at PND78-86, a different set of DEG was altered which was less suggestive of xenobiotic metabolism but still associated with altered energy metabolism. Thus, developmental exposure to EC mixtures may lead to short- and long-term effects mediated by different genes, but with the potential to alter the same set of functions related to energy metabolism.

Our data highlight some challenges for the risk assessment of environmentally relevant mixtures to humans. First, timing of tissue sampling might generate different lists of DEG and influence interpretation of potential outcomes. Indeed, changes in gene expression are dynamic (PND21 vs PND78-86), and at a given time and/or dose, DEG may reflect homeostasis, tissue adaptation, or toxicity. In that context, identification of consequential pathway(s) is difficult and might be confounded among more trivial pathways indirectly activated through cross-talk (Denison et al. 2011; Miao et al. 2006; Rezen et al. 2011; Tian 2009). To better understand the risk of developing a later adverse outcome, and identify DEG that are truly indicative of intermediary steps toward disease development, more low dose studies are needed with samples collected at different times following exposure to reflect immediate response, adaptation, compromised function and disease manifestation. Secondly, the identification of “driver chemicals” is another challenge when investigating the effects of mixtures. In the current study, the EC mixture was complex and included various organochlorines with different affinities for nuclear receptors [e.g.: 2,3,7,8-TCDD and co-planar PCBs with high affinity for AHR, but low affinity for CAR; or PCB153 and non-co-planar ortho-chlorinated PCBs with low affinity for AHR, but high affinity for CAR (Al-Salman and Plant 2012; Safe 1994)]. Based upon the dose-response observed on CYP expression and activities, the TCDD-TEQ system was reliable in identifying AhR-agonists as important “drivers” of short-term effects.

Current data, which suggests that cross-fostering acts as a form of non-chemical stress, strengthens previous observations that cross-fostering interacted with exposure to EC in inducing abnormal glucocorticoid stress-response at PND78-86 (Desaulniers, Xiao, and Cummings-Lorbetskie 2013). Humans are exposed to non-chemical stressors of various types that may modulate the effects of EC. Human exposures to EC (e.g.: PCDD, PCDF, PCB, pesticides) are associated with development of metabolic disorders such as obesity, diabetes and metabolic syndrome (Li et al. 2015; Marushka et al. 2018; Mostafalou and Abdollahi 2013; Singh and Chan 2018; Vafeiadi et al. 2016; Warner et al. 2013). Overall, transcriptomic clues are provided as to the origin of EC-induced metabolic dysfunction to guide future investigations. Our findings promote the need to identify and quantify the sources of non-chemical stressors that modulate the effects of chemical exposure. Finally, data suggest that consideration of non-chemical stressors is warranted to improve chemical risk assessments in humans.

Acknowledgments

The authors are grateful to Dr. Dalibor Breznan, Dr. Craig Parfett and Dr. Guillaume Pelletier for carefully reviewing this manuscript, to Dr. Mainul Husain and Andrea Rowan-Carroll for microarray hybridization, to Ming-Yan Yang, and the personnel from the Animal Resources Division for their precious assistance. Funded by Health Canada, the Northern Contaminants Program, and the Chemicals Management Plan of the Canadian government.

Supplemental data

Supplemental data for this article can be accessed here.